clear all;
clc;
%a1=input('please input the number of a1:');
%a2=input('please input the number of a2:');
a1=-1.5955;%structure AR model using a1 and a2;
a2=0.95;
t=1;%r(0) value
r(1)=-a1*t/(1+a2);
r(2)=(-a2+a1^2/(1+a2))*t;
P=(1-a2)/(1+a2)*((1+a2)^2-a1^2)*t;
for n=3:255;
    r(n)=-a1*r(n-1)-a2*r(n-2);
end
r=[r t r];%classical PSD method using BT.
i=0;
for w=0:pi/1000:2*pi;%sample the 2000 point to the curve of DFT
    i=i+1;
    m=-255:255;
    E=exp(-j*w*m);% generate a vectors,the element subscript range frome -255 to 255.
    PSD(i)=r*E.';% E.' indicate  transpose,and E' indicate transpose and conjugate.
end


i=1:2000;
subplot(4,1,1);
plot(i/1000,abs(PSD(i)),'r');
title('classical PSD method using BT edited by leishengmin');
num=[1 0 0];%PSD using AR model.
den=[1 a1 a2];
w=0:pi/1000:2*pi;
H=freqz(num,den,w);%get the frequence responce of system.
H=H.^2;% H.^2 indicate each element in the vector square.

subplot(4,1,2);

i=1:2000;
%plot(i/1000,(P*abs(H(i)))^2,'r');
plot(i/1000,abs(H(i)),'r');
title('PSD using AR model edited by leishengmin');


clear all;
clc;
N=1000;
a1=-0.195;
a2=0.95;
n=1:N;%gerenate noise signal.
V=sqrt(0.0965)*randn(1,N);
U=zeros(1,N);%initialize vectors.
U(1)=0;
U(2)=0;
for n=3:N;%obtain U(n) by using interative computation.
    U(n)=-a1*U(n-1)-a2*U(n-2)+V(n);
end


lags=255;%estimation correlation function
[C,lag]=xcorr(U,lags,'biased');
disp('C(256)=');
C(256)
disp('test if C(256) is equal to one');
a11=-C(257)*(C(256)-C(258))/(C(256)^2-C(257)^2)
a22=(C(257)^2-C(256)*C(258))/(C(256)^2-C(257)^2)
P1=C(256)+a11*C(255)+a22*C(254)% P1 stand for power of noise.

i=0;%PSD using BT method.
Z=zeros(1,2*N+1);
for w=0:pi/1000:2*pi;
    i=i+1;
    m=-255:255;
    Y=exp(-j*w*m);
    Z(i)=Y*C.';
end
subplot(4,1,3);
i=1:2000;
plot(i/1000,abs(Z(i)),'r');
title('classical PSD method using BT edited by leishengmin');


num=[1 0 0];%PSD using AR model.
den=[1 a1 a2];
w=0:pi/1000:2*pi;
H=freqz(num,den,w);
H=H.^2;
H=H*P1;
subplot(4,1,4);
i=1:2000;
plot(i/1000,abs(H(i)),'r');
title('PSD using AR model edited by leishengmin');