递推最小二乘法及模型阶次辨识

close all; clear all; %%%%%%%%%%%人机对话%%%%%%%%%%%%%%%%% lambda=input('噪声标准差为（0.0或0.1或0.5或1.0）? '); L=input('数据长度为（100或300或500）? '); miu=input('遗忘因子为（0,1之间）? '); %%%%%%%%设定模型阶次从2到6，即Nbeg=2,Nend=6(假定na=nb)%%%%%%%%% %%%%%%%%%%%%%%%赋初值并生成M序列和噪声序列%%%%%%%%%% %%%%%%%%%%%%%%生成a=1,P=4的M序列%%%%%%%%%%%%%%%%% M=[0 1 0 1 0 ]; Np=2^4-1; for k=1:534 if M(5)==0 u(k)=1;end; if M(5)==1 u(k)=-1;end; for i=5:-1:2 M(i)=M(i-1); end M(1)=M(2)+M(5); if M(1)==2 M(1)=0;end; end %%%%%%%%%%%%%%%生成白噪声v（k）%%%%%%%%%%%%%%%%%%%% A=179; xi=11; M=32768; for k=1:534 ksai=0; for i=1:12 xi=A*xi; xi=mod(xi,M); ksai=ksai+(xi/M); end v(k)=ksai-6; end %%%%%%%%%%%%%%%%%过程仿真%%%%%%%%%%%%%%%% z(1)=0; z(2)=0; for k=3:length(u) z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+lambda*v(k);% end %%%%%%%%%%%%%%%%%%模型参数估计，RLS算法%%%%%%%%%%%% for n=2:6 for i=1:2*n theta(i,n)=eps; end P=eye(2*n); for i=1:L for j=1:n if (i-j)<=0 H(j,i)=0; else H(j,i)=-z(i-j); end end for m=(n+1):2*n if (i-m+n)<=0 H(m,i)=0; else H(m,i)=u(i-m+n); end end end J(n)=0; for k=1:L K=P*H(:,k)*inv(H(:,k)'*P*H(:,k)+miu); J(n)=miu*(J(n)+((z(k)-H(:,k)'*theta(:,n))^2/(H(:,k)'*P*H(:,k)+miu))); theta(:,n)=theta(:,n)+K*(z(k)-H(:,k)'*theta(:,n)); P=((eye(2*n)-K*H(:,k)')*P)/miu; end end %%%%%%%%%%%%%%%模型阶次识别%%%%%%%%%%%%%%%% if L==100 for n=1:5 if (J(n)-J(n+1))*(L-2*n-2)/(2*J(n+1))<=3.09 N=n+1 break end end end if L==300 for n=1:5 if (J(n)-J(n+1))*(L-2*n-2)/(2*J(n+1))<=3.03 N=n+1 break end end end if L==500 for n=1:5 if (J(n)-J(n+1))*(L-2*n-2)/(2*J(n+1))<=3.01 N=n+1 break end end end %%%%%%%%%%%%%%%输出结果%%%%%%%%%%%%%%%% if N==2 A=theta(1:2,N) B=theta(3:4,N) end if N==3 A=theta(1:3,N) B=theta(4:6,N) end if N==4 A=theta(1:4,N) B=theta(5:8,N) end if N==5 A=theta(1:5,N) B=theta(6:10,N) end if N==6 A=theta(1:6,N) B=theta(7:12,N) end %%%%%%%%%%%%%%%%%%%计算性能指标%%%%%%%%%%%%%%%%%% lam1=sqrt(J(N)/(L-2*N))%%%%%%%%%%%%%%噪声标准差的估计 lam2=sum(B)/(1+sum(A))%%%%%%%%%%%%%模型静态增益估计 a0=1;a1=-1.5;a2=0.7;b0=0;b1=1;b2=0.5; thetazhen=[a1,a2,b1,b2]'; chazhi=thetazhen-theta(1:4,2); lam3=sqrt(sum((chazhi./thetazhen).^2))%%%%%%参数估计平方相对偏差 lam4=sqrt(sum(chazhi.^2)/sum(thetazhen.^2))%%%%%%%%%参数估计平方根偏差 Kzhen=(1+0.5)/(1-1.5+0.7); lam5=sqrt(abs(lam2-Kzhen)/Kzhen)%%%%%%%%%静态增益估计相对偏差 a02=a0;a12=a1;a22=a2; a01=(a02*a02-a22*a22)/a02; a11=(a02*a12-a22*a12)/a02; a00=(a01*a01-a11*a11)/a01; b02=b0;b12=b1;b22=b2; b11=(a02*b12-b22*a12)/a02; b01=(a02*b02-b22*a22)/a02; b00=(a01*b01-b11*a11)/a01; vary=b00^2/a00+b11^2/a01+b22^2/a02; xinzaobi=lambda/sqrt(vary)%%%%%%%%%%%计算信噪比