thd() in matlab workspace shows me an unexpected result, is anything wrong ?
14 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Tamura Kentai
le 4 Juin 2019
Réponse apportée : Jmv
le 28 Avr 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
//---------------//
0 commentaires
Réponse acceptée
Greg Dionne
le 4 Juin 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 commentaires
Plus de réponses (1)
Jmv
le 28 Avr 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.
0 commentaires
Voir également
Catégories
En savoir plus sur Spectral Measurements dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!