%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