ld.m

利用burg和Levinson-Durbin算法

%Levinson-Durbin算法
function [Pxx]=LD(M,N,x)

a=zeros(M,M);
%求自相关函数R
for m=1:N
    r=0;
     for n=1:(N-m+1)
         r=r+x(n)*x(n+m-1);
     end
     R(m)=r/N;
end

%计算一阶参数a(1,1)和d(1)
a(1,1)=-R(2)/R(1);
d(1)=R(1);
d(2)=(1-a(1,1)^2)*R(1);
%计算M阶参数
for k=2:M
    temp=0;
    for l=1:k-1
         temp=temp+a(k-1,l)*R(k-l+1);
         a(k,k)=-(R(k+1)+temp)/d(k);
    end
    for i=1:k-1
         a(k,i)=a(k-1,i)+a(k,k)*a(k-1,k-i);
    end
    d(k+1)=(1-abs(a(k,k))^2)*d(k);
end

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

%[a,E]=aryule(x,70);

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