Quantization SNR of sine wave doesn't match 1.761 + 6.02 * Q

Question

I am trying to show with numpy that the quantization noise of a sine wave matches the SNR formula of SNR = 1.761 + 6.02 * Q.

The numpy code is simple:

import numpy as np
import matplotlib
from matplotlib import pylab, mlab, pyplot
plt = pyplot

from pylab import *
from numpy import *
from scipy import signal

def quantization_noise(quant):
    N=8192
    freq = 128
    x = np.linspace(0., 1., N)

    y1 = 0.5 * np.sin(2 * np.pi * freq * x)

    y2 = (np.floor(quant * (y1)) / quant)
    diff = y2 - y1

    freqs = fftfreq(N)
    x_mask = freqs >= 0

    Y1 = np.fft.fft(y1)
    Y2 = np.fft.fft(y2)

    Y1db = 20 * np.log10(np.abs(Y1) / N * 4)[x_mask]
    Y2db = 20 * np.log10(np.abs(Y2) / N * 4)[x_mask]

    plt.plot(freqs[x_mask], Y1db, 'bx', label = "input")
    plt.plot(freqs[x_mask], Y2db, 'r-', label = "output")
    #plt.plot(freqs[x_mask], Y1db, 'bx')
    plt.ylim([-140, 5])
    plt.xlim([0, 0.5])

    snr = np.amax(Y2db[3*int(freq):])
    print(snr)

    plt.plot([0.0, 0.5], [snr, snr], 'm-.', linewidth=1.0)
    plt.text(0.3, snr+4, "SNR=%4.1fdB" % snr)

    plt.grid(True)
    plt.legend(loc=1)

if True:
    plt.figure(figsize=(10,6))
    quantization_noise(8)

    tight_layout()
    plt.savefig("quantization_noise_8.png")

    plt.figure(figsize=(10,6))
    quantization_noise(16)

    tight_layout()
    plt.savefig("quantization_noise_16.png")

When I look at the results, I get an SNR of 27.4dB for 3 bit of quantization. The theory predicts 19.8db.

Similarly, for 4 bits of quantization, I get an SNR of 36.1dB: ~9dB more than for 3 bits, where you'd a delta of 6dB.

Ultimately, I want to show how with 16 bits A/D conversion, you'd end up with 98dB, but as the quantization level increases, the output spectrum gets closer and closer to the input spectrum, which is a continuous downward slope, which raises the question at which point something is considered noise instead of part of the signal.

I used applied a hanning window to better isolate the sidelobes of the main signal, which, for 3 bits of quantization makes the SNR go up from the earlier 27.4dB to 33.3dB:

I'm trying to figure out where my understanding is lacking.

How can I numerically show demonstrate the validity of the 1.761 + 6.02Q theory?

Tom

Answer 1

Some issues here:

1. Your SNR formula only applies to full scale sine waves, your sine wave has -6dB amplitude so your SNR will be 6 dB lower
2. The formula also implies rounding, not truncation, that's another 6 dB
3. You use a frequency that's a small integer divider of the sample rate, that means you are just repeating the same samples over and over again and don't get enough sample coverage to get a statistically meaningful result.
4. Your SNR analysis in the frequency domain is needlessly complicated and prone do imprecision and masking errors. Just do it directly in the time domain.

Here is how this would look in Matlab

%% quantization noise of a 16-bit sine wave
fr = 975.3; % something odd
n = 8192;
quant = 2^15;
% make the sine wave
y0 = sin(2*pi*(0:n-1)'/n*fr);
% quantize
yq = round(quant*y0)/quant;
% noise
yNoise = yq-y0;
% SNR
fprintf('SNR = %6.2fdB\n', 10*log10(mean(y0.^2)/mean(yNoise.^2)));

Technically, you would also have to work around the fact that the positive max amplitude of the sine might clip, but for large quantization this makes no meaningful difference.

EDIT

It may be good to review where the formula comes from and what it actually means. It all starts with the quantization noise. If we quantize and round than the quantitation noise is uniformly distributed between $[-0.5 \delta ,-0.5 \delta ]$, where $\delta$ is the quantization step. For truncation it would be uniformly distributed $[0,\delta ]$. For rounding, the resulting noise powers are $$ P_{round} = \frac{\delta ^2}{12}, P_{trunc} = \frac{\delta ^2}{3}$$

If $B$ is the number of bits, then for a signed signal, we simply have $\delta = 2^{B-1}$ and hence we get the noise level for 16 bit rounding as

$$P_{16} = \frac{2^{-30}}{12} \Rightarrow -101.1 dB $$

That's true for ANY signal that has a reasonable wide distribution of samples. Since a full scale sine wave has a power of $-3dB$ the resulting SNR would be $-98.1 dB$

Answer 2

I was doing quite a bit wrong, but the key thing that I was missing was the fact that the SNR needs to be calculated over the whole Nyquist spectrum instead of only looking at the peaks.

This article explains everything very well: Taking the Mystery out of the Infamous Formula, "SNR = 6.02N + 1.76dB," and Why You Should Care.

Another issue was that the sample rate is an integer multiple of the frequency of my test sine wave. The same article talks about how this is an issue also in the analog domain when evaluating the performance of a real ADC. It can be avoided by either changing the frequency a little bit or by adding a dither noise to the input.

Answer 3

It's a partial answer to an old question. I tried to reproduce formulae 1.761 + 6.02 * N looping on N. It gives me this nice result with N ranging 24 to 48

But with N ranging 2 to 48 result is

I don't understand why results are bad with low value of N

My code is

import numpy as np
import matplotlib.pyplot

class Quantization:
    def __init__(self, u_min=-1, u_max=1, N=8):
        self.u_min = u_min
        self.u_max = u_max
        self.N = N
        self.q = ( self.u_max - self.u_min) / 2 ** self.N

    def dac(self, s):
        n = np.round((s  - self.u_min ) / self.q - 2 ** (self.N-1))
        n = np.clip(n, -2 ** (self.N-1), 2 ** (self.N-1) - 1)
        return n
        
    def adc(self, n):
        s = self.q * (n + 2 ** (self.N-1) ) +  self.u_min
        return s

snr = [] 
N_inter = [*range(2,48)]
f = 20.2554
Fe =80000
M = 32768 * 16
t = np.arange(0, M) / Fe
rng = np.random.default_rng()    
for N in N_inter:
    qac = Quantization(N=N,u_min=-1, u_max=1)   
    # y = t
    # y = np.linspace(qac.u_min, qac.u_max, M)
    y = qac.u_max * np.sin(2* np.pi * f * t)
    y_r = qac.adc(qac.dac(y))
    snr.append(10 * np.log10(np.mean(y ** 2) / np.mean(( y - y_r) ** 2)))

fig1, ax1 = matplotlib.pyplot.subplots(nrows=1, ncols=1)
ax1.scatter(N_inter,snr, marker='o')
coeff = np.polyfit(N_inter,snr,1)
print(coeff)
ax1.plot(N_inter, np.polyval(coeff, N_inter))
ax1.plot(N_inter, np.polyval([np.log10(2)*20, 10*np.log10(3/2)], N_inter))
ax1.set_xlabel('N-adc converter (N)')
ax1.set_ylabel('noise (dB) ')
leg_fit = str(int(coeff[0]*100)/100) + ' *N'
if coeff[1] >= 0:
    leg_fit = leg_fit + '+' + str(int(coeff[1]*100)/100)
else:
    leg_fit = leg_fit + str(int(coeff[1]*100)/100)
leg_theory = 'Theory ' + str(int(np.log10(2)*20*100)/100) + ' *N+' + str(int(10*np.log10(3/2)*100)/100)
ax1.legend(['SNR(N)', leg_fit, leg_theory])
ax1.grid()
matplotlib.pyplot.show()