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

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)

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','o');

ylabel('各次谐波相角/\circ');

if mdflag

for i=1:L1

text(x(i),zhi(2,i),dotstr{2,i});

end

end

程序二：从Excel表中读取数据，然后分析其频谱、相谱和总谐波畸变率，并画出图（plot），图示如下（测试数据含直流分量）：

%该程序从Excel表中读取数据，然后分析其频谱、相谱和总谐波畸变率

%一组数据+采样率阀值可以调小一些注意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(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和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;

else

Irms2=Irms2+0.5*mag1(1)^2;

end

Irms=sqrt(Irms2)

THD=sqrt(Irms2/I1rms^2-1)*100 %百分值