本文已参与「新人创造礼」活动，一同开启创造之路。

运用Matlab核算逆变器输出电压的THD

0.要求

本程序的作用是核算逆变器输出电压的THD（Total Harmonic Distortion，总谐波失真或谐波总畸变率），核算办法是除基波外各次谐波幅值的平方和再开根号，然后与基波（50Hz）幅值的比值即为THD。

要求如下：

  • 体系核算的最高次谐波是18次谐波；
  • 采样数据已供给，在data.mat中，采样频率5000Hz，共5万个点，可用load命令读入；
  • 规划一个低通滤波器（m言语编程或matlab自带工具箱），滤除采样数据中高于最高次谐波的频率，滤波后的数据再用来核算THD；
  • 请用MATLAB言语完结一个函数来核算电压的THD，输出参数为THD核算成果，函数的输入参数为电压采样值x；
  • 核算过程中注意结合课程中所讲过的信号处理办法，如去除直流重量、移动均匀滤波、分段求频谱再均匀等办法；
  • 要求有信号时域波形图（取200点）、滤波前信号的频谱图、滤波后信号的频谱图等；

1 运用load读取信号

loaddata.mat
Fs=5000;           
T=1/Fs;       
L=length(data);

2 运用fft函数制作信号的单边频谱

相关代码如下：

% 制作输入信号频谱
Yin=fft(data);
P2In=abs(Yin/L);
P1In=P2In(1:L/2+1);
P1In(2:end-1) =2*P1In(2:end-1);
fIn=Fs*(0:(L/2))/L;
plot(fIn,P1In)
title('滤波前信号的频谱图')
xlabel('f (Hz)')
ylabel('单侧幅值频谱|P1(f)|')

得到成果如下：

使用Matlab计算逆变器输出电压的THD

3 选择合适的滤波器

FIR滤波器是有限长单位冲激响应滤波器，又称为非递归型滤波器，是数字信号处理体系中最根本的元件，它可以在保证任意幅频特性的一起具有严厉的线性相频特性。一起其单位抽样响应是有限长的，因此滤波器是安稳的体系。

因为输入的信号是逆变器输出电压，需要满意线性相位，故这里运用FIR滤波器。

规划并运用滤波器的代码如下。

% 2
% 运用低通滤波器处理信号
% 最高次谐波是18次谐波，通带0-900Hz，阻带950Hz
Fpass=900;
Fstop=950;
Ap=1;
Ast=30;
lowpassfirFilter=designfilt('lowpassfir','PassbandFrequency',Fpass,...
'StopbandFrequency',Fstop,'PassbandRipple',Ap,...
'StopbandAttenuation',Ast,'SampleRate',Fs);
% fvtool(lowpassfirFilter,"Fs",Fs);
dataOut=filtfilt(lowpassfirFilter,data);

生成滤波器的幅频响应和相频响应图如下：

使用Matlab计算逆变器输出电压的THD

4. 制作输出信号的频谱

做图方式同过程2，得到滤波后信号频谱如下，不难发现19及更高次谐波已被过滤。

使用Matlab计算逆变器输出电压的THD

对应以下程序，为了便利比照，将滤波前函数的频谱和滤波后的频谱进行比照。

% 3
% 滤波后信号的频谱图
Yout=fft(dataOut);
P2Out=abs(Yout/L);
P1Out=P2Out(1:L/2+1);
P1Out(2:end-1) =2*P1Out(2:end-1);
fOut=Fs*(0:(L/2))/L;
figure(1)
subplot(1,2,1)
plot(fIn,P1In)
title('滤波前信号的频谱图')
xlabel('f (Hz)')
ylabel('单侧幅值频谱|P1(f)|')
subplot(1,2,2)
plot(fOut,P1Out)
title('滤波后信号的频谱图')
xlabel('f (Hz)')
ylabel('单侧幅值频谱|P1(f)|'

使用Matlab计算逆变器输出电压的THD

5. 用 MATLAB 言语编写一个函数来核算电压的 THD的函数

为了保证自己编写的thd函数的正确性，首要运用Matlab的函数thd求解。

使用Matlab计算逆变器输出电压的THD

接着介绍如何自己完成thd函数。

首要对输入的信号进行移动均匀值滤波。

定义滤波阶数为meanFitLeval。第n个点的数据为第n-meanFitLeval个点至第n点的均匀值。

% 移动均匀值滤波
 meanFitLeval=23;
 myDataMean=dataOut*0;
 myDataMean(1:meanFitLeval) =dataOut(1:meanFitLeval);
 i=meanFitLeval+1;
 whilei<length(dataOut)
   myDataMean(i) =sum(dataOut(i-meanFitLeval:i))/meanFitLeval;
   i=i+1;
 end

接着去除信号中的直流重量，将信号减去均值。

% 去除直流重量
 myDataMeanAC=myDataMean-mean(myDataMean);

做出去除直流重量并进行均匀值滤波后的信号与原始信号的比照图。

figure(3)
holdon;
plot((1:200)/Fs,vthSig(1:200));
plot((1:200)/Fs,data(1:200));
xlabel('时刻')
ylabel('输出电压')
legend("移动均匀值滤波+去除直流后的信号","原始信号")
holdoff;

使用Matlab计算逆变器输出电压的THD

运用分段求频谱再均匀的办法优化thd的核算。将输入信号按时刻分为10组，求出各组信号的fft成果，并取均匀值，最后依据公式核算thd值。

% 分段求频谱再均匀,输入信号分为10组
 lengthOfPart=length(myDataMeanAC)/10;
 fftPart=zeros(1,lengthOfPart);
 fori=1:10
   thisPart=myDataMeanAC(1+lengthOfPart*(i-1):lengthOfPart*i);
   fftPart=fftPart+abs(fft(thisPart));
 end
 fftPart=fftPart/10;
 P2fftPart=abs(fftPart/lengthOfPart);
 P1fftPart=P2fftPart(1:lengthOfPart/2+1);
 P1fftPart(2:end-1) =2*P1fftPart(2:end-1);
 ansfftPart=Fs*(0:(lengthOfPart/2))/lengthOfPart;
 % 求thd
 sumf218=0;
 fori=3:18
   sumf218= sumf218+P1fftPart(i*50+1)^2;
 end
 VthdMy=10*log10(sqrt(sumf218)/P1fftPart(50+1));

核算得关于输入信号的前18次谐波的THD值，Matlab thd函数成果为-12.059dB，自行完成的 thd函数成果为-12.067dB，自行完成的函数差错在合理范围内，满意要求。

"Matlab Thd=-12.059dB My Thd=-12.067dB"
"Matlab Thd=0.062 My Thd=0.062"