实验二递推最小二乘估计(RLS)及模型阶次辨识(F-Test）1实验方案设计1.1生成输入数据和噪声用M序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。生成白噪声时，首先利用乘同余法生成U[0,1]均匀分布的随机数，再利用U[0,1]均匀分布的随机数生成标准正态分布的白噪声。1.2过程仿真辨识模型的形式取，为方便起见，取即用M序列作为辨识的输入信号。1.3递推遗忘因子法数据长度L取534，初值1.4计算损失函数、噪声标准差损失函数噪声标准差1.6F-Test定阶法计算模型阶次统计量其中，为相应阶次下的损失函数值，为所用的数据长度，为模型的估计阶次。若，拒绝，若，接受，其中为风险水平下的阀值。这时模型的阶次估计值可取。1.6计算噪信比和性能指标噪信比参数估计平方相对偏差参数估计平方根偏差2编程说明M序列中，M序列循环周期取，时钟节拍=1Sec，幅度，特征多项式为。白噪声循环周期为。采样时间设为1Sec，。3源程序清单3.1正态分布白噪声生成函数functionv=noise(N)%生成正态分布N(0,sigma)%生成N个[01]均匀分布随机数A=179;x0=11;M=2^15;fork=1:Nx2=A*x0;x1=mod(x2,M);v1=x1/(M+1);v(:,k)=v1;x0=x1;endaipi=v;sigma=1;%标准差fork=1:length(aipi)ksai=0;fori=1:12temp=mod(i+k,length(aipi))+1;ksai=ksai+aipi(temp);endv(k)=sigma*(ksai-6);endend3.2M序列生成函数function[NprM]=createM(n,a)%生成长度为n的M序列,周期为Np，周期数为rx=[1111];%初始化初态fori=1:ny=x;x(2:4)=y(1:3);x(1)=xor(y(1),y(4));U(i)=2*y(4)-a;endM=U*a;lenx=length(x);Np=2^lenx-1;r=n/Np;end3.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);%初始Pe=0.0001;stop=1;%误差要求，循环停止信号n=N;Error=zeros(n_max,1);while(stop==1&&n<=n_max)H=[];%新的数据向量fori=1:naH=[H;-Z(n-i)];endforj=1:nbH=[H;U(n-j)];endK=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);%计算下次递推用到的Ptemp=abs((X-X_past)./X_past);%相对误差stop=sum(temp)>=e;%判断精度Error(n)=Z(n)-H'*X;n=n+1;endAes=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[Npru]=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;fork=3:L;%循环变量从3到Ly(k)=-aa1*z(k-1)-aa2*z(k-2)+bb1*u(k-1)+bb2*u(k-2);z(k)=y(k)+lamat*v(k);%给出辨识输出采样信号endov=fangcha(v);%计算噪声方差oy=fangch