Calculate THD using fft

I have a sampled sine wave of a 10kHz signal. It was sampled at 1Msps (1Mhz) so there is 100 samples per period. The total vector is 1024 samples long so its just over 10 periods total.

I have measured the THD of the analog signal using the fft function of a tektronix scope and I know the correct THD to be less than 1%. However, when I try to do this using MATLAB I am not getting the correct results and I expect this is because I am doing it incorrectly. I have pasted my current script below in hopes that somebody would be able to tell me what I am doing wrong. Thanks in advance!

%Use fft to find harmonics
    harmonics = abs(real(fft(data))).^2;
%Assume the fundamental frequency is the highest
    fundamental_bin = min(find(harmonics==max(harmonics)));
    fundamental_power = harmonics(power_bin);
%Find the sum of the rest of the harmonics
    harmonics_power = sum((harmonics((2:numHarmonics).*fundamental_bin)));
%Calculate THD as a ratio of square roots of powers
    THD = 100*sqrt(harmonics_power)/sqrt(fundamental_power);

When I run this script I get a THD of over 8%, which as I explained earlier is far too high. Can anyone point me in the right direction or give me an example of someone who has done this correctly?

[Anton]

Hello,

do you know how can I get 'power_bin'of signal ? I don't know what is it.

Thanks a lot.

[Shubhang Bandyopadhyay]

what is power_bin

Answer 1

[Walter Roberson]

Please go back and reformat your code. If you highlight it and click on the "{} Code" button, it will be formatted for display.

One problem that I can see is that your max(harmonics) is going to be looking at the first output from the fft, which is the DC component totaled over all of the samples. Although a sine wave taken over multiples of full periods would have a theoretical DC component of 0, you need to account for round-off error. And since you do not have an exact multiple of a period, you are going to have a non-zero DC shift, which in turn is going to affect the phase of everything.

At the very least, you should omit the first bin when looking for your max and your harmonics. Better would likely be to also trim to 1000 samples so that the 10 KHz tone and harmonics would lie at exact bins and the DC would be 0 to within round-off.

[Alex]

I see what you mean about the DC offset bin. I will try and omit that as well as trimming the number of samples so that my harmonics end up in exact bins. Thanks for both of those suggestions.

One more problem I am having is knowing how the bins will be spaced. According to another post, the bins will be spaced at Fs/N where N is the number of samples in the vector. Since I have a 1000 sample vector that has been sampled at 1Mhz, does this mean my harmonics will be spaced 1000 apart? And if so, how can I get them if my resulting fft vector is only 1000 values long?

[Walter Roberson]

The bins will be spaced at Fs/N *Hertz*. That's 1000000 / 1000 = 1000 Hz increments, so your 10 KHz primary signal would be at the 11th bin (including the DC bin.)

[Alex]

Great, thank you so much. You have been enormously helpful.

Answer 2

[Wayne King]

The Signal Processing Toolbox makes this easy with thd()