burg.m

计算太阳黑子的活动周期

%burg算法求解AR模型，并求功率谱
function  [a,e2]=burg(x,p)
N=length(x);
e(1,:)=x;
b(1,:)=x;
e2(1)=sum(x.^2)/N;
a=zeros(p,p);
kp=[];
for n=2:p+1;
    kp(n-1)=2*(e(n-1,[n:N])*b(n-1,[n-1:N-1])')/(e(n-1,[n:N])*e(n-1,[n:N])'+b(n-1,[n-1:N-1])*b(n-1,[n-1:N-1])');
    a(n-1,n-1)=-kp(n-1); 
    if n>2
        for i=1:n-2
            a(n-1,i)=a(n-2,i)-kp(n-1)*a(n-2,n-i-1); 
        end
    end
    for n2=2:N
        e(n,n2)=e(n-1,n2)-kp(n-1)*b(n-1,n2-1);
        b(n,n2)=-kp(n-1)*e(n-1,n2)+b(n-1,n2-1);
    end    
    
    e2(n)=(1-kp(n-1).^2)*e2(n-1);
end