zuye3.m

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

clc;  %用levinson算法求功率谱
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=3;
r=zeros(1,p);  %自相关函数
a=zeros(p,p);  %每阶的系数 定义成一个二维数组 
d=zeros(1,p); %方差值
w=[0:0.01:pi];
tem=0;
%w=[1:50];
for j=1:100;
    tem=tem+x(j)^2;
end
mm=tem/100;
nn=mm;
tem=0;
for j=1:99
    tem=tem+x(j)*x(j+1);
end
r(1)=tem/99;   %自相关函数的第一个值
tem=0;
a(1,1)=-r(1)/mm;
d(1)=(1-abs(a(1))^2)*mm;
tem1=0;
for j=2:p      %第二阶开始
    k=j;          
    for m=1:(100-j)
        tem=tem+x(m)*x(m+k);
    end
    r(j)=tem/(100-k);
    tem=0;
    for j=1:k-1;
        tem1=tem1+a(k-1,j)*r(k-j);
    end
    a(k,k)=-(r(k)+tem1)/d(k-1);
    tem1=0;
    for j=1:k-1;
        a(k,j)=a(k-1,j)+a(k,k)*a(k-1,k-j);
    end
    d(k)=(1-abs(a(k,k))^2)*d(k-1);
end
y=1;
for j=1:p
      y=y+a(p,j)*exp(-i*w*j);
 %    y=y+a(p,j)*exp(-i*j*w);
end
%for m=1:63
 %   y(m)=1+a(p,1)*exp(-i*w(m))+a(p,2)*exp(-i*2*w(m))+a(p,3)*exp(-i*3*w(m));
 %end
Pyy = y.* conj(y);
Pyy=d(p)./Pyy;
%Pyy=10*log10(Pyy);
plot(Pyy)