burg.m

功率普估计

close all
clear
clc
%%%%%%%%%%%%%%%%%%%    输入信号   %%%%%%%%%%%%%%%%%%%
N = 512;fs=100;IP=80;;
f1 = 0.2*fs;f2 = 0.3*fs;
fan1 = 0; fan2 = 0;
t = 0:1./fs:(N-1)/fs;
x = sin(2*pi*f1*t + fan1) + sin(2*pi*f2*t + fan2);
y = awgn(x,10);
plot(t,x);title('无噪声原始信号');
figure;plot(t,y);title('加噪输入信号');
%%%%%%%%%%%%%%%%%%%   周期图法求功率谱 %%%%%%%%%%%%%%%%
y_z = fftshift(abs(fft(y,N)).^2/N);
df=fs./N;
f=(0:df:df*(N-1))-fs./2;
figure;plot(f,10*log10(y_z));title('周期图法(dB)谱估计');

%%%%%%%%%%%%%%%%%%%%%%% p=0     %%%%%%%%%%%%%%%%%%%%%%
ef=zeros(IP,N);
eb=zeros(IP,N);
ef(1,:)=y;      % ef(0,:)
eb(1,:)=y;      % eb(0,:)
row(1)=y*y'/N;  %ρ=row
k(1)=0;         %k(0)
%%%%%%%%%%%%%计算a(p,i) power(p) p=1 2 3 ...IP%%%%%%%%%% 
a=zeros(IP+1,IP+1);
power(1)=row(1);
for p=2:IP+1
    numerator=0;
    denominator=0;
    for n=p:N
    numerator= numerator+ (-2)*ef(p-1,n)*eb(p-1,n-1);
    denominator=denominator+(ef(p-1,n))^2+(eb(p-1,n-1))^2;
    end
    k(p)=numerator./(denominator+0.0001);  
   row(p)=(1-abs(k(p))^2).*row(p-1);
    power(p)=row(p);
    a(p,p)=k(p);
     for n=p:N
        ef(p,n)=ef(p-1,n)+k(p).*eb(p-1,n-1);
        eb(p,n)=eb(p-1,n-1)+k(p).*ef(p-1,n);
     end
end
k=k(1,2:IP+1);
a=a(2:IP+1,2:IP+1);
for p=2:IP
    for i=1:p-1
        a(p,i)=a(p-1,i)+k(p)*a(p-1,p-i);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%% 输出burg波形 %%%%%%%%%%%%%%%%%%% 
NN=200;%0-2*pi 平分N份
df=fs/NN;
f=(0:df:df.*NN)-fs/2;
f1=-1:.01:1;
w=pi*f1;
%w=(-1:.001:1)*pi;
H=0;
 for i=1:p
 H1=a(p,i)*exp(-j*w*i);
 H=H+H1;
 end;
pxx=power(IP+1)./(abs(1+H).^2);
figure;plot(f,10*log10(pxx));title('Burg法(dB)谱估计')
%%%%%%%%%%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%