bg.m

利用burg和Levinson-Durbin算法

%burg算法
function [Pxx]=BG(M,N,x)

d=zeros(1,M+1);
K=zeros(1,M);
e=zeros(M+1,N);
a=zeros(M,M);
%确定初始条件
e(1,:)=x;
b(1,:)=x;
d(1)=sum(x.^2)/N;

%求反射系数K(P)，参数e、b，方差d
for p=1:M
    s1=0;
    s2=0;
    for n=p:N-1
        s1=e(p,n+1)*b(p,n)+s1;
        s2=e(p,n+1).^2+b(p,n).^2+s2;
    end
    K(p)=-2*(s1/s2);
    for n=2:N
        e(p+1,n)=e(p,n)+K(p)*b(p,n-1);
        b(p+1,n)=b(p,n-1)+K(p)*e(p,n);
    end
    d(p+1)=(1-K(p).^2)*d(p);
    a(p,p)=K(p);
    for i=1:p-1
    a(p,i)=a(p-1,i)+K(p)*a(p-1,p-i);
    end
end

%A是计算功率谱的参数a(n)
% A=zeros(1,M+1);
 A(1)=1;%根据Yule-Walker方程第一个值为1
 A(2:M+1)=a(M,:);

%计算功率谱
%因为只有100个值，若fft做256之类的则为100点以外补0，则结果不准确
F=abs(fft(A,100)).^2;
Pxx=d(M+1)./F;