N=1000;
r0=0;
v=sqrt(0.0965)*randn(N+998,1);
a(1)=0;
a(2)=0;
for k=3:1:N+1000
a(k)= 0.195*a(k-1)-0.95*a(k-2)+v(k-2);  
end; 
for n=1:1:N
u(n)=a(n+1000);   %取u(n)為a(n) 的后1000個值作為采樣值
end;
clear a;
r=zeros(N,1);
r0=0;
r1=zeros(N-1,1);   %將m不為0的N-1個值存入r1
for i=1:1:N
r0=r0+[u(i)^2]/N; 
end;           %計算r(0)
for m=1:N-1
   for n=m+1:N
   r1(m)=r1(m)+[u(n)*conj(u(n-m))]/N;
   end;
end;         %計算m不為0時的自相關函數
r=[r0;r1];   %將所有相關函數存入r，注意數組中的r(m)為實際中的r(m-1)
plot(r)
xlabel('r(m)');
title(['r(0)=',num2str(r(1))]);  %繪出相關函數
a1=-r(2)*[r(1)-r(3)]/[r(1)^2-r(2)^2]   %a1的估計值
a2=-[r(1)*r(3)-r(2)^2]/[r(1)^2-r(2)^2]  %a2的估計值
v2=(1-a2)*[(1+a2)^2-a1^2]/(1+a2)   %噪聲方差的估計值
figure;
w=0:2*pi/999:2*pi;
Pbt=1;
Par=0;   %為繪制PSD賦初值
for n=2:1:256
Pbt=Pbt+r(n)*[exp(-j*w*(n-1))+exp(j*w*(n-1))];
end;
subplot(2,1,1);
plot(w,Pbt)
xlabel('w');
title('Pbt by 張旭紅');    %在第一子圖繪制BT法估計的PSD
H=1./[1+a1*exp(-j*w)+a2*exp(-j*2*w)];
Par=v2.*abs(H).^2;
subplot(2,1,2);
plot(w,Par)
xlabel('w');
title('Par by 張旭紅');     %在第二子圖繪制AR法估計的PSD