zuye4.m

这是用matlab编的burg算法和lenvinson算法

clc;   %用burg算法求功率谱
clear;
x=[101 82 66 35 31 7 20 92 154 125 85 68 38 23 10 24 83 132 131 118 90 67 60 47 41 21 16 6 4 7 14 34 45 43 48 42 28 10 8 2 0 1 5 12 14 35 46 41 30 24 16 7 4 2 8 17 36 50 62 67 71 48 28 8 13 57 122 138 103 86 63 37 24 11 15 40 62 98 124 96 66 64 54 39 21 7 4 23 55 94 96 77 59 44 47 30 16 7 37 74];    
p=15;
w=[0:0.01:pi];
e=zeros(3,100);  %前向误差
b=zeros(3,100);   %后向误差
d=zeros(1,3);      %方差系数
app=zeros(1,3);
a=zeros(p,p);
qq=zeros(1,100);
xe=x;   %e0n
xb=x;   %b0n
%e(1,:)=x;     %e0n
%b(1,:)=x;     %b0n
fangcha=sum(x(:).^2)/100;
qq=[0 xb(1:99)];    %b0n
%k(1)=-2*sum(xe.*x)/(sum(xe(:).^2)+sum(xb(:).^2));
app(1)=-2*sum(qq.*x)/(sum(xe(2:100).^2)+sum(qq(2:100).^2));
a(1,1)=app(1);     %a11的值
d(1)=(1-app(1)^2)*fangcha;
e(1,:)=xe+app(1)*qq;
b(1,:)=qq+app(1)*xe;
qq=zeros(1,100);
for j=2:p
    qq=[0 b(j-1,1:99)];
    app(j)=-2*sum(qq.*e(j-1,:))/(sum(qq(j+1:100).^2)+sum(e(j-1,j+1:100).^2));
    a(j,j)=app(j);
    for k=1:j-1
        a(j,k)=a(j-1,k)+a(j,j)*a(j-1,j-k);
    end
    e(j,:)=e(j-1,:)+app(j)*qq;
    b(j,:)=b(j-1,:)+app(j)*qq;
    d(j)=(1-app(j)*app(j))*d(j-1);
end
y=1;
for j=1:p
      y=y+a(p,j)*exp(-i*w*j);
end
Pyy = y.* conj(y);
Pyy=d(p)./Pyy;
Pyy=10*log10(Pyy);
plot(Pyy)