How do I calculate THD in MATLAB?
chikha said
Wayne King
> How can I calculate THD (Total Harmonic Distortion) in a MATLAB file? Please note, I know I can do this in Simulink, I am looking for a way to do this from MATLAB directly.
Hi, You can certainly implement the algorithm of the Simulink block pretty easily using just fft() or using a spectrum object and the msspectrum method from the Signal Processing Toolbox.
For example:
% Creating a signal
t = linspace(0,1,1e3);
x = cos(2*pi*60*t)+0.05*sin(2*pi*120*t)+0.01*cos(2*pi*180*t);
% Constructing the spectral analysis object
hper = spectrum.periodogram;
% Getting the mean-square spectrum. These are square magnitudes
mspec = msspectrum(hper,x,'Fs',1e3,'NFFT',length(x));
% Creating a data vector with the fundamental and two harmonics
datavec = mspec.Data(61:60:181);
% Calculating the THD
THD = sqrt(sum(datavec(2:end)))/sqrt(datavec(1))
Multiply by 100 to express as a percentage.
Using fft
xdft = fft(x);
datavec = xdft(61:60:181);
THD = norm(datavec(2:end),2)/norm(datavec(1),2)
The trick is for you to correctly identify the DFT bins that contain your fundamental and harmonics, but that should not be difficult. You can manipulate the NFFT input to fft() or msspectrum() so that your fundamental and harmonics fall on a DFT bin.
Hope that helps,
Wayne
chikha said
thank you a lot
But I don't understand this:
Fs,NFFT,
how you make this vecter
datavect=mspec.Data(61:60:181);
i want change this data beceause the frréquence is 50 for my state. and the time change like that
t=0:h:tf
tf=0.08;
h=1/fs;
fs=1e4;
how i do for this new data
Wayne King
You have to locate the DFT bins that correspond to your frequencies of interest. The DFT bins are spaced at Fs/N where Fs is the sampling frequency and N is the length of the input vector. If you use NFFT, the bins are spaced at Fs/NFFT but you should be aware that doing this does not improve your frequency resolution, it can help you to identify peaks. If you're sampling at 10 kHz, you have to figure out given your N, or NFFT, where your fundamental and harmonics are located so that you can extract the correct Fourier coefficients.
Wayne
chikha said
Alex
Hey Wayne,
Sorry to try and revive a dead thread, but I have a specific question about finding the DFT bins. If I have a sampled signal that is sampled at 1Mhz but only 1000 samples long, does this mean that my DFT bins will be spaced 1000 apart (Fs/N)? And if so, doesn't this mean I'll only be able to get the fundamental frequency since my resulting fft vector will only have 1000 samples?
Wayne King
Yes, your DFT bins will be spaced 1 kHz apart. You can zero pad the data to interpolate the DFT between those bins, but your frequency resolution is determined by your N and your sampling frequency.
As far as:
"And if so, doesn't this mean I'll only be able to get the fundamental frequency since my resulting fft vector will only have 1000 samples?"
That depends on your data and you have not provided those details.
Wayne
Mohsen
I was wondering if you'd mind explaining these lines more:
"Fs" is sampling time or the number of samples per second. because my sampling time is 20e-6 s, but there are samples in my signal every 2e-6 s. (of course after every 10 samples, signal's value changes!)
what dose "(61:60:181)" mean? what for is it. in my case, time is "0.01 to 0.035". and the fundamental frequency is 60 Hz
and what is "NNFT" ?(in your answer to Chikha, I don't know what dose "DFT bins" mean?)
If it is needed, I can send the signal for you.
Best regards
"Wayne King" wrote in message <iibf29$6fq$1...@fred.mathworks.com>...
Stephane
fs = 6400; %sampling frequency
Duree=1; %duration
freq=50; %fundamental frequency of the signal
t = linspace(0,Duree,Duree*fs);
x = 100*sin(2*pi*50*t)+6*sin(2*pi*5*50*t)+5*sin(2*pi*7*50*t);
%So we have three frequencies in the signal : 50Hz, 250 Hz and 350 Hz
%The theoretical value of the THD is:
THD_teorikoa = sqrt(6^2+5^2)/100
% THD_teorikoa = 0.0781
%We will use FFT method in order to compute this THD
% The important element is to select a DFT_bin (DFT=Discrete Fourier Transform) %sufficiently accurate to be a multiple of the signal caracteristic frequencies (50, %250 and 350 Hz)
DFT_bin = fs/length(x); % = (sampling frequency)/(sample length)
% DFT_bin = 1, then doing the FFT, the frquency 50 Hz, 350Hz and 350 Hz will be multiple of this DFT_bin frequency.
% If you want more explanation you can put this words "Amplitude Estimation and %Zero Padding" in the research tool of matlab
%Then now we use fft to compute the DFT and its power:
m = length(x); % Window length (number of samples)
n = pow2(nextpow2(m)) ; % Transform length. This formula is used in order to accelerate the DFT computation, knowing that the computation is easiest if the %length is a multiple of two (prime number)
Y = fft(x,n)/n; %DFT computation
f = fs/2*linspace(0,1,n/2+1); % Frequency range (the maximum of the frequency %range is nyquist frequency (Fs/2)+1
% Plot single-sided amplitude spectrum.
plot(f,abs(Y(1:n/2+1)))
xlabel('Frequency (Hz)')
ylabel('|Y|')
title('{\bf Periodogram}')
% Selection of the bin where the DFT is computed
bin=fs/n; %The bin is computed taking into account that the frequency range is %divided in n bins. Then computation of the width of a bin will give us information, %where the interesting DFT at frequencies (50,250 and 350 Hz) are computed
% Finally the computatio is:
THD=sqrt(abs(Y(ceil((5*50)/bin+1)))^2+abs(Y(ceil((7*50)/bin+1)))^2)/sqrt(abs(Y(ceil(50/bin+1)))^2)
% It's possible that an other formula will be better to compute the THD.