harpacket_adap.m

自适应谐波小波包应用于电力系统谐波检测的程序

%%%%%%%%谐波小波变换的自适应二带小波包实现%%%%%%%%%%
clear
fs=2048;
fn=fs/2;%奈奎斯特频率
dett=1/fs;%采样周期
nx=1024;%采样点数
tp=nx*dett;%采样时间长度
detf=1/tp;%频率分辨率
t=0:dett:(nx-1)*dett;%时间向量
f=0:detf:(nx-1)*detf;%频率向量
%%%%信号部分%%%%
x=sin(100*pi*t);%+sin(1600*pi*t)
x(90:170)=0.98*x(90:170);
subplot(211);plot(t,x);
fft_x=fft(x);
%subplot(412);stem(f,abs(fft_x));
%%%%%构造二带分割树%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=fn/detf;%奈氏频率对应的点数
fft_fn=fft_x(1:k);%提取的奈氏频率下信号的频率序列,对应根节点%+1
J=7;%分割树的层数
for j=1:J
    L(j,:)=2^(j-1);%每层的节点数
    for jj=1:L(j,:)
        k1=k/L(j,:);%不同层下的子带对应的节点数
        B{j,jj}=fft_fn((jj-1)*k1+1:jj*k1);%不同层下各节点对应的频率序列
        b{j,jj}=k1*jj;%不同层每个频带截至点对应的频域点数，乘detf即为对应频率
    end
end
%%%%计算每个节点对应的谐波小波系数及对应的香农熵%%%
for m=1:J
    for mm=1:L(m,:)
        Bcoff{m,mm}=ifft(B{m,mm})/(sqrt(2^(m-1)));%谐波小波系数，sqrt(2^(m-1))用作对小波变换系数进行能量校正，因能量为平方项，故需要开方
        coff2=abs(Bcoff{m,mm}).^2;%每个节点对应的系数能量序列
        coff_sum=sum(x.^2);%原始信号能量
        Pi=coff2/coff_sum;
        HB{m,mm}=-sum(Pi.*log(Pi));%各频带下的香农熵
    end
end
y=sum(abs(Bcoff{1,1}).^2);%第一分解层能量
yy=sum(abs(Bcoff{2,1}).^2+abs(Bcoff{2,2}).^2);%第二分解层能量
yyy=sum(abs(Bcoff{3,1}).^2+abs(Bcoff{3,2}).^2+abs(Bcoff{3,3}).^2+abs(Bcoff{3,4}).^2);%第三分解层能量
yyyy=sum(x.^2);%原始信号能量
%%%%%%%%%%%%%搜索最优划分%%%%%%%%%%%%%%%
%%%最高层赋初值%%%
for i=1:L(J,:)
    a{J,i}=b{J,i};%a,b为划分的截至频率对应的频域点数，b原始划分，a最优划分
    HA{J,i}=HB{J,i};%最大层对应的熵值及最优熵值
end
%%%%%%开始搜索%%%%%
for n=J-1:-1:1%从次最大层直到第一层
    for nn=1:L(n,:)  %层中点数序列,次最大层开始
        T=HA{n+1,2*nn-1}+HA{n+1,2*nn};
        if HB{n,nn}>T
            a{n,nn}=[a{n+1,2*nn-1},a{n+1,2*nn}];
            HA{n,nn}=T;
        else
            a{n,nn}=b{n,nn};
            HA{n,nn}=HB{n,nn};
        end
    end
end
%%%%%%%%%%计算最优划分a{1,1}下的谐波小波变换系数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
long=length(a{1,1});
for d=1:long
    if d<2
        p=1;
        q=a{1,1}(1);
        c1{d}=zeros(1,p);
        c2{d}=ones(1,q-p);
        c3{d}=zeros(1,nx-q);
        c{d}=[c1{d},c2{d},c3{d}];
    else
        p=a{1,1}(d-1)+1;%
        q=a{1,1}(d);
        c1{d}=zeros(1,p);
        c2{d}=ones(1,q-p);
        c3{d}=zeros(1,nx-q);
        c{d}=[c1{d},c2{d},c3{d}];
    end
    har_mol(d,:)=1/((q-p)*detf*2*pi);%谐波小波频域幅值
    har_fre(d,:)=har_mol(d,:)*c{d};%谐波小波频域特性
    afreq(d,:)=fft_x.*har_fre(d,:);%频域小波系数
    acoff(d,:)=ifft(afreq(d,:));
    signal(d,:)=real(acoff(d,:))*4*pi*(q-p)*detf;
end

subplot(212);plot(t,signal(9,:));%plot(signal');