Calculating THD+N which one is the proper way?

[hiban]

So i have this personal software project, and i need to calculate the THD+N of a waveform starting from the FFT.
Some questions arise about how to exactly calculate the THD+N. Not about the formula, but about what exact data to pick for the calculations.

My FFT is stored in an array of Y values, so it's later plotted on a screen so the user can look at it. I'm using zero-padding to obtain a "higher resolution" FFT representation, in which each bin is composed of several sample points instead of only 1 sample point, this is so i can look at it in "higher detail" (although the number of samples from the input waveform is still the same, but it helps to better distinguish some details). I tell all this so there is some context.

Ok, so the FFT in that array has its lobes or bins, with the fundamental and harmonic bins being 2x the width of the rest of the lobes/bins. So when i had to calculate the THD+N, a question arised:
Do i need to handle the wider bin (fundamental or harmonic) as a single bin and use the highest value (center) as the "fundamental value"? or do i have to handle it like if it was 2 separate bins (because it's 2x the width and equals to 2 samples) and use the sum of the center samples of those imaginary "separate bins"? Of course, this would apply to both the fundamental and the harmonic bins, which are also 2x the width.

I have attached a picture for explanation. The FFT is the red representation at the bottom, and as you can see, the fundamental bin is 2x the width, so i can think about 2 options:
Option 1: handle it as a single bin, and use the value of the point marked in yellow for the fundamental.
Option 2: Since the fundamental is 2x as wide, handle it as 2 separate bins and use the sum of the 2 values marked in blue as the fundamental value.


For the rest of the bins, i calculate the range of samples covered by every bin and use the highest value point within that range.

As you can guess, the resulting THD+N value will be different depending on the option i choose, and each different option requires me to write different code in order to make it work, so how is it properly done? Thanks in advance.

P.S: ignore the values from the screenshot, they are placeholders.

[radiolistener]

The width of the fundamental peak in an FFT is influenced not only by the signal itself but also significantly by the window function applied before performing the FFT. Different window functions were designed specifically to balance trade-offs between frequency resolution and amplitude accuracy.

A window with a narrow main lobe (like a rectangular window) provides better frequency resolution but results in more spectral leakage (i.e., high sidelobes). Conversely, a window with better sidelobe suppression (like Nuttal, Blackman-Harris or Hann) sacrifices frequency resolution for dynamic range and cleaner separation between spectral components.

Regarding how to handle the wider fundamental or harmonic peaks, a common and recommended approach is to identify the full extent of the main lobe, which typically includes all bins from the peak down to where the magnitude stops decreasing (or starts increasing again). This region is considered part of the fundamental (or harmonic) tone.

Then integrate the energy across all these bins to compute the power of the tone (not just the center bin).

For THD+N, exclude these bins from the total spectrum when calculating the noise+distortion power. Usually, these bins are replaced or masked with an estimate of the average noise floor to avoid bias.

This method gives you a more accurate measure of both the fundamental and the remaining noise+harmonic content, especially when spectral leakage and windowing effects are significant.

So between your two options, neither Option 1 nor 2 is fully correct unless you're accounting for the full energy spread across the lobe. You should:

1) Not rely on a single bin.

2) Instead, sum the power (not magnitude) of all bins that are part of the main lobe.

Regarding the choice of window function, I personally prefer the Kaiser window, as it offers the flexibility to adjust the analysis parameters through the beta parameter. This allows you to directly control the level of sidelobe suppression and thereby optimize the trade-off between frequency resolution and amplitude dynamic range according to the specific requirements of your application. This makes the Kaiser window a versatile choice for both spectral visualization and precise measurements such as THD+N.

[Tation]

Comprehensive catalog of windows. Sure the Taylor window on page 132 will be of interest.

[SiliconWizard]

My usual go-to window for non-specific spectral analysis is Blackman-Harris (7-th order) - most often a good compromise.

[radiolistener]

My usual go-to window for non-specific spectral analysis is Blackman-Harris (7-th order) - most often a good compromise.

Yes, it is good. Blackman-Harris-7 often is better than Nuttall, but not always.

[hiban]

The width of the fundamental peak in an FFT is influenced not only by the signal itself but also significantly by the window function applied before performing the FFT. Different window functions were designed specifically to balance trade-offs between frequency resolution and amplitude accuracy.

A window with a narrow main lobe (like a rectangular window) provides better frequency resolution but results in more spectral leakage (i.e., high sidelobes). Conversely, a window with better sidelobe suppression (like Nuttal, Blackman-Harris or Hann) sacrifices frequency resolution for dynamic range and cleaner separation between spectral components.

Regarding how to handle the wider fundamental or harmonic peaks, a common and recommended approach is to identify the full extent of the main lobe, which typically includes all bins from the peak down to where the magnitude stops decreasing (or starts increasing again). This region is considered part of the fundamental (or harmonic) tone.

Then integrate the energy across all these bins to compute the power of the tone (not just the center bin).

For THD+N, exclude these bins from the total spectrum when calculating the noise+distortion power. Usually, these bins are replaced or masked with an estimate of the average noise floor to avoid bias.

This method gives you a more accurate measure of both the fundamental and the remaining noise+harmonic content, especially when spectral leakage and windowing effects are significant.

So between your two options, neither Option 1 nor 2 is fully correct unless you're accounting for the full energy spread across the lobe. You should:

1) Not rely on a single bin.

2) Instead, sum the power (not magnitude) of all bins that are part of the main lobe.

Regarding the choice of window function, I personally prefer the Kaiser window, as it offers the flexibility to adjust the analysis parameters through the beta parameter. This allows you to directly control the level of sidelobe suppression and thereby optimize the trade-off between frequency resolution and amplitude dynamic range according to the specific requirements of your application. This makes the Kaiser window a versatile choice for both spectral visualization and precise measurements such as THD+N.

Thanks for the reply. It will be very useful for improving my calculations.

You said that we add the averaged noise to the empty spaces left by the fundamental and the harmonics in order to account for the noise floor, but should we also subtract the averaged noise from the fundamental and harmonics when calculating those (to account for the existing noise  in the fundamental and harmonics)??

Also, i noticed that, when using most of the different windows, i get higher-than-zero (1-6% depending on the window) THD and THD+N values when testing with a pure sine wave (generated from the app itself using a sine function, so there is zero possible noise or distortion). And i noticed that, if i want to get the expected measurements (zero-point-zero THD/THD+N using a pure sine wave), i must choose a window which attenuates everything except the main fundamental lobe down to at least -85dB. So far, i've had best results with Nuttall and Blackman-Harris.
Flattop also works, but the main lobe gets way too wide, so i guess that could cause non-accurate measurements in certain situations, if it overlaps with the harmonics.

[radiolistener]

You said that we add the averaged noise to the empty spaces left by the fundamental and the harmonics in order to account for the noise floor, but should we also subtract the averaged noise from the fundamental and harmonics when calculating those (to account for the existing noise  in the fundamental and harmonics)??

That's an interesting question. However, I don't believe subtracting the averaged noise floor from the fundamental and harmonics would provide any significant benefit, as the power of the fundamental and its harmonics is much greater than that of the noise floor.

Also, i noticed that, when using most of the different windows, i get higher-than-zero (1-6% depending on the window) THD and THD+N values when testing with a pure sine wave (generated from the app itself using a sine function, so there is zero possible noise or distortion).

If you're obtaining such large results for a pure sine wave generated by the software, it suggests there may be an error in your calculations, as THD/THD+N cannot exceed 100% by definition.

Here is the test MATLAB code that demonstrates an exceptionally low THD+N value for a pure sine wave.

Code: [Select]

Hn = 10;        % harmonics number used in the calculations (including the fundamental)

Fs = 384000;    % sampel rate
N  = 2^22;      % sampel count
signal = cos(997 * 2*pi*(0:N-1)/Fs);    % test sine 997 Hz (the AES17-1998(r2004) standard recommends a 997 Hz frequency sine wave)

% normalize & adjust amplitude to 1 Vrms (just for better readability of power values)
signal = signal * sqrt(2) / max(abs(signal));

% Calculate SNR
[r,nois_pow] = snr(signal, Fs, Hn);
fprintf('SNR:   %.4f dB\n', r);

% Calculate SFDR
[r,spur_pow,spur_freq] = sfdr(signal, Fs);
fprintf('SFDR:  %.4f dB\n', r);

% Calculate SINAD
[r,totdist_pow] = sinad(signal,Fs);
fprintf('SINAD: %.4f dB\n', r);

% Calculate THD
[r,harm_pow,harm_freq] = thd(signal, Fs, Hn);
fprintf('THD:   %+.4f dBc = %g %%\n', r, 100*sqrt(10^(r/10)));


% Calculate THD+N
fund_rms = sqrt(10^(harm_pow(1)/10));       % RMS voltage of the fundamental
harm_rms = sqrt(10.^(harm_pow(2:end)/10));  % RMS voltage of the harmonics
nois_rms = sqrt(10^(nois_pow/10));          % RMS voltage of the noise
r = 20*log10(sqrt(sum(harm_rms.^2) + nois_rms^2) / fund_rms); % THD+N in dBc
fprintf('THD+N: %+.4f dBc = %g %%\n', r, 100*sqrt(10^(r/10)));
 
% Some more details
fprintf('Fundamental Freq:  %g Hz\n',        harm_freq(1));
fprintf('Fundamental Power: %+.4f dBW = %g Vrms\n', harm_pow(1), fund_rms);
fprintf('Noise Power:       %+.4f dBW = %g Vrms\n', nois_pow,    nois_rms);
fprintf('Spur %+.4f dBW = %+.4f dBc at %g Hz\n', spur_pow, spur_pow-harm_pow(1), spur_freq);

Results:

Code: [Select]

SNR:   232.7518 dB
SFDR:  242.0736 dB
SINAD: 232.7518 dB
THD:   -298.6906 dBc = 1.1627e-13 %
THD+N: -232.7518 dBc = 2.30361e-10 %
Fundamental Freq:  997 Hz
Fundamental Power: +0.0000 dBW = 1 Vrms
Noise Power:       -232.7518 dBW = 2.30361e-12 Vrms
Spur -242.0736 dBW = -242.0736 dBc at 135545 Hz

As observed, the THD+N value for the pure sine wave is 2.30361e-10 %, which is near the floating-point resolution limit.

[SiliconWizard]

1% or more for the calculated THD+N on a pure sine wave, unless you use a horrific sin() function to generate it, or use very limited precision, is definitely not normal.

First thing to check is how you determine the magnitude of the frequency components from the FFT.
As radiolistener mentioned earlier, you should not use a single FFT bin (such as the closest to your target frequency). This isn't correct and won't work unless your fundamental frequency is spot at the center of the bin, and even so - sidelobes will always bite you even with the "best" window.

So instead of using single bins, the usual approach (that I've used successfully) is to compute the magnitude for a given frequency from the FFT bin modules using the equivalent of a bandpass filter centered around the frequency of interest. The width of this filter should be proportional to the frequency of interest and is usually set as a fraction of octave, but it can be defined as a fraction of decade, doesn't matter.

[Tation]

To the OP, see pages 10-17 of this application note from NI. Discusses exactly this.

[hiban]

I changed the way i calculated the magnitudes. So now instead of selecting a single bin from the center of the sidelobes, i now just pick all bins from the desired range, as you guys suggested. So for example, when using a window, i assign the width of the main lobe as the range width of frequencies for both the fundamental and the harmonics. And in case they overlap (example: flattop window), i make sure to handle those exceptions so the bins are never counted twice, and i also make sure that the fundamental always takes priority over the harmonics.

Anyway, for me the problem became obvious when i looked at the FFT representations of the pure sine wave using different windows:


I noticed that a magnitude of about -87 dB (this is in reference to the signal peak value, so i guess the units should be dBc?) Anywayy, i noticed at -87 dB the magnitude is too small to even be noticeable when calculating THD+N with 2 decimals. The magnitude is so small you will not even get 0.01% in the final result. But if the magnitudes are higher than -87dB, then they start being big enough to sum at least 0.01% in the final result. Still super small, but they start being noticeable for the final result. So in the following images i placed the Y1 cursor at about -87dB, for reference.


If you take a look at image number 1, that is the FFT representation of a pure sine wave, using a Rectangular window (no-window). So all the FFT section marked in green is the harmonics+noise higher than -87dB. In fact, the highest sidelobe is -12.16 dB. All that portion of the signal is counted as THD+N, and that's the reason the calculated THD/THD+N values are so high.


Image number 2, same pure sine wave, now using a Hann window. As you can see, the magnitudes of the part of the FFT counted as THD/THD+N are now much smaller, and now the highest sidelobe is only -31.41 dB. Much smaller magnitudes, but there is still a portion over -87dB, which makes the zone in green to be noticeable for the THD/THD+N calculation. The THD/THD+N results are much smaller, but still high for a pure sine wave, which should be zero.


Image number 3, same pure sine wave, using a Blackman window. Highest sidelobe is now way lower, only -56.85 dB. However, the part of the FFT higher than -87dB is still big enough to be noticeable in the THD/THD+N calculations. The calculation results are now near zero, but still not zero.


Image number 4, same pure sine wave, using a Blackman-Harris window (Other windows like Nuttall also give similar results as Blackman-Harris). Now the highest sidelobe is -90 dB. So all the non-fundamental FFT magnitudes are now so small that they don't make even 0.01% in the THD/THD+N calculations. So the results are now 0.00%, as expected.


And to check if it's working properly, i added just a bit of noise to the pure sine wave, and i was welcomed with a proper randomized and everchanging noise floor, with the consequent THD/THD+N results reflecting the presence of noise. Keep in mind that i am not (yet) accounting for the noise in the fundamental or harmonic calculations, so if there is noise contained in the harmonic frequency ranges, this noise is still being counted as "harmonic distortion". I will modify the code to account for that aswell, but for now i was only focusing on the issue i explained in this post.

Ok, so after explaining all of my experiments, my question is: Is it ok to use a Blackman-Harris window? I ask because my understanding prior to this was that the Hann window was widely used for calculating THD/THD+N, but they must have some compensation algorithm implemented in their formulas, because the magnitudes of the FFT are always there unless you use a window with big suppression, like Blackman-Harris or Nuttall.

The problem with using a window with a wide main lobe like Blackman-Harris is that, if the main lobe frequency range (width) is higher than the frequency of the fundamental, then the fundamental and harmonic ranges will overlap. And this means that all existing bins will be used to calculate the fundamental and harmonics, leaving no bins left for calculating the noise.
And if i can't calculate the noise, then i cannot account for it, because it will always be calculated as "zero".
So if the fundamental frequency is 1KHz, then i need a window with a main lobe frequency range lower than 1KHz. Because a main lobe with a width of 1KHz will extend to -0.5KHz and +0.5KHz from the fundamental frequency, and the same will happen with the harmonic lobes, so for the first harmonic (2KHz) both the fundamental and harmonic lobes will meet at 1.5KHz, and this will extend for all the harmonics, leaving no bins left to calculate the noise. 

[radiolistener]

When referring to the Blackman-Harris window, it’s important to clarify that there are two common variants: Blackman-Harris-4 and Blackman-Harris-7. They are very different! Many sources and tools refer to the 4-term variant when simply saying “Blackman-Harris” but in your case, Blackman-Harris-7 is the more appropriate choice. It provides significantly better side lobe suppression, resulting in much cleaner THD/THD+N calculations with very low distortion leakage - far outperforming Blackman-Harris-4 in magnitude resolution.

The Nuttall window performs similarly because it is essentially a slight modification of the Blackman-Harris-7 design. Both are well-suited for high dynamic range analysis like THD/THD+N.

However, for the best results, especially when ultra-low THD+N measurements are needed, I would recommend using the Kaiser window with an appropriate β value. For example, MATLAB’s built-in thd function uses a Kaiser window with β = 38, which corresponds to approximately 350 dB side lobe attenuation - an exceptionally high suppression that effectively eliminates leakage artifacts from the FFT bins surrounding the fundamental and harmonics.

So yes, it is not only ok to use Blackman-Harris-7 or Nuttall - it is recommended for this application. Hann, Blackman, and especially Rectangular windows have much poorer side lobe suppression, which leads to leakage that contaminates the noise and harmonic estimates. If you must use these windows, you would need to compensate by significantly increasing your FFT size (sample length) to improve magnitude resolution - which is not always practical or possible.

[radiolistener]

The problem with using a window with a wide main lobe like Blackman-Harris is that, if the main lobe frequency range (width) is higher than the frequency of the fundamental, then the fundamental and harmonic ranges will overlap. And this means that all existing bins will be used to calculate the fundamental and harmonics, leaving no bins left for calculating the noise.

Yes, wide main lobes in windows like Blackman-Harris can indeed reduce frequency resolution, potentially causing the main lobe of the fundamental or harmonics to overlap adjacent bins. This makes it difficult to distinguish between harmonic energy and noise, especially at low frequencies or with short FFT lengths.

This is precisely why adjustable windows, such as the Kaiser window, were developed. The Kaiser window includes a tunable parameter β (beta) that allows you to control the trade-off between main lobe width (frequency resolution) and side lobe suppression (magnitude resolution). You can tailor it to your specific application - whether your priority is dynamic range or resolving closely spaced frequencies.

To determine the optimal β for your THD/THD+N measurements, you can use this simple MATLAB/Octave function:

Code: [Select]

% A - required suppression in dB
function bta = kaiserBeta(A)
    if (A >= 50)
        bta = 0.1102 * (A - 8.7);
    elseif (A < 50 && A > 21)
        bta = 0.5842 * (A - 21)^0.4 + 0.07886 * (A - 21);
    else
        bta = 0;
    end
end

Here is example of beta calculated for different suppression:

Code: [Select]

kaiserBeta(40 dB) = 3.39532
kaiserBeta(60 dB) = 5.65326
kaiserBeta(80 dB) = 7.85726
kaiserBeta(100 dB) = 10.0613
kaiserBeta(120 dB) = 12.2653
kaiserBeta(140 dB) = 14.4693
kaiserBeta(160 dB) = 16.6733
kaiserBeta(180 dB) = 18.8773
kaiserBeta(200 dB) = 21.0813
kaiserBeta(220 dB) = 23.2853
kaiserBeta(240 dB) = 25.4893
kaiserBeta(260 dB) = 27.6933
kaiserBeta(280 dB) = 29.8973
kaiserBeta(300 dB) = 32.1013
kaiserBeta(320 dB) = 34.3053
kaiserBeta(340 dB) = 36.5093
kaiserBeta(360 dB) = 38.7133
kaiserBeta(380 dB) = 40.9173
kaiserBeta(400 dB) = 43.1213

This way, you can simply specify the minimum required sidelobe suppression in dB, and the Kaiser window will provide the narrowest possible main lobe that still meets that suppression level.

[hiban]

In my case i'm using a Blackman-Harris-4.
I tried the Kaiser window with different beta values, and found out that using beta = 12 provided enough side lobe attenuation for obtaining 0.00% THD/THD+N in the calculations for a pure sine wave.

However, i also observed that, the higher the sidelobe attenuation, the wider the main lobe will be. So at the end of the day, a Kaiser window with a sidelobe attenuation roughly equal to a Blackman-Harris-4 window, will have a main lobe with roughly the same width as a Blackman-Harris-4 window. And if i increase the beta value of the Kaiser window, i will obtain a higher sidelobe attenuation, but at the cost of a wider main lobe.

Same happens with the Nuttall window, it has roughly the same sidelobe attenuation , and also roughly the same main lobe width. And i have the problem of the main lobe width potentially overlapping with harmonics, which is indeed the case with any window which attenuates the sidelobes as much as i need.

In my application i have some limitations:

The app i'm creating is completely written in plain javascript with no external libraries, because it's a web app and all the code is contained in a single HTML document. So i use HTML, CSS, and javascript in order to create all the graphics and make everything work from a web browser like Chrome or Microsoft Edge. So everything is being coded from scratch: FFT functions, interpolation algorithms, waveform measurement calculations... everything.

I started making this app because i was tinkering with a cheap "toy" oscilloscope i have, the "zeewei DSO2512G". A user created a modified firmware which added the capability of communicating with a PC via USB serial, so i decided to try to do a simple command console using HTML and javascript. Since the modded firmware also had some commands to request for waveform data, then i tried to request for waveform data and plot a waveform, and then i decided i could also compute the FFT from the waveform data, and then i started adding more and more features from real oscilloscopes like stacked mode, cursors (manual and track modes), a math channel, a configurable low pass input filter, averaging aquisition mode, waveform interpolation, and so on. So after a month i ended up with the abomination you can see in the attached picture.

However, i am limited by the quality of the data sent by the DSO2512G (among other things). So i must work with waveform data arrays of 1200 samples at most. That's the maximum resolution i can obtain on a raw waveform array. And the sample count can go down a lot when visualizing high frequency signals, because the DSO2512G max sample rate is 200MSa/s (only 100MSa/s when channel 2 is enabled). Of course i can interpolate the waveform data, and i do for visualization, but for performance reasons i decided to keep the data arrays at 1200 samples per waveform. Keep in mind this is not a precompiled executable, it's plain javascript code being executed in a browser, so it has its performance limitations. And for performance reasons, i also downsample the FFT input waveform used to calculate the THD/THD+N to 1024 samples, so i can compute the FFT with faster algorithms which require the input number of samples to be a power of 2.

I have thought about interpolating the input signal to higher power-of-2 values, like 2048. But i don't really know if increasing the input sample length will help with either spectral leakage or window main lobe width.
EDIT: I tried interpolating to higher values, but apparently it didn't help.

[Tation]

Apply the window to your 1200 points time-domain waveform, zero-pad to 2048 points and apply FFT, thus obtaining an interpolated spectrum.

[hiban]

Apply the window to your 1200 points time-domain waveform, zero-pad to 2048 points and apply FFT, thus obtaining an interpolated spectrum.

That's what i already do for computing the FFT i render to the screen. I pad with zeros (up to a higher power of 2) to achieve an horizontal zoom effect which the user can set. But all i get by doing that is a higher-horizontal-resolution FFT. The main lobe is still the same frequency width, all the other lobes have the same size and shape, and the FFT is exactly the same, regardless of the detail level of the FFT output or the number of samples of the input waveform.
So at the end of the day, the waveform or FFT resolution does not seem to be the issue here.

Of course if the resolution is way too low, then you start losing precision and problems arise. But once you got enough resolution (and 1200 samples in the input waveform is enough), adding more does not really improve anything, at least not if you're displaying your measurements with only 2 decimals.

[radiolistener]

In my case i'm using a Blackman-Harris-4.
I tried the Kaiser window with different beta values, and found out that using beta = 12 provided enough side lobe attenuation for obtaining 0.00% THD/THD+N in the calculations for a pure sine wave.

However, i also observed that, the higher the sidelobe attenuation, the wider the main lobe will be. So at the end of the day, a Kaiser window with a sidelobe attenuation roughly equal to a Blackman-Harris-4 window, will have a main lobe with roughly the same width as a Blackman-Harris-4 window. And if i increase the beta value of the Kaiser window, i will obtain a higher sidelobe attenuation, but at the cost of a wider main lobe.

Same happens with the Nuttall window, it has roughly the same sidelobe attenuation , and also roughly the same main lobe width. And i have the problem of the main lobe width potentially overlapping with harmonics, which is indeed the case with any window which attenuates the sidelobes as much as i need.

What you're encountering is not a limitation of a particular window function like Kaiser, Nuttall, or Blackman-Harris - it's a fundamental trade-off inherent to all window functions.

You can either achieve high frequency resolution (narrow main lobe) or high dynamic range (low sidelobe level), but not both at the same time. Each window function represents a different balance between these two competing objectives.

The advantage of the Kaiser window is that it lets you adjust this balance manually through the β parameter. But no window, including Kaiser, can overcome the fundamental limitation: if the required suppression demands a wider main lobe than your signal can tolerate (e.g., causing overlap between fundamental and harmonics), then the only solution is to increase the FFT size (i.e., use more data) to improve frequency resolution.

In short, if a Kaiser window with the minimum required suppression doesn’t provide a narrow enough lobe for your application, then no other window will either - the limitation lies in the not enough amount of input data, not the choice of window.

Apply the window to your 1200 points time-domain waveform, zero-pad to 2048 points and apply FFT, thus obtaining an interpolated spectrum.

Zero-padding cannot overcome the fundamental limitation caused by insufficient information in the input signal. The only effective way to improve frequency resolution is to acquire a longer segment of actual signal data - not artificial extension with zeros.

That's what i already do for computing the FFT i render to the screen. I pad with zeros (up to a higher power of 2) to achieve an horizontal zoom effect which the user can set. But all i get by doing that is a higher-horizontal-resolution FFT. The main lobe is still the same frequency width, all the other lobes have the same size and shape, and the FFT is exactly the same, regardless of the detail level of the FFT output or the number of samples of the input waveform.

You're right regarding the visual effect - zero-padding increases the apparent frequency resolution of the FFT display but does not improve the actual frequency resolution or reveal more information about the signal.

However, be cautious: while zero-padding is fine for visualization, using zero-padded data in THD+N calculations can distort the spectral content and lead to inaccurate or misleading results. For precise measurements, it's essential to rely solely on actual signal data.

[hiban]

In my case i'm using a Blackman-Harris-4.
I tried the Kaiser window with different beta values, and found out that using beta = 12 provided enough side lobe attenuation for obtaining 0.00% THD/THD+N in the calculations for a pure sine wave.

However, i also observed that, the higher the sidelobe attenuation, the wider the main lobe will be. So at the end of the day, a Kaiser window with a sidelobe attenuation roughly equal to a Blackman-Harris-4 window, will have a main lobe with roughly the same width as a Blackman-Harris-4 window. And if i increase the beta value of the Kaiser window, i will obtain a higher sidelobe attenuation, but at the cost of a wider main lobe.

Same happens with the Nuttall window, it has roughly the same sidelobe attenuation , and also roughly the same main lobe width. And i have the problem of the main lobe width potentially overlapping with harmonics, which is indeed the case with any window which attenuates the sidelobes as much as i need.

What you're encountering is not a limitation of a particular window function like Kaiser, Nuttall, or Blackman-Harris - it's a fundamental trade-off inherent to all window functions.

You can either achieve high frequency resolution (narrow main lobe) or high dynamic range (low sidelobe level), but not both at the same time. Each window function represents a different balance between these two competing objectives.

The advantage of the Kaiser window is that it lets you adjust this balance manually through the β parameter. But no window, including Kaiser, can overcome the fundamental limitation: if the required suppression demands a wider main lobe than your signal can tolerate (e.g., causing overlap between fundamental and harmonics), then the only solution is to increase the FFT size (i.e., use more data) to improve frequency resolution.

In short, if a Kaiser window with the minimum required suppression doesn’t provide a narrow enough lobe for your application, then no other window will either - the limitation lies in the not enough amount of input data, not the choice of window.

Apply the window to your 1200 points time-domain waveform, zero-pad to 2048 points and apply FFT, thus obtaining an interpolated spectrum.

Zero-padding cannot overcome the fundamental limitation caused by insufficient information in the input signal. The only effective way to improve frequency resolution is to acquire a longer segment of actual signal data - not artificial extension with zeros.

I still have some tricks under the sleeve. Since the harmonic magnitudes decrease as frequency increases, i can estimate the average noise level by using the magnitudes of all the bins from 3/4 Nyquist up to Nyquist (in which any existing harmonics should be at their smallest level compared with the noise), and then subtract the estimated average from all the fundamental and harmonic bins. So that way we could assume that whatever magnitudes are left in the fundamental and harmonic bins, those magnitudes should contain the fundamental and harmonics... mostly. It's not 100.00% accurate, but should be around 99.9% accurate, even when using windows with overlapping harmonic lobes. Well, that's the theory at least.

[radiolistener]

If acquiring a longer real data sample is not possible, it's better to present the measurement results as is, based solely on the available data, rather than extending the signal with artificial content. Padding with synthetic data can introduce subtle distortions that may lead to significantly inaccurate results, which might not be immediately noticeable. In fact, substantial errors can emerge under specific signal conditions, making the measurement unreliable.

[hiban]

I'm trying to do what i can do with what i have. However, i'm sure other people must have encountered the same problem before. So how did they solve it?
Is there any way to look at the raw code from the "thd" matlab function? i'm curious to know how do they exactly compute the calculations.

[radiolistener]

Is there any way to look at the raw code from the "thd" matlab function? i'm curious to know how do they exactly compute the calculations.

just click on the function with right mouse button and open source code.

[Tation]

Apply the window to your 1200 points time-domain waveform, zero-pad to 2048 points and apply FFT, thus obtaining an interpolated spectrum.

Zero-padding cannot overcome the fundamental limitation caused by insufficient information in the input signal. The only effective way to improve frequency resolution is to acquire a longer segment of actual signal data - not artificial extension with zeros.

Having more samples of the same spectrum can be useful in some calculations, apart from being useful to display/zoom. Some time ago a method to compute precise estimations of the frequency of harmonics relied on zero padding to get enough bins near the peak, then fit a parabole thru the highest bin and neighbours and identify the vertex of such parabole as the harmonic frequency. If I remember correctly, 5x-10x padding was used.

[hiban]

If acquiring a longer real data sample is not possible, it's better to present the measurement results as is, based solely on the available data, rather than extending the signal with artificial content. Padding with synthetic data can introduce subtle distortions that may lead to significantly inaccurate results, which might not be immediately noticeable. In fact, substantial errors can emerge under specific signal conditions, making the measurement unreliable.

Thanks for the explanations.

Now that you're around, i'm currently measuring SNR, SINAD, THD, THD+N and SFDR. I also considered calculating ENOB, but i'm not sure if it's something worth calculating when the oscilloscope i take the signal data from uses an 8-bit ADC.

[radiolistener]

ENOB is essentially an alternative representation of SINAD, expressed in bits rather than decibels:

ENOB = (SINAD - 1.76) / 6.02

[gf]

Zero-padding cannot overcome the fundamental limitation caused by insufficient information in the input signal. The only effective way to improve frequency resolution is to acquire a longer segment of actual signal data - not artificial extension with zeros.

The (infinite) windowed signal is considered to be zero prior to the window and past the window anyway, so zero padding the points falling into the window just makes some of the implied zeros outside the window explicit. So it is not really an "artificial" extension. The DTFT spectrum of the (infinite) windowed signal is continuous, and we use the DFT only as a tool to sample the continuous DTFT spectrum of the windowed signal at discrete frequency intervals. With zero-padding, we sample the same DTFT spectrum at narrower frequency intervals. In both cases, we calculate exact samples of the continuous DTFT spectrum. So IMO there is nothing wrong with zero padding per se, just as there is nothing wrong per se with sampling a continuous time domain signal at a higher sample rate. Some calculations (like approximating an integral with a sum) can even be more accurate when the samples have a narrower spacing.

If you have 1200 data points and only a power-of-2 FFT implementation available, which choice is less lossy? To discard 176 points and do a 1024 point DFT, or to keep all 1200 points, zero pad to 2048 (4096, 8192, etc.), and do a 2048 (4096, 8192, etc.) point DFT?

You are right, of course, that zero padding still does not increase the duration of the window, and the shape of the window function and the duration of the window determine the RBW and selectivity of the spectrum analysis filter.

ENOB is essentially an alternative representation of SINAD, expressed in bits rather than decibels:
ENOB = (SINAD - 1.76) / 6.02

Let me add, that THD+N is essentially yet another (reciprocal) representation of SINAD, too.

[gf]

Now that you're around, i'm currently measuring SNR, SINAD, THD, THD+N and SFDR. I also considered calculating ENOB, but i'm not sure if it's something worth calculating when the oscilloscope i take the signal data from uses an 8-bit ADC.

But this applies basically to all these measurements. You need an external notch filter (for the fundamental frequency) and LNA if you want to measure the SINAD of the signal and not the SINAD of the scope (unless the signal has a significantly worse SINAD than the scope).