limburge1.m

功率谱估计程序

%burge算法，选用电力系统谐波信号作为输入信号。参考书目：现代信号处理导论，陆传赉，北京邮电大学出版社。96－97页。
clc;clear all;close all;
f1=50;ts=1/1000;n=1:256;N=256;fs=1000;
x(1:128)=2*sin(2*pi*f1*(1:128)*ts);
x(129:256)=sin(2*pi*f1*(129:256)*ts)+0.1*sin(2*pi*350*(129:256)*ts)+0.2*sin(2*pi*400*(129:256)*ts);D=0;M=16;
for n=0:N-1
    ef(1,n+1)=x(n+1);
    eb(1,n+1)=x(n+1);
    D(1)=D+1/N*x(n+1)^2;
end
A=0;B=0;
for n=1:N-1
    A=A+2*ef(1,n+1)*eb(1,n);
    B=B+ef(1,n+1)^2+eb(1,n)^2;
    Z(1)=-A/B;
end
a(1,1)=Z(1);
D(2)=(1-Z(1))^2*D(1);
for n=1:N-1
ef(2,n+1)=ef(1,n+1)+Z(1)*eb(1,n);
eb(2,n+1)=eb(1,n)+Z(1)*ef(1,n+1);
end

for p=2:M
    A=0;B=0;
    for n=p:N-1
    A=A+2*ef(p,n+1)*eb(p,n);
    B=B+ef(p,n+1)^2+eb(p,n)^2;
    end
    Z(p)=-A/B;
    a(p,p)=Z(p);
    for k=1:p-1
        a(p,k)=a(p-1,k)+Z(p)*a(p-1,p-k);
    end
    D(p+1)=(1-Z(p))^2*D(p);
    for n=p:N-1
        ef(p+1,n+1)=ef(p,n+1)+Z(p)*eb(p,n);
        eb(p+1,n+1)=eb(p,n)+Z(p)*ef(p,n+1);
    end
end

for f=0:fs/2
    A1=0;
for k=1:M
A1=A1+a(M,k)*exp(-i*2*pi*f*k/fs);
end
S(f+1)=D(M)/(abs(1+A1))^2;
end
f=0:fs/2;
semilogy(f,S)