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

实验二递推最小二乘估计(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:N
x2=A*x0;
x1=mod(x2,M);
v1=x1/(M+1);
v(:,k)=v1;
x0=x1;
end
aipi=v;
sigma=1;%标准差
fork=1:length(aipi)
ksai=0;
fori=1:12
temp=mod(i+k,length(aipi))+1;
ksai=ksai+aipi(temp);
end
v(k)=sigma*(ksai-6);
end
end
3.2M序列生成函数
function[NprM]=createM(n,a)
%生成长度为n的M序列,周期为Np，周期数为r
x=[1111];%初始化初态
fori=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==1&&n<=n_max)
H=[];%新的数据向量
fori=1:na
H=[H;-Z(n-i)];
end
forj=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
[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到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=fangch