recursive_har.m

一种适用于电力系统谐波检测的类小波方法

%递归法谐波检测
% cgz   2008.11.22
% 参数已调整好，不必变化，检测结果很精确了
% fjc为频率检测结果；fuzhijc为幅值检测结果，此方法的缺陷是不能检测谐波的相位。

clear;
%%%%%%%%滤波器频率特性时序列%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a1=0.04;%系数a用来调整滤波器带宽
det=2*pi/(sqrt(3)*a1);%滤波器中参数det
dett1=0.0005;
nx1=20000;%采样点数
tp1=nx1*dett1;%采样时间长度
t1=0:dett1:(nx1-1)*dett1;%时间向量
detf1=1/tp1;%频率间隔
ff1=0:detf1:(nx1-1)*detf1;%频率向量
fc1=50;
omi1=2*pi*fc1;%中心角频率位置
yt1=((det.^4.*(t1.^4))./12-(det.^5.*(t1.^5))./30+(det.^6.*(t1.^6))./90).*exp(-det.*t1).*exp(j.*omi1.*t1);%得时频波形时滤波器函数
%%%%%%%%%递归滤波器的时频波形显示%%%%%%%%%
%plot(t1,real(yt1));
%hold on
%plot(t1,imag(yt1),'r');
fft_c=abs(fft(real(yt1)))/nx1*20;%
%plot(ff1,fft_c);
%hold on
%fft_d=abs(fft(imag(yt1)))/nx1*2000;
%plot(ff1,fft_d,'--');%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%正常计算时序列%%%%%%%%%
fs=2000;
dett=1/fs;%采样周期
nx=5000;%采样点数
tp=nx*dett;%采样时间长度
t=0:dett:(nx-1)*dett;%时间向量
detf=1/tp;%频率间隔
ff=0:detf:(nx-1)*detf;%频率向量

%%%%%%%被测信号%%%%%%%%%
hr0=0;f1=50.1;
hr(1)=25*sqrt(2);deg(1)=10;
hr(2)=0;deg(2)=0;
hr(3)=1.755*sqrt(2);deg(3)=40;
hr(4)=0;deg(4)=0;
hr(5)=0.885*sqrt(2);deg(5)=70;
hr(6)=0;deg(6)=0;
hr(7)=1.125;deg(7)=110;
M=7;f=[1:M]*f1; %设定频率

x=zeros(size(t));
for k=1:M
    x=x+hr(k)*cos(2*pi*k*f1*t+deg(k)*pi/180);
end

%%%%%%%%%%%%%%%%%%%%%%%%
fc=50:50:400;%频率遍历计算中的中心频率向量:Hz
for m=1:8    %%%%%%%%%

omi0=2*pi*fc(m);%中心角频率位置
yt=((det.^4.*(t.^4))./12-(det.^5.*(t.^5))./30+(det.^6.*(t.^6))./90).*exp(-det.*t).*exp(j.*omi0.*t);%滤波器函数

%%%%%%检测部分%%%%%

%%%%%%%%%%递归法检测结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=dett.*conv(x,yt);
%%%%幅值计算%%%%%
p=sqrt(real(w).^2+imag(w).^2);%计算幅值
most=max(fft_c);
fuzhi=p(4601:4680);
fuzhigu=fuzhi/most*2;
fuzhijc(m)=sum(fuzhigu)/80;%幅值取均值的检测结果

%%%%频率计算%%%%%%%%%%%%%
if fuzhijc(m)<0.1
    fjc(m)=0;
else
    theta=atan2(imag(w),real(w));
    k=1:6000;%小于卷积后的系数点数减1即可，由采样点数决定
    dtheta=theta(k+1)-theta(k);
    for i=1:6000   %由于tan函数的主值区间的关系，在相位突变点处进行校正，相差2*pi
        if dtheta(i)<0
            dtheta(i)=2*pi+dtheta(i);
        else
            dtheta(i)=dtheta(i);
        end
    end
    f1=dtheta/dett/2/pi;%f1为检测出的瞬时频率值
    fgu=f1(4601:4680);%7001为稳定的起始点，fgu为检测结果
    fguj=sum(fgu)/80;%频率取均值的检测结果
    fjc(m)=fguj;
end
%%%%%%%幅值校正%%%%%%%%%%%
if fjc(m)>10    %判断频率成分是否存在
    detfjc(m)=abs(fjc(m)-fc(m));%
else
    detfjc(m)=0;
end
h(m)=fc1/detf1+round(detfjc(m)/detf1);%hh(m);频率检测值的位置
if fjc(m)>10
    fuzhigu=fuzhi/fft_c(h(m)+1)*2;
    fuzhijc(m)=sum(fuzhigu)/80;
else
    fuzhijc(m)=0;
end

end   %总体结束
subplot(2,1,1);
stem(fc,fjc);%滤波器的中心频率对应被检测出的频率
axis([0,400,0,400]);
subplot(2,1,2);
stem(fc,fuzhijc);
axis([0,400,0,40]);