burg_kp.m

用Matlab编写的功率谱估计

function [K,D]=burg_Kp(x,P)%P为阶数
N=length(x);
e=zeros(P,N);
b=zeros(P,N);
e(1,:)=x;
b(1,:)=x;
D=zeros(1,P+1);
for i=1:N
    D(1)=D(1)+x(i).^2;
end
D(1)=D(1)/N;

%K=zeros(P,P);

for p=1:P
    K1=0;
    K2=0;
    for n=p:N-1
        K1=K1+e(p,n+1)*b(p,n);
        K2=K2+e(p,n+1).^2+b(p,n+1).^2;
    end
    K(p,p)=-2*K1/K2;
    for n=2:N
        e(p+1,n)=e(p,n)+K(p,p)*b(p,n-1);%(5-111)
        b(p+1,n)=b(p,n-1)+K(p,p)*e(p,n);%(5-114)
    end
    %D(p+1)=(1-K(p,p)*K(p,p))*D(p)*D(p);
    D(p+1)=sqrt((1-K(p,p).^2)*D(p)*D(p));
    for i=1:p-1
        K(p,i)=K(p-1,i)+K(p,p)*K(p-1,p-i);
    end
       
end
for i=2:p
    D(i)=D(i)*D(i);
end