Skip to main content
11 of 35
added 272 characters in body
Olli Niemitalo
  • 13.7k
  • 1
  • 35
  • 63

Theory

In the limit $f/f_s \to 0,$ where $f$ is the sine wave frequency and $f_s$ is the sampling frequency, the answer to the question will be found in the Fourier series of the continuous-time piece-wise constant quantized waveform. We can construct the waveform as a sum of components like the red curve here illustrated for sine wave amplitude $N = 7$:

enter image description here

We are working with cosine instead of a sine because the math is nicer this way. Because of the symmetries, each component has only odd harmonics and only cosine terms in its Fourier series. The non-zero coefficients of the Fourier series are given by:

$$a_n = \frac{4}{\pi}\int_{x_0}^{x_1}\cos(nx)\,k\,dx = \frac{4k}{\pi n}\big(\sin(x_1n) - \sin(x_0n)\big),\quad n\text{ odd},\tag{1}$$

with integer amplitude $k \in 1\ldots N$ in range $x_0 < x < x_1$ in the first quarter-period of the cosine. I did not bother to explicitly write the contributions from the symmetrically arranged pieces in the other quarters, because they will contribute identically to the odd harmonics. Instead I simply included the implicit factor 4 in the equation.

The Fourier series of the full piece-wise constant waveform is a sum of the Fourier series of the components. The boundaries of the pieces that we need to include in the sum are:

$$\begin{gather}0 < x < \operatorname{acos}\left(\frac{N - 0.5}{N}\right),\quad \text{if }k = N,\\ \operatorname{acos}\left(\frac{k + 0.5}{N}\right) < x < \operatorname{acos}\left(\frac{k - 0.5}{N}\right),\quad \text{if }k \in 0\ldots N-1.\end{gather}\tag{2}$$

The non-zero coefficients of the Fourier series of the full piece-wise constant waveform are then given by:

$$\begin{gather}a_n = \frac{4N}{\pi n}\sin\bigg(\operatorname{acos}\left(\tfrac{N - 0.5}{N}\right)n\bigg) +\\\frac{4}{\pi n}\sum_{k=1}^{N-1} \Bigg(k\sin\bigg(\operatorname{acos}\left(\tfrac{k - 0.5}{N}\right)n\bigg) - k\sin\bigg(\operatorname{acos}\left(\tfrac{k + 0.5}{N}\right)n\bigg)\Bigg),\quad n\text{ odd}.\end{gather}\tag{3}$$

16-bit quantization

It is not out of reach to compute Eq. 3 for $N=2^{23} - 1$ and $n \in {5, 7},$ as parameterized in the question, using something like Python's mpmath. Let's try that but with a 16-bit sine wave first:

import mpmath as mp
def a_n(N, n):
  return 4*N/(mp.pi*n)*mp.sin(mp.acos((N-0.5)/N)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/N)*n) - mp.sin(mp.acos((k+0.5)/N)*n) ), [1, N-1])

N = 2**15-1 # amplitude (integer)
n = 5  # harmonic number (odd integer)

mp.mp.prec = 53  # default precision
20*mp.log10(mp.fabs(a_n(N, n))/N) # Fourier series term amplitude in dB

mp.mp.prec = 106  # Compute again with double the default precision, see if we get the same result
20*mp.log10(mp.fabs(a_n(N, n))/N) # Fourier series term amplitude in dB

mp.mp.prec = 53  # Restore default precision

I'm only reporting as many digits as are agreed about by repeated computation using two different mpmath precision settings, including also the digits that got changed by $\pm1$ due to larger changes in later digits. For a 16-bit full-scale $N = 2^{15}-1$ sine wave and $f/f_s \to 0$, the result is -154.66449 dBFS for the 5th harmonic and, after setting n = 7 in the script, -154.66506 dBFS for the 7th harmonic, in about 15 seconds of computation on my PC for a single result with the higher precision setting. As a sanity check, I also computed the the amplitude of the 1st harmonic, which is interestingly 0.0000001605530 dBFS (or 1.605530E-7 dBFS in scientific notation) compared to the peak value $2^{15}-1$ at 0 dBFS. I think the rounding at the top of the sine wave is "pulling it up", giving it a higher amplitude of the fundamental.

Aliasing

As a side note, what we have learned from the results already is that for such a high $N$ the amplitudes of the harmonics decay very slowly as function of harmonic number, so with rational $f/f_s,$ higher harmonics aliased to a lower harmonic are likely to be a majority contributor to the amplitude of the lower harmonic. As an example of such aliasing, with $f_s =$ 48 kHz and $f =$ 1000 Hz, the 5th harmonic gets added to by the 43th, 53th, 91th, 101th, 139th, 149th, etc. harmonic:

$$\begin{eqnarray} && \ldots\\ &=& 48000\text{ Hz}\times -3 + 1000\text{ Hz}\times149\\ &=& 48000\text{ Hz}\times -2 + 1000\text{ Hz}\times101\\ &=& 48000\text{ Hz}\times -1 + 1000\text{ Hz}\times53\\ &=& 1000\text{ Hz}\\ &=& 48000\text{ Hz}\times 1 - 1000\text{ Hz}\times43\\ &=& 48000\text{ Hz}\times 2 - 1000\text{ Hz}\times91\\ &=& 48000\text{ Hz}\times 3 - 1000\text{ Hz}\times139\\ &=& \ldots\end{eqnarray}\tag{4}$$

The effect of the addition of the aliases depends also on the phase of the sinusoid.

24-bit quantization

Then back to limit $f/f_s \to 0$ free of aliasing. The following continuation of the earlier Python script computes in high precision your desired numbers for $N=2^{23} - 1$, in the limiting case $f/f_s \to 0$:

N = 2**23-1 # amplitude (integer)

n = 5  # harmonic number (odd integer)
mp.mp.prec = 70  # Precision 1
20*mp.log10(mp.fabs(a_n(N, n))/N) # in dBFS
mp.mp.prec = 140  # Precision 2
20*mp.log10(mp.fabs(a_n(N, n))/N) # in dBFS

n = 7  # harmonic number (odd integer)
mp.mp.prec = 70  # Precision 1
20*mp.log10(mp.fabs(a_n(N, n))/N) # in dB
mp.mp.prec = 140  # Precision 2
20*mp.log10(mp.fabs(a_n(N, n))/N) # in dB

mp.mp.prec = 53  # Restore default precision

After a total of a few hours of computation on my PC, the result is that with $N=2^{23} - 1$ and $f/f_s \to 0,$ the amplitude of the 5th harmonic is -226.91150085 dBFS and the amplitude of the 7th harmonic is -226.9115030 dBFS. Apparently, the amplitudes decay extremely slowly as function of harmonic number. I also computed the amplitude of the 1st harmonic, 3.9195785E-11 dBFS.

Olli Niemitalo
  • 13.7k
  • 1
  • 35
  • 63