使用Matlab的FFT分析六脉动整流的谐波

matlab谐波分析总结

一基本思路

为直观分析显示整流装置的谐波特性，使用matlab的simulink搭建整流电路，利用matlab的fft函数分析其电压与电流波形的谐波特性，并利用matlab的绘图工具，直观的显示谐波的相关参数。输出详细参数到文件。

包括以下想法：

1：用simulink搭建一个由多个不同幅值及相位的正弦波，输出到workspace的simout参数，主要是为了验证算法的正确性。

２：算出THD%

二算法及验证

1：Sine叠加输出

sine.mdl文件

其中含4个Sine Wave，其参数如下表所示。

Sinewave Amplitud

e bias Frequency(rad/sec

)

Phase(rad) Sample

Time

1 2 0.7 2*pi*50 0 -1

2 0.5 0 2*pi*50*5 Pi/180*90 -1

3 1 0 2*pi*50*9 pi/180*45 -1

4 0.3 0 2*pi*50*26 Pi/180/(-135) -1

表达的波形为

f(t)=2*sin(2*pi*50*t) +0.5*sin(2*pi*50*5*t+pi/2)

+1*sin(2*pi*50*9*t+pi/4) +0.3*sin(2*pi*50*26*t-pi*3/4)

为不同幅值与相位的50Hz的基波，5次、9次、26谐波的叠加。含基波、奇次、偶次、高次谐波。在基波上加了0.7的偏置，模拟直流分量。

示波器输出到workspace的参数名

仿真参数10个周波，每周波采样点2048个

使用1/50/2048的采样频率，是为了每个周波采2048个点，便于准确的FFT分析。理论上可以分析1024次以内的谐波。

simulink的scope的输出

simulink的workspace的输出

ScopeData.signals.values共10*2048个点。之所以采10个周波，是为了保证可以避开初始的过渡状态，虽然当前的仿真没有过渡状态，但六脉动整流如果负载有电容的话会有。实际进行FFT分析时，只抽取一个稳定状态周期的2048个点。

2：代码

程序名分fenxi_FFT.m

%基波频率为50Hz

N=2048*4;%每周波采2048个点

fp=30;%拟显示的倍频次数，需小于N/2 用于使前面的参数显示清晰

wq=5;%稳定周期个数，为保证下行程序数组不出错，workspce的数据周期必须大于wq值kh=0.01;%幅值比例系数，输出到文本时太小的谐波参数不输出

fid=fopen('output.txt','wt');%准备写入文件，文件在当前目录下

fprintf(fid,'每周波采样点数=%d\n',N);

y =ScopeData.signals.values ; %读取workspce数据,第一个点前面为示波器输出到workspce参数名

t=0:0.02/N:0.02/N*(size(y)-1);%每个周期为0.02秒，每周期采Ｎ个点，共size(y)-1个点

subplot(2,1,1);

plot(t,y);

xlabel('时间/s');

ylabel('幅值');

title('原始信号');

grid on;

subplot(2,1,2);

pn=wq*N;%取数据的基准点，主要目的是等系统稳定

plot(t(pn:pn+N-1),y(pn:pn+N-1));%数据周期必须小于wq，将出错

xlabel('时间/s');

ylabel('幅值');

title('抽取信号');

grid on;

figure(2);

subplot(2,1,1);

z=fft(y(pn:pn+N-1),N);

mag=abs(z)/N*2; %/N*2还原为真实的幅值

mag(1)=mag(1)/2;% 直流幅值要再除以2

f=0:N/2; %可分析N/2倍频率，N=128点时，可解析64倍频相关参数

bar(f(1:fp),mag(1:fp));%频谱直方图bar 火柴杆stem

%bar(f(2:fp),mag(2:fp));%除去直流分量显示

%bar(f(3:fp),mag(3:fp));%除去基波分量显示

grid on;

xlabel('频率/倍');

ylabel('幅值');

title('FFT分析结果');

fprintf(fid,'直流分量=%8.4g\n',mag(1));

fprintf(fid,'基波峰值=%8.4g\n',mag(2));

%求THDu

Uhf=0;

for i=3:(size(mag)-1)/2%去除直流与基波各次谐波平方加奈奎斯特频率为采样频率的一半

if((mag(i)/mag(2))>kh) %输出各次谐波电流

fprintf(fid,'%d次谐波峰值%8.3g \t %2.3g%% \n',i-1,mag(i),mag(i)/mag(2)*100);

end

Uhf =Uhf+(mag(i))^2;

end

Uhf=sqrt(Uhf);

THDu=Uhf*100/mag(2);%比上基波乘以100%

fprintf(fid,'总谐波畸变率THD=%6.3f\n',THDu);

subplot(2,1,2);

theta=angle(z)/pi*180+90;%不知为何要+90才对

%考虑matlab计算误差，幅值接近0的地方需将相位强制赋0，因为实部误差小/虚部误差小求出的角度不为0，但无意义

for i=1:size(mag)-1

if((mag(i)/mag(2))

theta(i)=0;

end

end

bar(f(1:fp),theta(1:fp));%频谱直方图bar 火柴杆stem

grid on;

xlabel('频率/倍');

ylabel('相位/度');

title('FFT分析结果');

YT=0:45:460;

set(gca,'YTickMode','manual','YTick',YT)

set(gca,'YMinorTick','on');

fclose(fid);%关闭文件

3：说明

（1）程序显示4幅图，原始波形、抽取的一个周波、FFT后的幅值、FFT后的相位。（2）t为每个采样点对应的时间，从０开始，步长为１/50//N秒

（3）参数wq=7。从第7个周期开始去N个点。该参数是为了保证取一个稳定的周期，避开过渡周期。该参数可根据观测波形后设定。为保证下行程序数组不出越线错误，workspce 的数据周期必须大于wn值。

（4）fp=30为拟显示的倍频次数。理论上每周期采N点fp可设为N/2。为清晰显示，故设为30。该参数可更改为≦N/2的值。

（5）程序运行过程中发现，当某次谐波的幅值很小时，发现其相位变化很大，分析其原因认为是FFT结果实部虚部都很小，误差大，故加了一判断，当幅值很小时相位强制置0。考虑阅读习惯，相位以角度显示。