%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;