thd() in matlab workspace shows me an unexpected result, is anything wrong ?

Hi, everyone:

Nice to meet you, this is the first time I ask question here, so if there are something improper, I beg your pardon.

I found something strange, unexpecting result of thd() function in Matlab workspace, below is the sample code:

//-------------------------------------//

Fs = 10000;
t=0:1/Fs:1-1/Fs;
y=sin(2*pi*5*t);
[r,harmpow,harmfreq]=thd(y,Fs,40);
plot(t,y,'Color','r','LineWidth',2);
title( ['y=sin(2*pi*5*t)']);
xlabel( 't-axis') , ylabel( 'sin(2*pi*5*t)' );
grid on;
r_percent = 100*10^(r/20);
r
r_percent
harmpow
harmfreq

//----------------------------------//

The unexpecting thing is the result from Matlab, shows the following:

Since the y=sin(...) should be a pure sine wave with the fundamental frequency=5, so the f(1) = 5, f(2) = 10, ...(here f(n), n=1 to 40, represents the 'harmfreq' parameters inside the code), but interestingly the harmfreq = 5, followed the next =9 (not 10 ,unexpectedly), the other issue (problem ?) is the calculated THD value (which is 'r' inside the code, in dBc), is I am not so wrong, the pure sine wave (whole integer periods data), THD should be approx to '0', at least , should not be so large, 12.4978%, from the plot window, I didn't see any distorted waveform.

So that, could anyone know what did I do anything worng ?

Thank you very much.

result:

//----------------//

r =

-18.0633

r_percent =

12.4978

harmpow =

-3.0103

-21.0736

-109.0255

-317.2018

-325.7079

-332.6195

-321.1158

-319.8021

-324.0444

-324.6667

-326.2448

-324.2538

-326.3275

-332.3429

-336.4523

-336.5048

-325.3973

-332.0665

-325.9998

-335.6720

-321.5481

-345.2346

-328.1761

-322.4158

-331.6181

-326.1659

-331.8427

-323.3940

-324.8792

-320.6322

-319.6517

-318.2413

-325.1491

-318.5866

-326.6455

-315.2376

-322.3228

-329.5249

-321.7790

-326.6853

harmfreq =

5.0000

9.0000

14.0000

19.0000

25.0000

31.0000

35.0000

41.0000

46.0000

49.0000

56.0000

60.0000

64.0000

71.0000

76.0000

80.0000

86.0521

89.0000

94.0000

99.0000

106.0000

109.0000

116.0000

120.0000

126.0000

131.0000

136.0000

141.0000

144.0000

151.0000

154.0000

160.5032

164.0000

170.0878

176.0000

180.5322

184.0000

189.0000

194.0000

201.0000

//---------------//

Answer 1

[Greg Dionne]

The default window used when running the THD function has a wide rolloff which masks the second and third harmonics in your signal. You can see this by running the THD function without output arguments and zooming in.

To fix, this, give it a few more samples to work with:

Fs = 10000;
t=0:1/Fs:20-1/Fs;
y=sin(2*pi*5*t);
thd(y,Fs,40);

When I do this, I get around -300 dB or so... which is a reasonable value for double-precision arithmetic.

Alternatively, you can try computing the spectrum with a different window (e.g. rectangular) since you are testing a sinusoid constrained to be periodic in the input record:

Fs = 10000;
t=0:1/Fs:20-1/Fs;
y=sin(2*pi*5*t);

w = rectwin(numel(y));
[Sxx, F] = periodogram(y,w,numel(y),Fs,'power');
  
% compute thd via a power spectrum
rbw = enbw(w,Fs);
[sineTHD, hPower, hFreq] = thd(Sxx,F,rbw,'power')
  
% plot and annotate the spectrum
thd(Sxx,F,rbw,'power')

Hope this helps

-Greg

[Tamura Kentai]

Hi, Greg:

Thank you for the answers.

In fact, the pure sine wave example is only an demonstration to descrbe...

The real needs is I have some waveform data from scope in reality, the file size is approx to 22Mb in *.CSV format, 1e6 x 2 sized array of [time , data_value], and I need to analyze the THD value from it. I have successfully gotton the THD value from 'simulink' by THD block, though the first time what I thought was thd() in workspace , as described previous, not very good news.

So that, though I could run the sine wave example of previous then ran out the reasonable result (low dB value), but when I go back to the *.CSV , it seems the same, the result of thd from workspace is still higher than from others (simulink or some other tools out of Matlab), I don't know why ? Thoug I would give the *.CSV file , but I don't know where to upload it, and let you try it (if you will...), the *.CSV from scope is stored in only 5~8 cycles (of main fundamental) for the reasons to fitting the screen of scope and gain higher resolution (at least, I thought), though from you description, it seems more periods would be helpful (about rolloff, etc), I don't know the detail so I don't know why.

I appreciate you help. I have solved the problem of getting thd value by 'simulink', so actually I don't have any problem now, but if about how to gain THD value from numerical data by the thd() in workspace , if this is interesting, and there are places to upload the *.CSV , maybe could let me know.

Best regards.

[Greg Dionne]

I'm glad you have a version of THD that is working for you.

If you follow the second example in the help, you can experiment with the value of "beta" in the window (corresponding to the "38" in the help text). Setting beta to zero is equivalent to a rectangular window. Using larger values will gradually broaden the main peak (and recover more of the noise floor).

    % Example 2:
    %   Generate the periodogram of a 2.5 kHz distorted sinusoid sampled
    %   at 48 kHz and compute the thd
 
    load('sineex.mat','x','Fs');
    beta = 38; % 0 is rectangular.
    w = kaiser(numel(x),beta);
    [Sxx, F] = periodogram(x,w,numel(x),Fs,'power');
  
    % compute thd via a power spectrum
    rbw = enbw(w,Fs);
    [sineTHD, hPower, hFreq] = thd(Sxx,F,rbw,'power')
  
    % plot and annotate the spectrum
    thd(Sxx,F,rbw,'power')

I'm not certain what the limitation is on uploading files, but saving as a .mat file may be less. If you just save the data_value as a .mat file, it will take ~9MB instead of 22 MB.

Hope this helps.

-Greg

[Tamura Kentai]

Hi, Greg:

It's nice to hear from you again.

I have just tried your suggestion, and I found the max uploading file size is 5MB, so I strip some tail of the array to make it smaller. The file 'M_array_frm_workspace.mat' I have uploaded, I have tried as possible as I can to make the data in whole integer period of the fundamental frequency which is approx. to 110.1Hz, it seems only 2 cycles, since greater would make it's size larger than 5MB.

The code is as below, but unfortunately I haven't got good result of THD value, as described before, it's still unexpected. My calculation result is -17.4377 in dB, about 13.43%.

For your information.

Thank you, I appreciate your help.

Best regards.

% The loaded array of M in (time, data_value) format
% M(:,1) as 'time' array, M(:,2) as 'data_value' array
load('M_array_frm_workspace.mat','M');
plot(M( :,1),M( : ,2),'Color','b','LineWidth',2);
title( ['Plot waveform from mat-file']);
xlabel( 'time-axis') , ylabel( 'Iu Current of Motor-axis' );   
grid on;

Fs = 1 / ( M(2,1) - M(1,1) );
beta = 38; % 0 is rectangular.
w = kaiser(numel( M(:,2) ),beta);
[Sxx, F] = periodogram( M(:,2) ,w,numel( M(:,2) ),Fs,'power');
  
% compute thd via a power spectrum
rbw = enbw(w,Fs);
[sineTHD, hPower, hFreq] = thd(Sxx,F,rbw,'power')
  
% plot and annotate the spectrum
thd(Sxx,F,rbw,'power')

[Tamura Kentai]

Hi, Greg:

I'm sorry, after another try of changing the 'beta = 38' to 'beta=0', the result is quite reasonable to approx. 1.9448%, very similar to other analysis result from else tools.

I guess I need to know kaiser() function first.

Thank you for the grateful help.

[Greg Dionne]

Sounds good! Glad I could help.

Best,

-Greg

Answer 2

[Jmv]

Hi Greg,

Thanks for the answer provided here, it has helped me as I had a similar problem.

Fs = 10000;
t=0:1/Fs:20-1/Fs;
y=sin(2*pi*5*t);
w = rectwin(numel(y));
[Sxx, F] = periodogram(y,w,numel(y),Fs,'power');
  
% compute thd via a power spectrum
rbw = enbw(w,Fs);
[sineTHD, hPower, hFreq] = thd(Sxx,F,rbw,'power')
  
% plot and annotate the spectrum
thd(Sxx,F,rbw,'power')

If I may ask, the code provided is performing calculation in DB. How can I instead plot in terms of power measurement and percentage. I am follwoing this https://au.mathworks.com/help/signal/ref/db.html but can't seem to get it working.

Thanks for any more assistance you can provide.