附录4 递推最小二乘算法在白噪声和有色噪声下工作的源程序

1、nsr=0.1;v=randn(1,1024);v1=nsr*v;b2=1 0.5;a2=1 -0.9 0.95;v2=filter(b2,a2,v1);y1=x+v1;y2=x+v2;附录4递推最小二乘算法在白噪声和有色噪声下工作的源程序%定义噪信比的取值范围%产生指定方差的白噪声信号%在此基础上生成有色噪声%分别将白噪声和有色噪声叠加在理想输出上c10=0.001 0.001 0.001 0.001;c20=0.001 0.001 0.001 0.001;p10=10A6*eye(4);p20=10A6*eye(4);E=0.000000005;c1=c10,zeros(4,1023);c

2、2=c20,zeros(4,1023);e1=zeros(4,1024);e2=zeros(4,1024);for k=3:1024h1=-y1(k-1) -y1(k-2) u(k-1) u(k-2);h2=-y2(k-1) -y2(k-2) u(k-1) u(k-2);z1=h1*p10*h1+1; z11=inv(z1);z2=h2*p20*h2+1; z21=inv(z2);k1=p10*h1*z11;k2=p20*h2*z21;d1=y1(k)-h1*c10; c11=c10+k1*d1;d2=y2(k)-h2*c20; c21=c20+k2*d2;p11=p10-k1*k1*h1*p

3、10*h1+1;p21=p20-k2*k2*h2*p20*h2+1;e11=(c11-c10)./c10;e21=(c21-c20)./c20;e1(:,k)=e11;e2(:,k)=e21;c1(:,k)=c11;c2(:,k)=c21;c10=c11;c20=c21;p10=p11;p20=p21;if (e11E)&(e21E) j=k,break;end%准备求k的值%求出k的值%求出theta的估计值%求出p的值%求出两次估计值之间的相对误差%保存相对误差%保存估计值%将本次估计值作为下次估计的初值enda11=c1(1,:);a12=c1(2,:);b11=c1(3,:);b12=

4、c1(4,:);a21=c2(1,:);a22=c2(2,:);b21=c2(3,:);b22=c2(4,:);ea11=e1(1,:);ea12=e1(2,:);eb11=e1(3,:);eb12=e1(4,:);ea21=e2(1,:);ea22=e2(2,:);eb21=e2(3,:);eb22=e2(4,:);figure% 绘图i=1:j; plot(i,a11(i),r,i,a12(i),b,i,b11(i),m,i,b12(i),k);title(NSR=0.1时白噪声影响下各参数辨识结果随递推次数的变化)； xlabel(递推次数);ylabel(各参数辨识结果)；axis t

5、ight;grid on;figureplot(i,a21(i),r,i,a22(i),b,i,b21(i),m,i,b22(i),k);title(NSR=0.1时有色噪声影响下各参数辨识结果随递推次数的变化);xlabel(递推次数);ylabel(各参数辨识结果)；axis tight;grid on;figureplot(i,ea11(i),r,i,ea12(i),b,i,eb11(i),m,i,eb12(i),k);title(NSR=0.1时白噪声影响下各参数辨识误差随递推次数的变化)； xlabel(递推次数);ylabel (各参数辨识误差)；axis tight;grid o

6、n;figureplot(i,ea21(i),r,i,ea22(i),b,i,eb21(i),m,i,eb22(i),k);title(NSR=0.1时有色噪声影响下各参数辨识误差随递推次数的变化);xlabel(递推次数);ylabel (各参数辨识误差)；axis tight;grid on;nsr1=0.4;v1=nsr1*v;v2=nsr1/nsr*v2;y1=x+v1;y2=x+v2;%重新定义噪信比的取值范围%产生指定方差的白噪声信号%将原有色噪声扩大4倍%分别将白噪声和有色噪声叠加在理想输出上c10=0.001 0.001 0.001 0.001;c20=0.001 0.001

7、0.001 0.001;p10=10A6*eye(4);p20=10A6*eye(4);E=0.000000005;c1=c10,zeros(4,1023);c2=c20,zeros(4,1023);e1=zeros(4,1024);e2=zeros(4,1024);for k=3:1024h1=-y1(k-1) -y1(k-2) u(k-1) u(k-2);h2=-y2(k-1) -y2(k-2) u(k-1) u(k-2);z1=h1*p10*h1+1; z11=inv(z1);z2=h2*p20*h2+1; z21=inv(z2);% 准备求 k 的值k1=p10*h1*z11;k2=p

8、20*h2*z21;% 求出 k 的值d1=y1(k)-h1*c10; c11=c10+k1*d1;%求出theta的估计值%求出p的值%求出两次估计值之间的相对误差%保存相对误差%保存估计值%将本次估计值作为下次估计的初值d2=y2(k)-h2*c20; c21=c20+k2*d2; p11=p10-k1*k1*h1*p10*h1+1;p21=p20-k2*k2*h2*p20*h2+1;e11=(c11-c10)./c10;e21=(c21-c20)./c20;e1(:,k)=e11;e2(:,k)=e21;c1(:,k)=c11;c2(:,k)=c21;c10=c11;c20=c21;p1

9、0=p11;p20=p21;if (e11E)&(e21E) j1=k,break;endenda11=c1(1,:);a12=c1(2,:);b11=c1(3,:);b12=c1(4,:);a21=c2(1,:);a22=c2(2,:);b21=c2(3,:);b22=c2(4,:);ea11=e1(1,:);ea12=e1(2,:);eb11=e1(3,:);eb12=e1(4,:);ea21=e2(1,:);ea22=e2(2,:);eb21=e2(3,:);eb22=e2(4,:);figure% 绘图i=1:j1;plot(i,a11(i),r,i,a12(i),b,i,b11(i)

10、,m,i,b12(i),k);title(NSR=0.4时白噪声影响下各参数辨识结果随递推次数的变化)；xlabel(递推次数);ylabel(各参数辨识结果)；axis tight; grid on;figureplot(i,a21(i),r,i,a22(i),b,i,b21(i),m,i,b22(i),k);title(NSR=0.4时有色噪声影响下各参数辨识结果随递推次数的变化);xlabel(递推次数);ylabel(各参数辨识结果)；axis tight; grid on;figureplot(i,ea11(i),r,i,ea12(i),b,i,eb11(i),m,i,eb12(i),k);title(NSR=0.4时白噪声影响下各参数辨识误差随递推次数的变化)；xlabel(递推次数);ylabel (各参数辨识误差)；axis ti