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

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=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: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; end zhi(1,1)=mag1(1); zhi(2,1)=ang1(1); zhi(1,L1)=mag1(L1); zhi(2,L1)=ang1(L1); dotstr{1,1}=num2str(mag1(1),%8.3f); dotstr{2,1}=num2str(ang1(1),%8.3f); dotstr{1,L1}=num2str(mag1(L1),%8.3f); dotstr{2,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); dotstr{1,i}=num2str(mag1(i),%8.3f); dotstr{2,i}=num2str(ang1(i),%8.3f); end %计算Irms和THD I1rms=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; dotstr{2,1}=; else Irms2=Irms2+0.5*mag1(1)^2; end Irms=sqrt(Irms2) THD=sqrt(Irms2/I1rms^2-1)*100 if mdflag for i=1:L1 text(x(i),zhi(1,i),dotstr{1,i}); end end subplot(212); hdl2=stem(f1(1:L1),ang1(1:L1),m:o); set(hdl2,MarkerFaceColor,red,Marker,