burg_hz.m

此代码是高级数字信号处理中用到的谱估计伯格法

clear all
clc
% 产生信号
sample=32;    % 采样点
p=10;         % AR模型的阶数
m=sqrt(-1);
f1=0.05;f2=0.30;f3=0.42;
a(1)=-0.850848;delta=0.101043;
ur(1)=sqrt(delta/2)*randn(1);
ui(1)=sqrt(delta/2)*randn(1);
u(1)=ur(1)+j*ui(1);z(1)=u(1);
x(1)=2*cos(2*pi*f1)+2*cos(2*pi*f2)+2*cos(2*pi*f3)+z(1);
for n=2:sample
    ur(n)=sqrt(delta/2)*randn;
    ui(n)=sqrt(delta/2)*randn;
    u(n)=ur(n)+j*ui(n);z(n)=-a(1)*z(n-1)+u(n);
    x(n)=2*cos(2*pi*f1*n)+2*cos(2*pi*f2*n)+2*cos(2*pi*f3*n)+z(n);
end 
%伯格法
%初值
rxx=0;p0=zeros(1,11);ef=zeros(11,sample);eb=zeros(11,sample);a=zeros(10,10);
for n=1:sample
    rxx=rxx+(abs(x(n)))^2;
end
rxx=(1/sample)*rxx;p0(1)=rxx;
ef(1,:)=x;eb(1,:)=x;ef(1,1)=0;eb(1,32)=0;
%算法
p=10;kk=zeros(1,10);
for k=1:p
    e1=0;e2=0;
    for n=(k+1):sample
        e1=e1+ef(k,n)*(conj(eb(k,n-1)));
        e2=e2+(abs(ef(k,n))^2+abs(eb(k,n-1))^2);   
    end
    kk(k)=(-2)*e1/e2;
    for n=(k+1):(sample-1)
    ef(k+1,n)=ef(k,n)+kk(k)*eb(k,n-1);
    eb(k+1,n)=eb(k,n-1)+conj(kk(k))*ef(k,n);
    end
    for i=1:(k-1)
        a(k,i)=a(k-1,i)+kk(k)*conj(a(k-1,k-i));
    end
       a(k,k)=kk(k);
       p0(k+1)=(1-abs(kk(k))^2)*p0(k);
end
%求PARf
f=-0.5:0.001:0.5;
f0=length(f);
   for z=1:f0
       s3=0;
       for k=1:p
       s3=s3+a(10,k)*exp(-m*2*pi*f(z)*k);
       end
       PARf(z)=delta/(abs(1+s3))^2;
       PARf1(z)=10*log10(PARf(z));
   end
   plot(f,PARf1);
   title('伯格法f2=0.30')