clear
clc
%LD
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=3;
N=length(x);
a=zeros(p,p);
dert=zeros(1,N);
%求自相關(guān)函數(shù)
fai=zeros(1,p+1);
for (i=0:p)
    x1=[zeros(1,i),x];
    x2=[x,zeros(1,i)];
    fai(i+1)=sum(x1.*x2)/N;
end
%fai
%求AR模型參數(shù)
a(1,1)=-fai(2)/fai(1);
dert(1)=(1-(abs(a(1)))^2)*fai(1);
for (k=2:p)
     a(k,k)=fai(k+1);
    for(l=1:k-1)
        a(k,k)=a(k,k)+a(k-1,l)*fai(k-l+1);
    end 
    a(k,k)=-a(k,k)/dert(k-1);
    dert(k)=(1-(abs(a(k,k)))^2)*dert(k-1);
    for (i=1:k-1)
        a(k,i)=a(k-1,i)+a(k,k)*a(k-1,k-i);
    end
end
%求功率譜
w=0:0.01:2*pi;
M=length(w);
P=zeros(1,M);
for (i=1:M)
    P(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)/P(i);
    P(i)=10*log(P(i));
end
%f=(1:length(w));
%繪制功率譜圖像
%figure(2)
%plot(w,P);title('Levinson-Durbin算法功率譜');grid on;xlabel('w');ylabel('P(w)');axis([0,pi,0,130]);
plot((1:length(w)),P);title('Levinson-Durbin 算法功率譜');grid on;xlabel('');ylabel('P(w)');axis([0,length(w)/2,0,130]);