Return to Answer

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

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 $N=2^{23} - 1$, 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, the amplitude of the 1st harmonic is, interestingly, 0.0000001605530 dBFS, 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.

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.

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:

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

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, which is in the order of 4E-11 dBFS.

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

16-bit

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

24-bit

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

16-bit quantization

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

24-bit quantization

Then back to my calculations that concern the limit $f/f_s \to 0$ in which case there is no aliasing. I will try to compute overnight your desired numbers for $N=2^{23} - 1$ in high precision, in the limiting case $f/f_s \to 0$, using the following:

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. Here, the amplitudes decay extremely slowly as function of harmonic number.

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

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.