first1.asv

作业的第二部分 数字信号处理的基本思想rls算法

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