递推最小二乘估计及模型阶次辨识.doc

实验 递推最小二乘估计(RLS) 及模型阶次辨识(F-Test） 1.1 生成输入数据和噪声 用M序列作为辨识的输入信号即 用M序列作为辨识的输入信号 3.1 正态分布白噪声生成函数 function v=noise(N) %生成正态分布N(0,sigma) %生成N个[0 1]均匀分布随机数 A=179; x0=11; M=2^15; for k=1:N x2=A*x0; x1=mod(x2,M); v1=x1/(M+1); v(:,k)=v1; x0=x1; end aipi=v; sigma=1; %标准差 for k=1:length(aipi) ksai=0; for i=1:12 temp=mod(i+k,length(aipi))+1; ksai=ksai+aipi(temp); end v(k)=sigma*(ksai-6); end end 3.2 M序列生成函数 function [Np r M]=createM(n,a) %生成长度为n的M序列,周期为Np，周期数为r x=[1 1 1 1]; %初始化初态 for i=1:n y=x; x(2:4)=y(1:3); x(1)=xor(y(1),y(4)); U(i)=2*y(4)-a; end M=U*a; lenx=length(x); Np=2^lenx-1; r=n/Np; end 3.3 加权最小二乘递推算法函数 function [Aes,Bes,Error]=RLS(na,nb,Z,U,f) %Aes、Bes为参数估计值，na、nb为模型阶次，Z、U为输出输入数据，f为加权因子 N=na+nb; n_max=length(Z); X=0.001.*ones(N,1); %初始估计值 P=10^5.*eye(N); %初始P e=0.0001; stop=1; %误差要求，循环停止信号 n=N; Error=zeros(n_max,1); while(stop==1n=n_max) H=[]; %新的数据向量 for i=1:na H=[H;-Z(n-i)]; end for j=1:nb H=[H;U(n-j)]; end K=P*H*inv(H*P*H+f); %计算增益矩阵 X_past=X; X=X+K*(Z(n)-H*X); %计算新的估计值 P=P-K*K*(H*P*H+f); %计算下次递推用到的P temp=abs((X-X_past)./X_past); %相对误差 stop=sum(temp)=e; %判断精度 Error(n)=Z(n)-H*X; n=n+1; end Aes=X(1:na); Bes=X(na+1:N); 3.4 主函数 clear %清理工作间变量 L=534; %M序列的周期,四级移位寄存器生成M序列,作为输入信号u(k) ex=60; %在图像中展示的数据个数 a=1; aa1=-1.5; aa2=0.7; bb1=1; bb2=0.5; %提前规定的a,b,c,d [Np r u]=createM(L,a); %生成M序列 figure(1); %画第1个图形：u(k) stem(u(1:ex)),grid; %以径的形式显示出部分输入信号并给图形加上网格 xlabel(k) %标注横轴变量 ylabel(输入信号) %标注纵轴变量 title([四级移位寄存器生成M序列输入信号(前,int2str(ex),位)]) %图形标题 z(2)=0;z(1)=0; %取z的前两个初始值为零 y=z; v=noise(L);% 生成白噪声 lamat=0.1; for k=3:L; %循环变量从3到L y(k)=-aa1*z(k-1)-aa2*z(k-2)+bb1*u(k-1)+bb2*u(k-2); z(k)=y(k)+lamat*v(k); %给出辨识输出采样信号 end ov=fangcha(v); %计算噪声方差 oy=fangcha(y); %计算信号方差 yita=sqrt(oy\ov); %计算噪信比 %用最小二乘递推算法辨识参数：a,b,c,d e0=[0.001 0.001 0.001 0.001]; %被辨识参数的初始值采用直接取方式，取一个充分小的实向量 p0=10^7*eye(4,4); %初始状态P0也采用直接取