Thread Subject:

How do I calculate THD in MATLAB?

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.

"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

"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

"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

"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

"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?

"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

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