clc;
clear;

fs=1;     
N=10000;
t=0:1/fs:N;

A=20^0.5;
B=A;
C=1;
x1=A*sin(2*0.15*pi*t);
x2=B*sin(2*0.2*pi*t);

noise=C*randn(size(t)); 

x=x1+x2+noise;

M=15;     %階數(shù)
M1=20;     %階數(shù)
%%%%%%%%求自相關(guān)函數(shù)
for m=0:M-1,            % 對r賦初值
    temp=0;
    for n=1:N-M,         
    temp=x(n)*x(n+m)+temp;
    end
    r(m+1)=(1/N)*temp;
end

a(m+1,m+1)=0;
a(1,1)=-r(2)/r(1);
xgm(1)=(1-(a(1,1))^2)*r(1);    %xgm為平方值

p=2;
  while   p<M
p=p+1;
    temp=0;
    for i=1:p-2
    temp=a(p-2,i)*r(p-i)+temp;
    end
    a(p-1,p-1)=-(r(p)+temp)/xgm(p-2);
     for k=1:p-2
        a(p-1,k)=a(p-2,k)+a(p-1,p-1)*a(p-2,p-k-1);
    end
    xgm(p-1)=(1-(a(p-1,p-1))^2)*xgm(p-2);   
  end

%%%%%%%%%%%%%根據(jù)數(shù)學(xué)模型計算功率譜  

q=1; 
for f=0:0.001:0.5-0.001;
w=2*pi*f;
t=0;
for i=1:p-1
t=t+a(p-1,i)*exp((-1)*j*w*i);
end
pxx(q)=(1/abs(1+t))^2;
q=q+1;
end

plot(f/(2*pi),pxx,'r')
f=0.001:0.001:0.5;
subplot(2,1,1);
axis([0 0.8 -20 120]);
plot(f,10*log(pxx));  
title('自編程序')

%%%%%%%%%%%%%%%%%%%%%%%
nfft=100;
subplot(2,1,2);
[pxx1,f1]=pyulear(x,M1,nfft,fs);
axis([0 0.4999 -20 120]);
plot(f1,10*log(pxx1));
title('用MATLB函數(shù)')