MATLAB Answers

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

5 views (last 30 days)
Tamura Kentai
Tamura Kentai on 4 Jun 2019
Answered: Jmv on 28 Apr 2020
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
//---------------//

Accepted Answer

Greg Dionne
Greg Dionne on 4 Jun 2019
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
  5 Comments

Sign in to comment.

More Answers (1)

Jmv
Jmv on 28 Apr 2020
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.

Products


Release

R2013b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by