rels.m

最小二乘递推算法的Matlab仿真 v(k)是服从N(0,1)分布的不相关随机噪声

u0=idinput(15,'prbs');
u0=u0';
u=[];
for i=1:4
    u=[u,u0];
end                        %取长度为60的M序列
figure(1);
subplot(211);
stem(u),grid on  %画出M序列输入信号
v=randn(1,60);%产生一组60个正态分布的随机噪声
subplot(212);
plot(v),grid on %画出随机噪声信号
u,v  %显示输入信号和噪声信号
%建模
a1=-1.5;a2=0.7;b1=1;b2=0.5;d1=-1;d2=0.2;
z=zeros(1,60);
z(1)=v(1);
z(2)=-a1*z(1)+b1*u(1)+v(2)+d1*v(1);
for k=3:60
    z(k)=-a1*z(k-1)-a2*z(k-2)+b1*u(k-1)+b2*u(k-2)+v(k)+d1*v(k-1)+d2*v(k-2);
end
%辨识
theta0=[0.001 0.001 0.001 0.001 0.001 0.001]';
p0=10^6*eye(6,6);
%E=0.00000000005;%相对误差E=0.00000000005
theta=[theta0,zeros(6,14)];
%e=zeros(6,15);
for k=3:60
    h1=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k-1),v(k-2)]';
    x=h1'*p0*h1+1;
    x1=inv(x);
    k1=p0*h1*x1;%求出K值
    m1=z(k)-h1'*theta0;
    theta1=theta0+k1*m1;%辨识参数theta
    v(k)=z(k)-h1'*theta1;
   % e1=theta1-theta0;
    %e2=e1./theta0;
   % e(:,k)=e2;
    theta0=theta1;
    theta(:,k)=theta1;%把递推出的辨识参数C的列向量加入辨识参数矩阵
    p1=p0-k1*h1'*p0;%find p(k)
    p0=p1;%给下次使用
   % if e2<E break;%收敛满足要求，中止计算
   % end  %判断结束
end      %循环结束
theta,v  %显示被辨识参数及收敛情况
%分离变量
a1=theta(1,:);a2=theta(2,:);b1=theta(3,:);b2=theta(4,:);
d1=theta(5,:);d2=theta(6,:);
%ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:)%分离出a1,a2,b1,b2的收敛情况
%ed1=e(5,:);ed2=e(6,:);%分离出d1,d2的收敛情况
%figure(2);
%for i=1:60
%plot(i,a1,'r',i,a2,'r:',i,b1,'b',i,b2,'b:',i,d1,'g',i,d2,'g:')%画出各个被辨识参数
%end
%title('Parameter Identificantion with Receursive Least Squares Method')%标题
%figure(3);i=1:60;
%plot(i,ea1,'r',i,ea2,'r:',i,eb1,'b',i,eb2,'b:',i,ed1,'g',i,ed2,'g:')%画出各个参数收敛情况
%title('Identification Precision')%标题