FFT谐波频谱相谱分析Matlab程序二则.doc

1、FFT谐波频谱相谱分析Matlab程序二则by ggihhimm程序一：从Excel表中读取数据，然后分析其频谱、相谱和总谐波畸变率，并画出茎状图（stem），图示如下（测试数据含直流分量）：%该程序从Excel表中读取数据，然后分析其频谱、相谱和总谐波畸变率%一组数据+采样率 显示频谱 下面的阀值应适当调大 注意I1rms的定义clear all;close all;clc;%y=xlsread(HCwithDI_2.xls); %一列数据：excel表格的第一列要知道数据的采样率，为%保证结果的准确性，最好为整数倍周期L=length(y);mdflag=1; %是否标点,1是，0否del

2、=1e-1; %删除的阀值，根据数据幅值大小来设置Fs=50*L; %采样率f0=50; %信号基波频率%f=0:L-1;f=f*Fs/L/f0; yfft=fft(y);yfft(L/2+1:L)=;f(L/2+1:L)=;mag=abs(yfft)/L*2;ang=angle(yfft)*180/pi;mag(1)=mag(1)*0.5; mag1=mag;mag1(mag=del)=; ang1=ang;ang1(mag=del)=;f1=f;f1(mag=del)=; L1=length(f1);subplot(211);hdl1=stem(f1(1:L1),mag1(1:L1),b:

3、o);set(hdl1,MarkerFaceColor,red,Marker,square);ylabel(各次谐波幅值/A);tem=get(gca,XLim);xwidth=tem(end)-tem(1);x=zeros(1,L1);zhi=zeros(2,L1);dotstr=cell(2,L1);xright=1.5;x(1)=f1(1)+xwidth/100*xright; if f1(L1)xwidth-1 x(L1)=f1(L1)+xwidth/100*xright;else x(L1)=f1(L1)-xwidth/100*10;endzhi(1,1)=mag1(1);zhi(2

4、,1)=ang1(1);zhi(1,L1)=mag1(L1);zhi(2,L1)=ang1(L1);dotstr1,1=num2str(mag1(1),%8.3f);dotstr2,1=num2str(ang1(1),%8.3f);dotstr1,L1=num2str(mag1(L1),%8.3f);dotstr2,L1=num2str(ang1(L1),%8.3f);for i=2:L1-1 x(i)=f1(i)+xwidth/100*xright; zhi(1,i)=mag1(i); zhi(2,i)=ang1(i); zhi(1,i)=mag1(i); zhi(2,i)=ang1(i);

5、 dotstr1,i=num2str(mag1(i),%8.3f); dotstr2,i=num2str(ang1(i),%8.3f);end%计算Irms和THDI1rms=max(mag1(abs(f1-1)0.5)/sqrt(2);Irms2=0.5*sum(mag1(2:L1).2);if f1(1)=0 Irms2=Irms2+mag1(1)2; dotstr2,1=;else Irms2=Irms2+0.5*mag1(1)2;endIrms=sqrt(Irms2)THD=sqrt(Irms2/I1rms2-1)*100 if mdflag for i=1:L1 text(x(i),

6、zhi(1,i),dotstr1,i); endendsubplot(212);hdl2=stem(f1(1:L1),ang1(1:L1),m:o);set(hdl2,MarkerFaceColor,red,Marker,o);ylabel(各次谐波相角/circ);if mdflag for i=1:L1 text(x(i),zhi(2,i),dotstr2,i); endend程序二：从Excel表中读取数据，然后分析其频谱、相谱和总谐波畸变率，并画出图（plot），图示如下（测试数据含直流分量）：%该程序从Excel表中读取数据，然后分析其频谱、相谱和总谐波畸变率%一组数据+采样率 阀值

7、可以调小一些 注意I1rms的定义clear all;close all;clc;%y=xlsread(HCwithDI_2.xls); %一列数据L=length(y);del=1e-1; %删除的阀值，根据数据幅值大小来设置Fs=1e4; %采样率f0=50; %信号基波频率%f=0:L-1;f=f*Fs/L/f0; yfft=fft(y);yfft(L/2+1:L)=;f(L/2+1:L)=;mag=abs(yfft)/L*2; ang=angle(yfft)*180/pi;mag(1)=mag(1)*0.5; mag1=mag;mag1(mag=del)=;ang1=ang;ang1(

8、mag=del)=;f1=f;f1(mag=del)=; L1=length(f1);subplot(211);plot(f1(1:L1),mag1(1:L1),b:o);ylabel(各次谐波幅值/A);subplot(212);plot(f1(1:L1),ang1(1:L1),m:s);ylabel(各次谐波相角/circ);%计算Irms和THDI1rms=max(mag1(abs(f1-1)0.5)/sqrt(2);Irms2=0.5*sum(mag1(2:L1).2);if f1(1)=0 Irms2=Irms2+mag1(1)2;else Irms2=Irms2+0.5*mag1(1)2;endIrms=sqrt(Irms2)THD=sqrt(Irms2/I1rms2-1)*100