clear
clc

x=[101,82,66,35,31,7,20,92,154,125,85,68,38,23,10,24,83,132,131,118,90,67,60,47,41,21,16,6,4,7,14,34,45,43,48,42,28,10,8,2,0,1,5,12,14,35,46,41,30,24,16,7,4,2,8,17,36,50,62,67,71,48,28,8,13,57,122,138,103,86,63,37,24,11,15,40,62,98,124,96,66,64,54,39,21,7,4,23,55,94,96,77,59,44,47,30,16,7,37,74];
p=20;
N=length(x);
e=zeros(p+1,N);%前預測誤差
e=x;
b=zeros(p+1,N);%后向預測誤差
b=x;
dert=zeros(1,p+1);%方差
dert(1)=mean(x.*x);
K=zeros(1,p);%部分相關系數
a=zeros(p,p);%AR模型參數
for (k=1:p)
    K(k)=-(2*sum([e(k,k+1:N)].*[b(k,k+2:N),0])/(sum([e(k,k+1:N)].^2+[b(k,k+2:N),0].^2)));
    a(k,k)=K(k);
    if(k>1)
       for (i=1:k-1)
           a(k,i)=a(k-1,i)+a(k,k)*a(k-1,k-i);
       end
    end
    e(k+1,:)=e(k,:)+K(k)*[b(k,2:N),0];
    b(k+1,:)=[b(k,2:N),0]+K(k)*e(k,:);
    dert(k+1)=(1-(K(k)^2))*dert(k);
end
%計算功率譜
w=0:0.01:2*pi;%頻率向量
P=zeros(1,length(w));%功率譜向量
for (i=1:length(w))
    P(1,i)=1;
    for (k=1:p)
        P(i)=P(i)+a(p,k)*exp(-j*w(i)*k);
    end
    P(i)=P(i)*conj(P(i));
    P(i)=dert(p+1)/P(i);
    P(i)=10*log(P(i));
end
%繪制功率譜圖像
plot((1:length(w)),P);title('Burg算法功率譜');grid on;xlabel('');ylabel('P(w)');axis([0,length(w)/2,0,130]);