where $f$ is the frequency of the sinusoid, $f_s$ is the sampling frequency, $L$ is the fundamental period of the quantized sequence, and $A$ is the amplitude of the sinusoid. Decibels full-scale (dBFS) shall be referenced to the sinusoid amplitude $A$ at 0 dBFS.
Let the zero-to-peak amplitude of the continuous-time prototypical sinusoid, which is being sampled and quantized, be an integera positive real number $A$. In the limit $f/f_s \to 0,$ where $f$ is the sinusoid frequency and $f_s$ is the sampling frequency, the amplitudes of the harmonics are found in the Fourier series representation of the continuous-time piece-wise constant quantized waveform. The terms of the Fourier series are directly the harmonics of the fundamental frequency. We can construct the waveform as a sum of components like the red curve here, illustrated for sinusoid amplitude $A = 7$:
Figure 1. prototypical sinusoid, of amplitude $A=7$, and its quantization in blue, with a component highlighted in red.
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 of a component of integer amplitude $k \in 1\ldots \operatorname{round}(A)$ are given by:
with integerthe component having the said amplitude $k \in 1\ldots A$ in rangeover the interval $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 non-zero pieces of the component in the other quarters, because they will contribute identically to the odd harmonic cosine terms in the Fourier series. Instead I simply included the implicit factor 4 in the equation.
$$\begin{gather}0 < x < \operatorname{acos}\left(\frac{A - 0.5}{A}\right),\quad \text{if }k = A,\\ \operatorname{acos}\left(\frac{k + 0.5}{A}\right) < x < \operatorname{acos}\left(\frac{k - 0.5}{A}\right),\quad \text{if }k \in 1\ldots A-1,\end{gather}\tag{2}$$$$\begin{gather}0 < x < \operatorname{acos}\left(\frac{\operatorname{round}(A) - 0.5}{A}\right),\quad \text{if }k = \operatorname{round}(A),\\ \operatorname{acos}\left(\frac{k + 0.5}{A}\right) < x < \operatorname{acos}\left(\frac{k - 0.5}{A}\right),\quad \text{if }k \in 1\ldots \operatorname{round}(A)-1,\end{gather}\tag{2}$$
where $\operatorname{acos}\left(\frac{s}{A}\right)$ comes from solving $A\cos(x) = s$ in the first quadrantquarter-period. The non-zero coefficients of the Fourier series of the full piece-wise constant waveform are then given by:
$$\begin{gather}a_n = \frac{4A}{\pi n}\sin\bigg(\operatorname{acos}\left(\tfrac{A - 0.5}{A}\right)n\bigg) +\\\frac{4}{\pi n}\sum_{k=1}^{A-1} \Bigg(k\sin\bigg(\operatorname{acos}\left(\tfrac{k - 0.5}{A}\right)n\bigg) - k\sin\bigg(\operatorname{acos}\left(\tfrac{k + 0.5}{A}\right)n\bigg)\Bigg),\quad n\text{ odd}.\end{gather}\tag{3}$$$$\begin{gather}a_n = \frac{4\operatorname{round}(A)}{\pi n}\,\sin\bigg(\operatorname{acos}\left(\tfrac{\operatorname{round}(A) - 0.5}{A}\right)n\bigg) +\\\frac{4}{\pi n}\sum_{k=1}^{\operatorname{round}(A)-1} \Bigg(k\sin\bigg(\operatorname{acos}\left(\tfrac{k - 0.5}{A}\right)n\bigg) - k\sin\bigg(\operatorname{acos}\left(\tfrac{k + 0.5}{A}\right)n\bigg)\Bigg),\quad n\text{ odd}.\end{gather}\tag{3}$$
The rounding convention does not matter as long as it is to a nearest integer.
It is not out of reach to compute Eq. 3 for $A=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 8-bit sinusoid first, because we can then evaluate $a_n$ at multiple $n$ in a reasonable time and perhaps learn something from the results. In Python:
import numpy as np
import mpmath as mp
def a_n(A, n):
if int(n) & 1:
return 4*A4*mp.floor(A + 0.5)/(mp.pi*n)*mp.sin(mp.acos((mp.floor(A + 0.5)-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, mp.floor(A + 0.5)-1])
else:
return 0
A = 2**7-1 # amplitude (integerpositive real)
for n in np.arange(1, 10**4, 2):
print(str(n)+","+str(20*mp.log10(mp.fabs(a_n(A, n))/A))) # n, a_n (dBFS)
Figure 2. Amplitudes of harmonics of the quantized prototypical sinusoid, with of amplitude $A=2^{7}-1$. Only odd harmonics are present.
It took some minutes on my PC to get the results (Fig. 2). The graph looks jagged. I think this is because the waveform is similar to a sawtooth riding on a sinusoid, and near the zero crossing the sawtooth frequency stays approximately the same over a large portion of the waveform, being $2\pi A \approx$ 797.964534 times the fundamental frequency at the zero crossing as an asymptotic approximation in the limit $A\to\infty$ of the frequency $\pi/\operatorname{asin}\big((k + 0.5)/A\big)$ calculated based on the distance between the discontinuities of the saw wave just around the zero crossing. Indeed, the largest amplitude is found nearby at $|a_{787}| =$ -68.30 dBFS. This is much higher than the amplitude of a "random" early harmonic. Such peaks are repeated later at approximate harmonics of the sawtooth frequency but less clearly and with a roll off arising from the factor $1/n$ in Eq. 3, making the asymptotic decay -20 dB/decade (~-6 dB/octave) as $n\to\infty$. The amplitude of the fundamental frequency is 0.00066 dBFS, I think because the rounding at the top of the sinusoid is "pulling it up" from the 0 dBFS amplitude of the prototypical cosine.
A = 2**23-1 # amplitude (integerpositive real)
n = 5 # number of the harmonic (odd positive integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
n = 7 # number of the harmonic (odd positive integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 53 # Restore default precision
A = 2**23-1 # amplitude (integerpositive real)
for m in range(1, 20):
n = mp.mpf(10)**m+1 # number of the harmonic
mp.mp.prec = 70 # precision
a_n_low_prec = a_n(A, n)
mp.mp.prec = 90 # precision
a_n_high_prec = a_n(A, n)
print(str(n)+","+str(20*mp.log10(mp.fabs(a_n_low_prec)/A))+","+str(20*mp.log10(mp.fabs(a_n_high_prec)/A))) # n, a_n
What we have learned from the results so far is that the amplitudes of the harmonics of the quantized prototypical sinusoid do not really start to decayroll off until after the approximately $2\pi A$th harmonic. This means that with rational $f/f_s$, the amplitude of an in-band harmonic frequency of the quantized sequence is typically almost fully determined by aliased frequencies. As an example of such aliasing, with $f_s =$ 48 kHz and $f =$ 1000 Hz, the 5th harmonic gets summed with the aliases of the 43th, 53th, 91th, 101th, 139th, 149th, etc. harmonic:
where $a_n$ is given by Eq. 3 for odd $n$ and $a_n$ is zero for even $n$. With $f = $ 1000 Hz, $f_s = $ 48 kHz, on my PC it would take perhaps thousands of years to compute Eq. 6 to reasonable accuracy for $A=2^{23}-1$. I also have difficulties evaluating the series even with $A=2^3-1$, with the sequence of partial sums oscillating. I think this difficulty arises from that with $A = 2^3 - 1$ and $f/f_s = 1/48$, the sequence before quantization contains values $-3.5$ and $3.5$ exactly at a discontinuity in the quantized prototypical sinusoid. In any case, with rational $f/f_s$, the real-DFT-based approach to determine the amplitudes of the harmonics is much more practical. It would also work with non-integer amplitude $A$, which Eq. 2 cannot handle in its current form.
import numpy as np
from sympy import Rational
def bins(A, c, d):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (integerpositive real)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)] # real DFT magnitudes (dBFS)
We have encountered here a troublesome combination of $L = 48$ and $A = 2^{23}-1$. We have samples $A\cos(8\times2\pi/L) =$ $A\cos(40\times2\pi/L) =$ $A/2$ and $A\cos(16\times2\pi/L) =$ $A\cos(32\times2\pi/L) =$ $-A/2$, and as $A$ ishappens to be an odd integer, those samples will always have a fractional part of $1/2$. The above reported amplitudes of the harmonics may thus be dependent on the rounding mode, or, the computation of the phase and/or evaluation of the cosine may have given a small numerical error which has determined the direction of rounding. This may even break the symmetry properties of the sequence, affecting which amplitudes of the harmonics of the period-$L$ frequency are zero.
import numpy as np
from sympy import Rational
def bins(A, c, d):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
return np.fft.rfft(waveform)/(A*L/2)
A = 8 # amplitude (integerpositive reals)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)] # real DFT magnitudes (dBFS)
import mpmath as mp
def a_n(A, n):
if int(n) & 1:
return 4*A4*mp.floor(A + 0.5)/(mp.pi*n)*mp.sin(mp.acos((mp.floor(A + 0.5)-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, mp.floor(A + 0.5)-1])
else:
return 0
mp.mp.prec = 100
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
if c == np.gcd(c, d):
for n in range(L/2):
bin = a_n(A, n) + mp.nsum(lambda k: a_n(A, k*L - n) + a_n(A, k*L + n), [1, int(10000*mp.pi*A)])
float(20*mp.log10(mp.fabs(bin)/A)) # harmonic amplitude (dBFS)
else:
print("Error: fundamental frequency must be in bin_1")
mp.mp.prec = 53
import numpy as np
from sympy import Rational
def bins(A, c, d, phi):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (integerpositive real)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
phi = 0.123 # phase shift of the cosine (radians)
[20*np.log10(float(abs(x))) for x in bins(A, c, d, phi)] # real DFT magnitudes (dBFS)
import numpy as np
from sympy import Rational
def E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
num_bins = int((L + 2)/2)
sum_sq_bins = np.zeros(num_bins)
max_sq_bins = np.zeros(num_bins)
for n in range(N):
phi = float((0.5 + n/N)*2*np.pi/L)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
sq_bins = np.square(np.absolute(np.fft.rfft(waveform)/(A*L/2)))
sum_sq_bins = sum_sq_bins + sq_bins
max_sq_bins = np.maximum(max_sq_bins, sq_bins)
return [sum_sq_bins/N, max_sq_bins]
A = 2**23-1 # amplitude (integerpositive real)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
N = 2**18 # size of statistical sample (positive integer)
results = E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N)
[10*np.log10(float(x)) for x in results[0]] # est. expected value of real DFT squared magnitudes (dBFS)
[10*np.log10(float(x)) for x in results[1]] # est. max real DFT squared magnitudes (dBFS)
import numpy as np
from sympy import Rational
def bins_max_err_seq(A, c, d, n):
L = int(Rational(d, np.gcd(c, d))) # Fundamental period (samples)
waveform = 0.5*np.sign(np.sin(np.dot(range(L), 2*np.pi*n/L)))
for k in range(L):
if 2*k*n % L == 0:
waveform[k] = 0.5*np.cos(k*2*np.pi*n/L)
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (integerpositive real)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
if c == np.gcd(c, d) and int(Rational(d, np.gcd(c, d))) % 1 == 0:
for n in np.arange(1, int(d/c/2 + 1), 2):
print(str(n) + "," + str(20*np.log10(float(abs(bins_max_err_seq(A, c, d, n)[n]))))) # real DFT magnitude of bin n for the maximum quantization error sequence (dBFS)
else:
print("Error: Fundamental frequency must be in bin 1 and L must be even\n")
This answer has reached the maximum allowed length of 30000 characters.
Let the zero-to-peak amplitude of the continuous-time prototypical sinusoid, which is being sampled and quantized, be an integer $A$. In the limit $f/f_s \to 0,$ where $f$ is the sinusoid frequency and $f_s$ is the sampling frequency, the amplitudes of the harmonics are found in the Fourier series representation of the continuous-time piece-wise constant quantized waveform. The terms of the Fourier series are directly the harmonics of the fundamental frequency. We can construct the waveform as a sum of components like the red curve here, illustrated for sinusoid amplitude $A = 7$:
Figure 1. prototypical sinusoid, of amplitude $A=7$, and its quantization in blue, with a component highlighted in red.
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:
with integer amplitude $k \in 1\ldots A$ 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 harmonic cosine terms in the Fourier series. Instead I simply included the implicit factor 4 in the equation.
$$\begin{gather}0 < x < \operatorname{acos}\left(\frac{A - 0.5}{A}\right),\quad \text{if }k = A,\\ \operatorname{acos}\left(\frac{k + 0.5}{A}\right) < x < \operatorname{acos}\left(\frac{k - 0.5}{A}\right),\quad \text{if }k \in 1\ldots A-1,\end{gather}\tag{2}$$
where $\operatorname{acos}\left(\frac{s}{A}\right)$ comes from solving $A\cos(x) = s$ in the first quadrant. The non-zero coefficients of the Fourier series of the full piece-wise constant waveform are then given by:
$$\begin{gather}a_n = \frac{4A}{\pi n}\sin\bigg(\operatorname{acos}\left(\tfrac{A - 0.5}{A}\right)n\bigg) +\\\frac{4}{\pi n}\sum_{k=1}^{A-1} \Bigg(k\sin\bigg(\operatorname{acos}\left(\tfrac{k - 0.5}{A}\right)n\bigg) - k\sin\bigg(\operatorname{acos}\left(\tfrac{k + 0.5}{A}\right)n\bigg)\Bigg),\quad n\text{ odd}.\end{gather}\tag{3}$$
It is not out of reach to compute Eq. 3 for $A=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 8-bit sinusoid first, because we can then evaluate $a_n$ at multiple $n$ in a reasonable time:
import numpy as np
import mpmath as mp
def a_n(A, n):
if int(n) & 1:
return 4*A/(mp.pi*n)*mp.sin(mp.acos((A-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, A-1])
else:
return 0
A = 2**7-1 # amplitude (integer)
for n in np.arange(1, 10**4, 2):
print(str(n)+","+str(20*mp.log10(mp.fabs(a_n(A, n))/A))) # n, a_n (dBFS)
Figure 2. Amplitudes of harmonics of the quantized prototypical sinusoid, with $A=2^{7}-1$. Only odd harmonics are present.
It took some minutes on my PC to get the results (Fig. 2). The graph looks jagged. I think this is because the waveform is similar to a sawtooth riding on a sinusoid, and near the zero crossing the sawtooth frequency stays approximately the same over a large portion of the waveform, being $2\pi A \approx$ 797.964534 times the fundamental frequency at the zero crossing. Indeed, the largest amplitude is found nearby at $|a_{787}| =$ -68.30 dBFS. This is much higher than the amplitude of a "random" early harmonic. Such peaks are repeated later at approximate harmonics of the sawtooth frequency but less clearly and with a roll off arising from the factor $1/n$ in Eq. 3, making the asymptotic decay -20 dB/decade (~-6 dB/octave) as $n\to\infty$. The amplitude of the fundamental frequency is 0.00066 dBFS, I think because the rounding at the top of the sinusoid is "pulling it up" from the 0 dBFS amplitude of the prototypical cosine.
A = 2**23-1 # amplitude (integer)
n = 5 # number of the harmonic (odd integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
n = 7 # number of the harmonic (odd integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 53 # Restore default precision
A = 2**23-1 # amplitude (integer)
for m in range(1, 20):
n = mp.mpf(10)**m+1 # number of the harmonic
mp.mp.prec = 70 # precision
a_n_low_prec = a_n(A, n)
mp.mp.prec = 90 # precision
a_n_high_prec = a_n(A, n)
print(str(n)+","+str(20*mp.log10(mp.fabs(a_n_low_prec)/A))+","+str(20*mp.log10(mp.fabs(a_n_high_prec)/A))) # n, a_n
What we have learned from the results so far is that the amplitudes of the harmonics of the quantized prototypical sinusoid do not really start to decay until after the approximately $2\pi A$th harmonic. This means that with rational $f/f_s$, the amplitude of an in-band harmonic frequency of the quantized sequence is typically almost fully determined by aliased frequencies. As an example of such aliasing, with $f_s =$ 48 kHz and $f =$ 1000 Hz, the 5th harmonic gets summed with the aliases of the 43th, 53th, 91th, 101th, 139th, 149th, etc. harmonic:
where $a_n$ is given by Eq. 3 for odd $n$ and $a_n$ is zero for even $n$. With $f = $ 1000 Hz, $f_s = $ 48 kHz, on my PC it would take perhaps thousands of years to compute Eq. 6 to reasonable accuracy for $A=2^{23}-1$. I also have difficulties evaluating the series even with $A=2^3-1$, with the sequence of partial sums oscillating. I think this difficulty arises from that with $A = 2^3 - 1$ and $f/f_s = 1/48$, the sequence before quantization contains values $-3.5$ and $3.5$ exactly at a discontinuity in the quantized prototypical sinusoid. In any case, with rational $f/f_s$, the real-DFT-based approach to determine the amplitudes of the harmonics is much more practical. It would also work with non-integer amplitude $A$, which Eq. 2 cannot handle in its current form.
import numpy as np
from sympy import Rational
def bins(A, c, d):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (integer)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)] # real DFT magnitudes (dBFS)
We have encountered here a troublesome combination of $L = 48$ and $A = 2^{23}-1$. We have samples $A\cos(8\times2\pi/L) =$ $A\cos(40\times2\pi/L) =$ $A/2$ and $A\cos(16\times2\pi/L) =$ $A\cos(32\times2\pi/L) =$ $-A/2$, and as $A$ is odd, those samples will always have a fractional part of $1/2$. The above reported amplitudes of the harmonics may thus be dependent on the rounding mode, or, the computation of the phase and/or evaluation of the cosine may have given a small numerical error which has determined the direction of rounding. This may even break the symmetry properties of the sequence, affecting which amplitudes of the harmonics of the period-$L$ frequency are zero.
import numpy as np
from sympy import Rational
def bins(A, c, d):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
return np.fft.rfft(waveform)/(A*L/2)
A = 8 # amplitude (integer)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)] # real DFT magnitudes (dBFS)
import mpmath as mp
def a_n(A, n):
if int(n) & 1:
return 4*A/(mp.pi*n)*mp.sin(mp.acos((A-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, A-1])
else:
return 0
mp.mp.prec = 100
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
if c == np.gcd(c, d):
for n in range(L/2):
bin = a_n(A, n) + mp.nsum(lambda k: a_n(A, k*L - n) + a_n(A, k*L + n), [1, int(10000*mp.pi*A)])
float(20*mp.log10(mp.fabs(bin)/A)) # harmonic amplitude (dBFS)
else:
print("Error: fundamental frequency must be in bin_1")
mp.mp.prec = 53
import numpy as np
from sympy import Rational
def bins(A, c, d, phi):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (integer)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
phi = 0.123 # phase shift of the cosine (radians)
[20*np.log10(float(abs(x))) for x in bins(A, c, d, phi)] # real DFT magnitudes (dBFS)
import numpy as np
from sympy import Rational
def E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
num_bins = int((L + 2)/2)
sum_sq_bins = np.zeros(num_bins)
max_sq_bins = np.zeros(num_bins)
for n in range(N):
phi = float((0.5 + n/N)*2*np.pi/L)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
sq_bins = np.square(np.absolute(np.fft.rfft(waveform)/(A*L/2)))
sum_sq_bins = sum_sq_bins + sq_bins
max_sq_bins = np.maximum(max_sq_bins, sq_bins)
return [sum_sq_bins/N, max_sq_bins]
A = 2**23-1 # amplitude (integer)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
N = 2**18 # size of statistical sample (positive integer)
results = E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N)
[10*np.log10(float(x)) for x in results[0]] # est. expected value of real DFT squared magnitudes (dBFS)
[10*np.log10(float(x)) for x in results[1]] # est. max real DFT squared magnitudes (dBFS)
import numpy as np
from sympy import Rational
def bins_max_err_seq(A, c, d, n):
L = int(Rational(d, np.gcd(c, d))) # Fundamental period (samples)
waveform = 0.5*np.sign(np.sin(np.dot(range(L), 2*np.pi*n/L)))
for k in range(L):
if 2*k*n % L == 0:
waveform[k] = 0.5*np.cos(k*2*np.pi*n/L)
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (integer)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
if c == np.gcd(c, d) and int(Rational(d, np.gcd(c, d))) % 1 == 0:
for n in np.arange(1, int(d/c/2 + 1), 2):
print(str(n) + "," + str(20*np.log10(float(abs(bins_max_err_seq(A, c, d, n)[n]))))) # real DFT magnitude of bin n for the maximum quantization error sequence (dBFS)
else:
print("Error: Fundamental frequency must be in bin 1 and L must be even\n")
Let the zero-to-peak amplitude of the continuous-time prototypical sinusoid, which is being sampled and quantized, be a positive real number $A$. In the limit $f/f_s \to 0,$ where $f$ is the sinusoid frequency and $f_s$ is the sampling frequency, the amplitudes of the harmonics are found in the Fourier series representation of the continuous-time piece-wise constant quantized waveform. The terms of the Fourier series are directly the harmonics of the fundamental frequency. We can construct the waveform as a sum of components like the red curve here, illustrated for sinusoid amplitude $A = 7$:
Figure 1. prototypical sinusoid, of amplitude $A=7$, and its quantization in blue, with a component highlighted in red.
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 of a component of integer amplitude $k \in 1\ldots \operatorname{round}(A)$ are given by:
with the component having the said amplitude over the interval $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 non-zero pieces of the component in the other quarters, because they will contribute identically to the odd harmonic cosine terms in the Fourier series. Instead I simply included the implicit factor 4 in the equation.
$$\begin{gather}0 < x < \operatorname{acos}\left(\frac{\operatorname{round}(A) - 0.5}{A}\right),\quad \text{if }k = \operatorname{round}(A),\\ \operatorname{acos}\left(\frac{k + 0.5}{A}\right) < x < \operatorname{acos}\left(\frac{k - 0.5}{A}\right),\quad \text{if }k \in 1\ldots \operatorname{round}(A)-1,\end{gather}\tag{2}$$
where $\operatorname{acos}\left(\frac{s}{A}\right)$ comes from solving $A\cos(x) = s$ in the first quarter-period. The non-zero coefficients of the Fourier series of the full piece-wise constant waveform are then given by:
$$\begin{gather}a_n = \frac{4\operatorname{round}(A)}{\pi n}\,\sin\bigg(\operatorname{acos}\left(\tfrac{\operatorname{round}(A) - 0.5}{A}\right)n\bigg) +\\\frac{4}{\pi n}\sum_{k=1}^{\operatorname{round}(A)-1} \Bigg(k\sin\bigg(\operatorname{acos}\left(\tfrac{k - 0.5}{A}\right)n\bigg) - k\sin\bigg(\operatorname{acos}\left(\tfrac{k + 0.5}{A}\right)n\bigg)\Bigg),\quad n\text{ odd}.\end{gather}\tag{3}$$
The rounding convention does not matter as long as it is to a nearest integer.
It is not out of reach to compute Eq. 3 for $A=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 8-bit sinusoid first, because we can then evaluate $a_n$ at multiple $n$ in a reasonable time and perhaps learn something from the results. In Python:
import numpy as np
import mpmath as mp
def a_n(A, n):
if int(n) & 1:
return 4*mp.floor(A + 0.5)/(mp.pi*n)*mp.sin(mp.acos((mp.floor(A + 0.5)-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, mp.floor(A + 0.5)-1])
else:
return 0
A = 2**7-1 # amplitude (positive real)
for n in np.arange(1, 10**4, 2):
print(str(n)+","+str(20*mp.log10(mp.fabs(a_n(A, n))/A))) # n, a_n (dBFS)
Figure 2. Amplitudes of harmonics of the quantized prototypical sinusoid of amplitude $A=2^{7}-1$. Only odd harmonics are present.
It took some minutes on my PC to get the results (Fig. 2). The graph looks jagged. I think this is because the waveform is similar to a sawtooth riding on a sinusoid, and near the zero crossing the sawtooth frequency stays approximately the same over a large portion of the waveform, being $2\pi A \approx$ 797.964534 times the fundamental frequency at the zero crossing as an asymptotic approximation in the limit $A\to\infty$ of the frequency $\pi/\operatorname{asin}\big((k + 0.5)/A\big)$ calculated based on the distance between the discontinuities of the saw wave just around the zero crossing. Indeed, the largest amplitude is found nearby at $|a_{787}| =$ -68.30 dBFS. This is much higher than the amplitude of a "random" early harmonic. Such peaks are repeated later at approximate harmonics of the sawtooth frequency but less clearly and with a roll off arising from the factor $1/n$ in Eq. 3, making the asymptotic decay -20 dB/decade (~-6 dB/octave) as $n\to\infty$. The amplitude of the fundamental frequency is 0.00066 dBFS, I think because the rounding at the top of the sinusoid is "pulling it up" from the 0 dBFS amplitude of the prototypical cosine.
A = 2**23-1 # amplitude (positive real)
n = 5 # number of the harmonic (odd positive integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
n = 7 # number of the harmonic (odd positive integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 53 # Restore default precision
A = 2**23-1 # amplitude (positive real)
for m in range(1, 20):
n = mp.mpf(10)**m+1 # number of the harmonic
mp.mp.prec = 70 # precision
a_n_low_prec = a_n(A, n)
mp.mp.prec = 90 # precision
a_n_high_prec = a_n(A, n)
print(str(n)+","+str(20*mp.log10(mp.fabs(a_n_low_prec)/A))+","+str(20*mp.log10(mp.fabs(a_n_high_prec)/A))) # n, a_n
What we have learned from the results so far is that the amplitudes of the harmonics of the quantized prototypical sinusoid do not start to roll off until after the approximately $2\pi A$th harmonic. This means that with rational $f/f_s$, the amplitude of an in-band harmonic frequency of the quantized sequence is typically almost fully determined by aliased frequencies. As an example of such aliasing, with $f_s =$ 48 kHz and $f =$ 1000 Hz, the 5th harmonic gets summed with the aliases of the 43th, 53th, 91th, 101th, 139th, 149th, etc. harmonic:
where $a_n$ is given by Eq. 3 for odd $n$ and $a_n$ is zero for even $n$. With $f = $ 1000 Hz, $f_s = $ 48 kHz, on my PC it would take perhaps thousands of years to compute Eq. 6 to reasonable accuracy for $A=2^{23}-1$. I also have difficulties evaluating the series even with $A=2^3-1$, with the sequence of partial sums oscillating. I think this difficulty arises from that with $A = 2^3 - 1$ and $f/f_s = 1/48$, the sequence before quantization contains values $-3.5$ and $3.5$ exactly at a discontinuity in the quantized prototypical sinusoid. In any case, with rational $f/f_s$, the real-DFT-based approach to determine the amplitudes of the harmonics is much more practical.
import numpy as np
from sympy import Rational
def bins(A, c, d):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (positive real)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)] # real DFT magnitudes (dBFS)
We have encountered here a troublesome combination of $L = 48$ and $A = 2^{23}-1$. We have samples $A\cos(8\times2\pi/L) =$ $A\cos(40\times2\pi/L) =$ $A/2$ and $A\cos(16\times2\pi/L) =$ $A\cos(32\times2\pi/L) =$ $-A/2$, and as $A$ happens to be an odd integer, those samples will always have a fractional part of $1/2$. The above reported amplitudes of the harmonics may thus be dependent on the rounding mode, or, the computation of the phase and/or evaluation of the cosine may have given a small numerical error which has determined the direction of rounding. This may even break the symmetry properties of the sequence, affecting which amplitudes of the harmonics of the period-$L$ frequency are zero.
import numpy as np
from sympy import Rational
def bins(A, c, d):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
return np.fft.rfft(waveform)/(A*L/2)
A = 8 # amplitude (positive reals)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)] # real DFT magnitudes (dBFS)
import mpmath as mp
def a_n(A, n):
if int(n) & 1:
return 4*mp.floor(A + 0.5)/(mp.pi*n)*mp.sin(mp.acos((mp.floor(A + 0.5)-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, mp.floor(A + 0.5)-1])
else:
return 0
mp.mp.prec = 100
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
if c == np.gcd(c, d):
for n in range(L/2):
bin = a_n(A, n) + mp.nsum(lambda k: a_n(A, k*L - n) + a_n(A, k*L + n), [1, int(10000*mp.pi*A)])
float(20*mp.log10(mp.fabs(bin)/A)) # harmonic amplitude (dBFS)
else:
print("Error: fundamental frequency must be in bin_1")
mp.mp.prec = 53
import numpy as np
from sympy import Rational
def bins(A, c, d, phi):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (positive real)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
phi = 0.123 # phase shift of the cosine (radians)
[20*np.log10(float(abs(x))) for x in bins(A, c, d, phi)] # real DFT magnitudes (dBFS)
import numpy as np
from sympy import Rational
def E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N):
L = Rational(d, np.gcd(c, d)) # Fundamental period (samples)
num_bins = int((L + 2)/2)
sum_sq_bins = np.zeros(num_bins)
max_sq_bins = np.zeros(num_bins)
for n in range(N):
phi = float((0.5 + n/N)*2*np.pi/L)
waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
sq_bins = np.square(np.absolute(np.fft.rfft(waveform)/(A*L/2)))
sum_sq_bins = sum_sq_bins + sq_bins
max_sq_bins = np.maximum(max_sq_bins, sq_bins)
return [sum_sq_bins/N, max_sq_bins]
A = 2**23-1 # amplitude (positive real)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
N = 2**18 # size of statistical sample (positive integer)
results = E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N)
[10*np.log10(float(x)) for x in results[0]] # est. expected value of real DFT squared magnitudes (dBFS)
[10*np.log10(float(x)) for x in results[1]] # est. max real DFT squared magnitudes (dBFS)
import numpy as np
from sympy import Rational
def bins_max_err_seq(A, c, d, n):
L = int(Rational(d, np.gcd(c, d))) # Fundamental period (samples)
waveform = 0.5*np.sign(np.sin(np.dot(range(L), 2*np.pi*n/L)))
for k in range(L):
if 2*k*n % L == 0:
waveform[k] = 0.5*np.cos(k*2*np.pi*n/L)
return np.fft.rfft(waveform)/(A*L/2)
A = 2**23-1 # amplitude (positive real)
c = 1000 # numerator of f/f_s (positive integer)
d = 48000 # denominator of f/f_s (positive integer)
if c == np.gcd(c, d) and int(Rational(d, np.gcd(c, d))) % 1 == 0:
for n in np.arange(1, int(d/c/2 + 1), 2):
print(str(n) + "," + str(20*np.log10(float(abs(bins_max_err_seq(A, c, d, n)[n]))))) # real DFT magnitude of bin n for the maximum quantization error sequence (dBFS)
else:
print("Error: Fundamental frequency must be in bin 1 and L must be even\n")
This answer has reached the maximum allowed length of 30000 characters.
Figure 2. Amplitudes of a selection of odd harmonics of the quantized prototypical sinusoid, with $A=2^{7}-1$. Only odd harmonics are present.
It took some minutes on my PC to get the results (Fig. 2). The graph looks jagged. I think this is because the waveform is similar to a sawtooth riding on a sinusoid, and near the zero crossing the sawtooth frequency stays approximately the same over a large portion of the waveform, being $2\pi A \approx$ 797.964534 times the fundamental frequency at the zero crossing. Indeed, the largest amplitude is found nearby at $|a_{787}| =$ -68.30 dBFS. This is much higher than the amplitude of a "random" early harmonic. Such peaks are repeated later at approximate harmonics of the sawtooth frequency but less clearly and with a roll off arising from the factor $1/n$ in Eq. 3, making the asymptotic decay -20 dB/decade (~-6 dB/octave) as $n\to\infty$. The amplitude of the fundamental frequency is 0.00066 dBFS, I think because the rounding at the top of the sinusoid is "pulling it up" from the 0 dBFS amplitude of the prototypical cosine.
Figure 2. Amplitudes of a selection of odd harmonics of the quantized prototypical sinusoid, with $A=2^{7}-1$.
It took some minutes on my PC to get the results (Fig. 2). The graph looks jagged. I think this is because the waveform is similar to a sawtooth riding on a sinusoid, and near the zero crossing the sawtooth frequency stays approximately the same over a large portion of the waveform being $2\pi A \approx$ 797.964534 times the fundamental frequency at the zero crossing. Indeed, the largest amplitude is found nearby at $|a_{787}| =$ -68.30 dBFS. This is much higher than the amplitude of a "random" early harmonic. Such peaks are repeated later at approximate harmonics of the sawtooth frequency but less clearly and with a roll off arising from the factor $1/n$ in Eq. 3, making the asymptotic decay -20 dB/decade (~-6 dB/octave) as $n\to\infty$. The amplitude of the fundamental frequency is 0.00066 dBFS, I think because the rounding at the top of the sinusoid is "pulling it up" from the 0 dBFS amplitude of the prototypical cosine.
Figure 2. Amplitudes of harmonics of the quantized prototypical sinusoid, with $A=2^{7}-1$. Only odd harmonics are present.
It took some minutes on my PC to get the results (Fig. 2). The graph looks jagged. I think this is because the waveform is similar to a sawtooth riding on a sinusoid, and near the zero crossing the sawtooth frequency stays approximately the same over a large portion of the waveform, being $2\pi A \approx$ 797.964534 times the fundamental frequency at the zero crossing. Indeed, the largest amplitude is found nearby at $|a_{787}| =$ -68.30 dBFS. This is much higher than the amplitude of a "random" early harmonic. Such peaks are repeated later at approximate harmonics of the sawtooth frequency but less clearly and with a roll off arising from the factor $1/n$ in Eq. 3, making the asymptotic decay -20 dB/decade (~-6 dB/octave) as $n\to\infty$. The amplitude of the fundamental frequency is 0.00066 dBFS, I think because the rounding at the top of the sinusoid is "pulling it up" from the 0 dBFS amplitude of the prototypical cosine.