MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

# Thread Subject: How do I calculate THD in MATLAB?

Subject: How do I calculate THD in MATLAB?

From: chikha said

### chikha said

Date: 2 Feb, 2011 06:58:03

Message: 1 of 9

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.

Subject: How do I calculate THD in MATLAB?

From: Wayne King

### Wayne King

Date: 2 Feb, 2011 11:27:05

Message: 2 of 9

"chikha said" <chikhasaid@gmail.fr> wrote in message <iiav9r$461$1@fred.mathworks.com>...
> 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

Subject: How do I calculate THD in MATLAB?

From: chikha said

### chikha said

Date: 2 Feb, 2011 14:36:03

Message: 3 of 9

"Wayne King" <wmkingty@gmail.com> wrote in message <iibf29$6fq$1@fred.mathworks.com>...
> "chikha said" <chikhasaid@gmail.fr> wrote in message <iiav9r$461$1@fred.mathworks.com>...
> > 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

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

Subject: How do I calculate THD in MATLAB?

From: Wayne King

### Wayne King

Date: 2 Feb, 2011 17:09:04

Message: 4 of 9

"chikha said" <chikhasaid@gmail.fr> wrote in message <iibq4j$8pl$1@fred.mathworks.com>...
> "Wayne King" <wmkingty@gmail.com> wrote in message <iibf29$6fq$1@fred.mathworks.com>...
> > "chikha said" <chikhasaid@gmail.fr> wrote in message <iiav9r$461$1@fred.mathworks.com>...
> > > 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
>
> 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

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

Subject: How do I calculate THD in MATLAB?

From: chikha said

### chikha said

Date: 2 Feb, 2011 19:11:21

Message: 5 of 9

"Wayne King" <wmkingty@gmail.com> wrote in message <iic33g$483$1@fred.mathworks.com>...
> "chikha said" <chikhasaid@gmail.fr> wrote in message <iibq4j$8pl$1@fred.mathworks.com>...
> > "Wayne King" <wmkingty@gmail.com> wrote in message <iibf29$6fq$1@fred.mathworks.com>...
> > > "chikha said" <chikhasaid@gmail.fr> wrote in message <iiav9r$461$1@fred.mathworks.com>...
> > > > 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
> >
> > 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
>
> 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
thank you thank you relay you help me. now I'm try what you said

Subject: How do I calculate THD in MATLAB?

From: Alex

### Alex

Date: 28 Feb, 2011 06:18:04

Message: 6 of 9

"Wayne King" <wmkingty@gmail.com> wrote in message <iic33g$483$1@fred.mathworks.com>...
> "chikha said" <chikhasaid@gmail.fr> wrote in message <iibq4j$8pl$1@fred.mathworks.com>...
> > "Wayne King" <wmkingty@gmail.com> wrote in message <iibf29$6fq$1@fred.mathworks.com>...
> > > "chikha said" <chikhasaid@gmail.fr> wrote in message <iiav9r$461$1@fred.mathworks.com>...
> > > > 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
> >
> > 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
>
> 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

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?

Subject: How do I calculate THD in MATLAB?

From: Wayne King

### Wayne King

Date: 28 Feb, 2011 11:48:28

Message: 7 of 9

"Alex " <specialedster@gmail.com> wrote in message <ikfems$ksk$1@fred.mathworks.com>...
> "Wayne King" <wmkingty@gmail.com> wrote in message <iic33g$483$1@fred.mathworks.com>...
> > "chikha said" <chikhasaid@gmail.fr> wrote in message <iibq4j$8pl$1@fred.mathworks.com>...
> > > "Wayne King" <wmkingty@gmail.com> wrote in message <iibf29$6fq$1@fred.mathworks.com>...
> > > > "chikha said" <chikhasaid@gmail.fr> wrote in message <iiav9r$461$1@fred.mathworks.com>...
> > > > > 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
> > >
> > > 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
> >
> > 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
>
> 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?

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

Subject: How do I calculate THD in MATLAB?

From: Mohsen

### Mohsen

Date: 11 Aug, 2012 11:17:07

Message: 8 of 9

Dear sir,
I was wondering if you'd mind explaining these lines more:
>mspec = msspectrum(hper,x,'Fs',1e3,'NFFT',length(x));
> datavec = mspec.Data(61:60:181);
"1e3" is total time or sth else?
"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>...
> "chikha said" <chikhasaid@gmail.fr> wrote in message <iiav9r$461$1@fred.mathworks.com>...
> > 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

Subject: How do I calculate THD in MATLAB?

From: Stephane

### Stephane

Date: 10 Oct, 2013 13:02:08

Message: 9 of 9

%Following the proposed solution by Wayne, I will explain a method in order to %compute THD using the FFT.

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.