# Precise 5th and 7th harmonics of a sampled sine wave

Does anyone know in decibels (to 1/100th of a dB) what the theoretical 3rd, 5th and 7th harmonic of a 0dB fs 24-bit (i.e. full-level; 0dB = -8,388,607 to 8,388,607) sampled sine wave without dither will be? Where the harmonics are a product of the sampling process (caused by the effect of the discrete sampled steps). Better still, how I calculate the theoretical values for the harmonic components.

The fundamental is around 1,000Hz (happy to use exactly 1kHz). The sampling frequency is 48kHz. Asynchronous sampling.

• This is kind of a "do my homework for me" question – but that's not too bad, but I might warn you that we'll with a high chance answer based on an assumption that you need to understand what to do more than a need to get a number down to 0.01 dB (a scale which almost certainly makes no sense at all). Mar 20 '19 at 6:32
• Problem: these harmonics won't be the same for different sine frequencies – their power isn't just a linear (or any) function of the original sine's amplitude, but also its frequency. You're not giving the frequency, but demanding a high-accuracy number – makes no sense, if you ask me. I think your question is really only worth answering if you explain what you'd need the harmonic-power-over-frequency curve for, based on which we could deduct, together, how to approach this. Mar 20 '19 at 6:34
• Hi. I am running a test and can use any sampling frequency, fundamental resolution etc. I have written an algorithm Mar 20 '19 at 6:41
• OK, then we can probably just pick an arbitrary set of powers for the harmonics you ask for and just let you look for exactly the frequency that generates these. So, please explain your test case in your question (edit it, please), and also, the algorithm! Mar 20 '19 at 6:44
• Edit your question to include crucial information. Don't just add it in comments. Mar 20 '19 at 6:58

This answer discusses the harmonic spectra of the quantized sequence in five cases:

• limit $$f/f_s \to 0$$,
• synchronous sampling of a cosine with rational $$f/f_s$$,
• synchronous sampling of a sinusoid of arbitrary phase with rational $$f/f_s$$,
• asynchronous sampling of a sinusoid with rational $$f/f_s$$, and
• asynchronous sampling of a sinusoid with rational $$f/f_s$$ in the limit $$L/A \to 0$$,

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.

# Limit $$f/f_s \to 0$$

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:

$$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 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.

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{\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.

## 8-bit quantization

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.

## 24-bit quantization

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

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


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. After a total of a few hours of computation on my PC, the result is that with $$A=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. I also computed the amplitude of the fundamental frequency, 3.9195785E-11 dBFS. I also computed in a full day the amplitudes of a number of higher harmonics, in Python:

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 Figure 3. Amplitudes of a selection of odd harmonics of the quantized prototypical sinusoid, with $$A=2^{23}-1$$.

The corner frequency before the roll-off is again approximately the sawtooth frequency at the zero crossing of the quantized prototypical sinusoid, $$2\pi A = 2\pi(2^{23}-1) \approx$$ 5.270717225E7 times the fundamental frequency $$f$$.

# Synchronous sampling of a cosine with rational $$f/f_s$$

Let the prototypical sinusoid, which is being sampled and quantized, be a cosine without a phase shift. That is to say, one of the samples is at the peak of the sinusoid. Let $$f/f_s$$ be equal to a rational number $$c/d$$, with $$c$$ and $$d$$ integer. The fundamental period $$L$$ of the periodic sequence of quantized samples is:

$$L = \frac{d}{\operatorname{gcd}(c, d)},\tag{4}$$

where gcd denotes the greatest common divisor. A real discrete Fourier transform (real DFT) of length $$L$$ has bin frequencies that are the harmonic frequencies of the period-$$L$$ frequency. The period-$$L$$ frequency is equivalent to the frequency of the prototypical sinusoid only if the numerator of the irreducible fraction representation of $$f/f_s$$ equals 1. Otherwise the real DFT bin frequencies contain also other frequencies than those equivalent to in-band harmonics of the sinusoid.

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:

$$\begin{eqnarray} && \ldots\\ &=& 48000\text{ Hz}\times -3 + 1000\text{ Hz}\times 149\\ &=& 48000\text{ Hz}\times -2 + 1000\text{ Hz}\times 101\\ &=& 48000\text{ Hz}\times -1 + 1000\text{ Hz}\times 53\\ &=& 1000\text{ Hz}\times 5\\ &=& 48000\text{ Hz}\times 1 - 1000\text{ Hz}\times 43\\ &=& 48000\text{ Hz}\times 2 - 1000\text{ Hz}\times 91\\ &=& 48000\text{ Hz}\times 3 - 1000\text{ Hz}\times 139\\ &=& \ldots\end{eqnarray}\tag{5}$$

The effect of summation of the harmonics aliased to a frequency depends also on the phase of the sinusoid being quantized, because the phases of the aliased harmonics control whether there is constructive or destructive interference. Aliasing does not change the phase of a harmonic. A phase shift of the quantized prototypical sinusoid is equivalent to shifting it in time, which accounts to a phase shift of each unaliased harmonic proportional to its number. In this answer we will only be sampling and quantizing a prototypical cosine without a phase shift, so to compute a single bin of the real DFT of the discrete-time periodic sequence of samples, it would suffices to sum those $$a$$ coefficients of the Fourier series of the quantized prototypical sinusoid that alias to the equivalent frequency of the bin. We'll use the normalization convention that $$|\operatorname{bin}_1| = 1$$ for the real DFT of a sinusoid of amplitude $$1$$ and of angular frequency $$2\pi/L$$, where $$L$$ is the length of the real DFT in time domain. If the numerator of $$f/f_s$$ expressed as an irreducible fraction equals 1, we have:

$$\operatorname{bin}_n = a_n + \sum_{k=1}^{\infty}\left(a_{kL-n} + a_{kL+n}\right),\tag{6}$$

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.

## 24-bit quantization

The following Python script computes using the real-DFT-based method the amplitudes of the harmonics of a synchronously sampled cosine for $$A = 2^{23}-1$$, $$f =$$ 1000 Hz, and $$f_s =$$ 48 kHz:

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)


While not required by the script, in this case $$c/d = 1/48,$$ so each real DFT bin corresponds to a harmonic of the prototypical sinusoid. The fundamental frequency's amplitude measures as -1.75E-7 dBFS, the 5th harmonic as -160.90 dBFS and the 7th as -160.75 dBFS. Unfortunately I could not find a suitable multiple precision Fast Fourier Transform (FFT) library, but I think those three numbers came out in sufficient precision using NumPy's rfft. The precision is not sufficient to calculate the amplitude of say the 3rd harmonic seemingly at about -333 dBFS (Fig. 4). Figure 4. Magnitudes of the real DFT bins of the quantized sequence of a synchronously sampled cosine with $$A=2^{23}-1$$, $$f =$$ 1000 Hz, and $$f_s =$$ 48 kHz. These results are only applicable when one of the samples is at the peak of the sinusoid, so not applicable to arbitrary phase shifts or asynchronous sampling.

With this $$f/f_s$$ the fundamental period is $$L = 48$$. The computed amplitudes of the harmonics are also correct for phase shifts of the prototypical cosine that are multiples of $$2\pi/L$$, because the cosine is sampled also at those phases, or such a multiple plus $$\pi$$, because of the symmetry properties of cosine. In the presented case the compatible phase shifts of the cosine are multiples of $$7.5°$$.

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.

# Method validation

We can validate the two approaches against each other by choosing a pair of $$A$$ and $$f/f_s$$ that does not result in sampling of the quantized prototypical cosine at any of its discontinuities. Firstly, this removes having to choose arbitrarily which way a number exactly half-way between integers should be rounded, and secondly, should help in convergence of the infinite series in Eq. 6. Let's try $$A = 8$$ and $$f/f_s = 1/48$$. First, the real-DFT-based method, in Python:

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)


This results in the following amplitudes of the harmonics, printed in precision beyond their accuracy, in dBFS:

-inf, 0.005622747208892056, -inf, -331.19944825653926, -inf, -50.41376796795221, -inf, -35.41599672115829, -inf, -327.2878768807577, -inf, -38.14783193548751, -inf, -45.42857299685606, -inf, -324.8188873214554, -inf, -40.41625528549065, -inf, -36.12873038111221, -inf, -337.55682735433885, -inf, -33.04816436790989, -inf


Then, aliasing the harmonics of the quantized prototypical cosine by evaluating Eq. 6, with the series truncated to the first approximately $$10^4\pi A$$ terms, continuing with the imports and constants of the previous Python script:

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


which results in the following amplitudes, printed in precision beyond their accuracy, in dBFS:

-inf, 0.005622452922818016, -inf, -165.35811090601203, -inf, -50.41370472648242, -inf, -35.41599394709451, -inf, -155.34167605811237, -inf, -38.14782298030069, -inf, -45.42859376902997, -inf, -164.6586753088176, -inf, -40.41623886311205, -inf, -36.12872396082294, -inf, -149.78722963456818, -inf, -33.04815743325275


The results seem to agree within their numerical accuracy (Fig 6.). Figure 6. Comparison of amplitudes, in dBFS, of harmonics of a quantized cosine sampled synchronously, with $$A = 8$$ and $$f/f_s = 1/48$$, calculated by the real-DFT-based method (blue squares) and by summing coefficients $$a_n$$ of harmonics that alias to the same frequency (orange crosses).

The 3rd harmonic amplitude calculated by truncating the sum in Eq. 6 converges very slowly as function of the number of terms (Fig. 6). This is probably the case also with the other low-amplitude harmonics. It does not seem feasible to reach the accuracy of the real-DFT-based method using the harmonic aliasing method. Figure 7. $$\operatorname{bin}_3$$ as function of the number of terms in truncation of the sum in Eq. 6.

# Synchronous sampling of a sinusoid of arbitrary phase with rational $$f/f_s$$

## 24-bit quantization

The real-DFT-based method is simple to modify for arbitrary phase shifts of the prototypical cosine. For $$A = 2^{23}-1$$, $$f =$$ 1000 Hz, $$f_s =$$ 48 kHz, and a phase shift $$\phi = 0.123$$ in radians, the amplitudes of the harmonics of the fundamental period can be computed by the following Python script:

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) Figure 8. Magnitudes of the real DFT bins of the quantized sequence of a synchronously sampled cosine with phase shift $$\phi = 0.123$$, $$A=2^{23}−1$$, $$f=$$ 1000 Hz, and $$f_s$$= 48 kHz.

# Asynchronous sampling of a sinusoid with rational $$f/f_s$$

Asynchronous sampling with rational $$f/f_s$$ can be modeled by assuming that $$f/f_s$$ is very close to a rational number. Then, locally the sequence of samples appears periodic with a computable harmonic spectrum. After sufficient time has elapsed, the phase of the prototypical sinusoid will have drifted to a different value and with a different local harmonic spectrum. For a random such sufficient time, the drift is a random variable between $$0$$ and $$2\pi$$. A reasonable probabilistic analog of the amplitude of a real DFT bin of the periodic sequence is the expected value of the square of its absolute value, $$\operatorname{E}[|\operatorname{bin}_n|^2]$$, which can be expressed in the same dBFS scale. We can estimate its value using uniform statistical sampling over phase shifts $$\phi$$ in a non-redundant range $$0 < \phi < 2\pi/L$$, where $$L$$ is the fundamental period of the locally periodic sequence. The end points of the range shall be deliberately left out so that we don't have the previously encountered problem with rounding modes and such.

## 24-bit quantization

$$\operatorname{E}[|\operatorname{bin}_n|^2]$$ is estimated for $$A = 2^{23}-1$$, $$f =$$ 1000 Hz, and $$f_s =$$ 48 kHz, by the following Python script, which also estimates $$\max|\operatorname{bin}_n|$$ over random phase shifts, representing the worst-case non-fundamental harmonic amplitudes in the local spectra:

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]  # est. expected value of real DFT squared magnitudes (dBFS)
[10*np.log10(float(x)) for x in results]  # est. max real DFT squared magnitudes (dBFS)


The results convergence quite slowly as the size of the statistical sample N is increased, and the accuracy requirements stated in the question are almost but not quite fulfilled by a day's calculation on my PC for $$\operatorname{E}[|\operatorname{bin}_n|^2]$$ and not at all for $$\max|\operatorname{bin}_n|$$. In the results obtained (Fig. 9) $$\operatorname{E}[|\operatorname{bin}_n|^2]$$ is about -157.1 dB and $$\max|\operatorname{bin}_n|$$ about -147.2 to -146.0 dBFS for non-fundamental harmonics. For the fundamental frequency, $$\operatorname{E}[|\operatorname{bin}_1|^2]$$ is in the order of -0.95E-11 dBFS, and $$\max|\operatorname{bin}_1|$$ is in the order of 3.4E-7 dBFS or higher. Figure 9. Estimated $$\operatorname{E}[|\operatorname{bin}_n|^2]$$ (blue squares) and $$\max|\operatorname{bin}_n|$$ (red lines) for local $$\operatorname{bin}_n$$ in slow-drifting asynchronous sampling, for $$A = 2^{23}-1$$, $$f =$$ 1000 Hz, and $$f_s =$$ 48 kHz, estimated by a uniform statistical sample of phase shifts with a sample size of $$2^{18}$$.

The estimates show that asynchronous sampling gives a smoothing effect on the statistics of the local amplitudes of the harmonics, compared to synchronous sampling with a fixed alignment of the prototypical cosine with respect to the sampling grid. However, with asynchronous sampling, sometimes the worst case local harmonic peaks occur, whereas with synchronous sampling the harmonic amplitudes are well controlled and the worst cases are avoided by a suitable choice of phase shift.

# Asynchronous sampling of a sinusoid with rational $$f/f_s$$ in the limit $$L/A\to 0$$

In this section of the answer, we will consider only the question's case that $$L$$ is even and $$f/f_s = 1/L$$, that is to say that $$\operatorname{bin}_1$$ corresponds to the frequency of the prototypical sinusoid.

If the length of the fundamental period of the locally periodic sequence of samples $$L$$ calculated by Eq. 4 is very small compared to the amplitude $$A$$ of the prototypical sinusoid, then for different phase shifts of the sinusoid, the quantization error can perhaps be considered a sequence of independent pseudorandom numbers distributed uniformly from -0.5 to 0.5, with an antiperiodicity constraint if $$L$$ is even (and some other constraint if $$L$$ is odd). Independence of the quantization errors would require that at a number of certain small time lags the autocorrelation of the quantized prototypical sinusoid approaches zero, but this is not something we know to be a fact. If it was true, then in the limit $$L/A \to 0$$ any sequence consisting of ±0.5's could be found as a quantization error sequence at some phase shift $$\phi$$ of the prototypical cosine, and within those sequences there would be sequences that maximize $$|\operatorname{bin}_n|$$ of the real DFT of the quantization error sequence. Due to orthogonality of the real DFT basis vectors, if $$n$$ does not correspond to the fundamental frequency of the sinusoid, then $$\operatorname{bin}_n$$ of the quantized sequence equals $$\operatorname{bin}_n$$ of the quantization error. If $$n$$ corresponds to the fundamental frequency, then $$\operatorname{bin}_n$$ of the quantized sequence equals $$\operatorname{bin}_n$$ of the quantization error plus $$Ae^{i\phi}$$.

Even if we do not assume independence of quantization errors, we can use the line of thinking of the previous paragraph to find an upper bound of $$\max|\operatorname{bin}_n|$$, the maximum absolute bin value for a harmonic of the prototypical sinusoid. A worst-case sequence of ±0.5's that maximizes $$|\operatorname{bin}_n|$$ could be found by an exhaustive search in the binary sequences, but can also obtained directly. A worst-case sequence is any sequence $$\operatorname{seq}_k = 0.5\operatorname{sgn}\big(\cos(2\pi nk/L + \chi)\big)$$, where $$\operatorname{sgn}$$ is the sign function, and with a phase shift $$\chi$$ that does not result in a zero crossing of the expression's sinusoid being at any sample point (the integers). One worst-case sequence is found in the limit of approaching from one side the transition between two worst-case sequences:

$$\begin{eqnarray}\operatorname{seq}_k &=& 0.5\times\begin{cases}\operatorname{sgn}\big(\sin(2\pi nk/L)\big),&\text{if }\sin(2\pi nk/L) \ne 0,\\ \cos(2\pi nk/L),&\text{if }\sin(2\pi nk/L) = 0,\end{cases}\\ &=& 0.5\times\begin{cases}\operatorname{sgn}\big(\sin(2\pi nk/L)\big),&\text{if }2kn \text{ mod } L \ne 0,\\ \cos(2\pi nk/L),&\text{if }2kn \text{ mod } L = 0,\end{cases}\end{eqnarray}\tag{7}$$

where $$\text{ mod }$$ denotes remainder in integer division in a numerically safe expression. An upper bound of $$\max|\operatorname{bin}_n|$$ is then found by real DFT of the worst-case sequence.

## 24-bit quantization

For odd $$n$$ the sequence given by Eq. 7 satisfies the antiperiodicity constraint for even $$L$$, and we can calculate an upper bound of $$\max|\operatorname{bin}_n|$$ of the maxima of local harmonic amplitudes $$|\operatorname{bin}_n|$$ for $$A = 2^{23}-1$$, $$f =$$ 1000 Hz and $$f_s =$$ 48 kHz, in Python:

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")


For the odd harmonics, for $$A = 2^{23}-1$$, $$f =$$ 1000 Hz and $$f_s =$$ 48 kHz, we get locally $$\max|\operatorname{bin}_n| \le$$ -142.34 dBFS when the number of the harmonic $$n$$ is divisible by 3, and $$\max|\operatorname{bin}_n| \le$$ -142.39 dBFS otherwise which applies to the 5th and to the 7th harmonic. This seems compatible with Fig. 9, with an about 4 dB higher upper bound of the maximum local amplitudes of harmonics compared to the maximum local harmonic amplitude found by numerical search in actual quantization error sequences. A willing eye can also see an elevation in $$\max|\operatorname{bin}_n|$$ of Fig. 9 when $$n$$ is divisible by 3.

This answer has reached the maximum allowed length of 30000 characters.

• Those alias-reflected harmonics are so important. It explains why small phase changes make huge differences to the measured harmonic values in practice. I know that I am not supposed to express approval; but all praise to you Olli. I wish I had half your mental agility. Mar 21 '19 at 7:15
• @Richard, You are most definitely supposed to express approval. That's what the little up arrows are for in voting for an answer which can be done by anyone. The ultimate approval is acceptance of an answer which is done by the OP, that's you in this case. Mar 21 '19 at 14:23
• i don't think this answer is sufficiently exhaustive .  :-) Mar 25 '19 at 20:57
• @robertbristow-johnson that's true... :D hold on. Mar 26 '19 at 5:40
• you know what they said when lasers were invented: "a very good solution to an unknown problem" ;-)) Mar 27 '19 at 22:23

The harmonics will be in the range of -150dB to -170dB. The exact value will depend on exact frequency, phase relationship to sampling frequency, phase locked or unlocked, rate of phase drift, integration interval, index of harmonics, etc.

You need to specify this all in excruciating detail to get anywhere near the resolution that you are asking for.

I also seem rather pointless: it's very hard to imagine a 24-bit system where the actual noise floor is dominated by the 24-bit quantization noise. Even if it were: why would you need to now this with that extreme level of detail. It's hard to imagine an application where a harmonic at -154.02 dB would be meaningfully different from a harmonic at -154.01 dB

• even though the sampling is "asychronous", the ratio of $f_\mathrm{s}$ to $f_0$ is an integer (48). how will the amplitude at the frequencies of the harmonics be any different than that of the other frequencies that are not $f_0$? since it's 24-bits, i would expect something like -144 dB everywhere. Mar 20 '19 at 19:45
• With 24 bits, the total quantization noise power should be $(2^{-23})^2/12$, which is around -149 dB (for rounding). That power gets spread out over all the harmonics, so each individual harmonic must be less than that. If sampling is assynchronous, the ratio is not an integer, but a number that's very close to an integer, such 48.0001 for a 100ppm clock. As a result the phase will very slowly drift and hence move the energy around between the different harmonics Mar 20 '19 at 21:52
• you're right about the -4.77 dB offset. otherwise it's -6.02 dB per bit. for async, that power should be spread out over not just the harmonic frequencies, but all frequencies below Nyquist nearly equally. for synchronous sampling with an integer ratio, it's a perfectly periodic function and that same noise power should exist only at the harmonic frequencies. but if 149 dB is a large number, it the rounding, even if synchronous with other cycles, still should be uniformly random for individual samples. it just won't be broadbanded quantization noise. Mar 21 '19 at 1:33
• I think you are right: if it's "async" the quantization spectrum will be white. If it's "sync" you get a harmonic spectrum and you get something in between for "almost sync" depending on the amount of clock difference/drift AND the window of observation. Mar 21 '19 at 12:40

Welcome to dsp.stackexchange.

I was kind of curious, so here goes.

If I am understanding correctly, you are trying to tease apart harmonic distortion caused by your ADC vs quantization error. That is, get an idea of the size of each. Luckily, the latter can be measured numerically. It's not theoretical, it's numerical, but sometimes that's just as good for practical purposes.

I've whipped up some code, shown below, in a language called Gambas, which I recently found and am already very fond of. Unfortunately, it doesn't have the DSP support in it (yet) that would allow me to make a broad recommendation for DSP processing. For folks who like VB6, you'll want to install it ASAP.

PPA: gambas-team/gambas3

Notes: Float is a double precision value, Floor(-0.5) = -1, etc.

The calculations are the Fundamental, Third, Fifth, and Seventh Harmonics 1/N normalized DFT bins then divided by the amplitude.

The results are:

8388607

0.50000000212064
-3.97347889201056E-11

3.58733863253808E-10
3.97364158530641E-11

-6.2388133833524E-10
-3.97364423993961E-11

-2.44162068960777E-9
3.97363696177629E-11


Error free would be 0.5, then zeros. You can 0.01 decibelize them if you want. Different parameters will give slightly different results. This is for one second of your data.

Since the sampling rate is commensurate with the frequencies, you should expect the quantization errors to be cumulative, so you should divide the resulting errors by the cycles per frame to get your per cycle values.

I used the cosine fot the signal. It is interesting how the imaginary bins have very similar values for all the harmonics.

Floats are good to about 10^(-16), so only the first 4 or 5 digits are worth anything.

Hope this helps.

Ced

' Gambas module file

'=============================================================================
Public Sub Main()

Dim n, theSamplesPerFrame, theCyclesPerFrame As Integer
Dim theAmplitude, theNormalizer As Float

Dim theRadiansPerSample, theThirdRps, theFifthRps, theSeventhRps As Float

Dim theRealSum, theThirdRealSum, theFifthRealSum, theSeventhRealSum As Float
Dim theImagSum, theThirdImagSum, theFifthImagSum, theSeventhImagSum As Float

Dim theSignalValue, theSignalQuant As Float

theAmplitude = 2 ^ 23 - 1

theSamplesPerFrame = 48000
theCyclesPerFrame = 48

theNormalizer = 1 / theSamplesPerFrame

theRadiansPerSample = Pi(2) * theCyclesPerFrame / theSamplesPerFrame

theRealSum = 0.0
theImagSum = 0.0

theThirdRealSum = 0.0
theThirdImagSum = 0.0

theFifthRealSum = 0.0
theFifthImagSum = 0.0

theSeventhRealSum = 0.0
theSeventhImagSum = 0.0

For n = 0 To theSamplesPerFrame - 1
theSignalValue = theAmplitude * Cos(theRadiansPerSample * n)
theSignalQuant = Floor(theSignalValue)

theRealSum += theSignalQuant * Cos(theRadiansPerSample * n)
theImagSum += theSignalQuant * Sin(theRadiansPerSample * n)

theThirdRealSum += theSignalQuant * Cos(theThirdRps * n)
theThirdImagSum += theSignalQuant * Sin(theThirdRps * n)

theFifthRealSum += theSignalQuant * Cos(theFifthRps * n)
theFifthImagSum += theSignalQuant * Sin(theFifthRps * n)

theSeventhRealSum += theSignalQuant * Cos(theSeventhRps * n)
theSeventhImagSum += theSignalQuant * Sin(theSeventhRps * n)
Next

Print theAmplitude
Print
Print theNormalizer * theRealSum / theAmplitude
Print theNormalizer * theImagSum / theAmplitude
Print
Print theNormalizer * theThirdRealSum / theAmplitude
Print theNormalizer * theThirdImagSum / theAmplitude
Print
Print theNormalizer * theFifthRealSum / theAmplitude
Print theNormalizer * theFifthImagSum / theAmplitude
Print
Print theNormalizer * theSeventhRealSum / theAmplitude
Print theNormalizer * theSeventhImagSum / theAmplitude

End
'=============================================================================


Here are some more numeric results.

I was curious about your "why small phase changes make huge differences to the measured harmonic values in practice". So, I upped my code (not shown) and ran a bunch of trials with random phases, then measured the errors.

These are the ( Magnitude - Expected Magnitude ) values:

Errors:
Minimum                Maximum                 Average              Standard Deviation

Trials = 10

1st   -2.24787882752153E-9     1.98062277689814E-9    -1.28245553066009E-10    1.44083332364943E-9
3rd    5.84740316427721E-10    2.25501205370262E-9     1.30720619887375E-9     5.22116141270907E-10
5th    3.12618630073034E-10    3.50159759485387E-9     1.40715363571154E-9     1.11117800026694E-9
7th    3.00033566816365E-10    2.65474965884586E-9     1.41006537724231E-9     6.90489229891695E-10

Trials = 100

1st   -2.64673338801202E-9     2.49246612238352E-9     3.92873139576366E-11    1.03351131455717E-9
3rd    1.25076964468878E-10    3.54534685452094E-9     1.27463466771551E-9     6.71882224324839E-10
5th    1.54350057280304E-10    3.36439759301702E-9     1.46553038670545E-9     7.38531583568459E-10
7th    2.43528172304884E-10    3.73979952276669E-9     1.29326190868848E-9     7.09356107751877E-10

Trials = 1000

1st   -3.67113961452503E-9     2.97127311643663E-9     6.26527135727528E-12    1.05723843099991E-9
3rd    2.34467663801663E-11    4.5534373569163E-9      1.38396033029499E-9     7.26924843417221E-10
5th    1.65873381161482E-11    4.07156873679203E-9     1.32658045132769E-9     6.854053318804E-10
7th    4.59417305079006E-11    3.99110453111402E-9     1.34706083510439E-9     7.13102512232999E-10



Olli's treatments are works of art.