levinson.m

计算太阳黑子的活动周期

%Levinson算法求AR模型参数，并求得功率谱
function [a,e2]=Levinson(x,p)
R=r(x);
a=zeros(p,p);                        %用来存放AR模型参数
e2=zeros(1,p);                       %用来存放误差功率
a(1,1)=-R(2)/R(1);                   %计算一阶模型参数a11；
e2(1)=(1-(abs(a(1,1)))^2)*R(1);      %计算一阶模型参数误差功率；
for k=2:p
   s=0;
    for l=1:k-1
        s=s+a(k-1,l)*R(k-l+1)';
    end
    a(k,k)=-(R(k+1)+s)/e2(k-1);
    for i=1:k-1
        a(k,i)=a(k-1,i)+a(k,k)*a(k-1,k-i);
    end
    e2(k)=(1-abs(a(k,k))^2)*e2(k-1);  
end