fun.m

在MATLAB下

%SI3.1_RLS.m
clear all
close all
clc
miu=0.95
bata=0.95
%产生N(0,1)正态分布的随机噪声
randn('seed',100);
v=randn(1,60); 

%产生M序列
L=60;% M序列的周期
y1=1;y2=1;y3=1;y4=0;%四个移位积存器的输出初始值
for i=1:L;
x1=xor(y3,y4);
x2=y1;
    x3=y2;
    x4=y3;
    y(i)=y4;
    if y(i)>0.5,u(i)=-5;%M序列的幅值为5
    else u(i)=5;
    end
    y1=x1;y2=x2;y3=x3;y4=x4;
end 
%figure(1);
%stem(u),grid on
% 递推最小二乘辨识程序
z(2)=0;z(1)=0;
%观测值由理想输出值加噪声
for k=3:60;%循环变量从3到15   
    z(k)=1.6*z(k-1)-0.7*z(k-2)+u(k-1)+5.5*u(k-2)+v(k); 
end
%RLS递推最小二乘辨识
c0=[0.001 0.001 0.001 0.001]';
c1=[0.001 0.001 0.001 0.001]';

p0=10^3*eye(4,4);

E=0.0000000000005;%相对误差
c=[c0,c1,zeros(4,58)];%被辨识参数矩阵的初始值及大小
e=zeros(4,60);%相对误差的初始值及大小
lamt=1;
for k=3:60; 
    h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';
    p1=inv(miu*eye(4)-bata*miu*eye(4)+bata*inv(p0)+h1*h1');
    
    c2=c1+bata*miu*p1*(c1-c0)+p1*h1*(z(k)-h1'*c1);
     
    %k1=p0*h1*inv(h1'*p0*h1+1*lamt);%求出K的值
    %new=z(k)-h1'*c0; 
    %c1=c0+k1*new;%求被辨识参数c
    %p1=1/lamt*(eye(4)-k1*h1')*p0;
    e1=(c2-c1)./c1;%求参数当前值与上一次的值的差值
    e(:,k)=e1; %把当前相对变化的列向量加入误差矩阵的最后一列    
    c(:,k)=c2;%把辨识参数c 列向量加入辨识参数矩阵的最后一列 
    c0=c1;%新获得的参数作为下一次递推的旧参数
    c1=c2;
    p0=p1;
    if norm(e1)<=E 
        break;%若参数收敛满足要求，终止计算
    end
end
%分离参数
a1=c(1,:);
aa1=c(1,60)
a2=c(2,:);
aa2=c(2,60)
b1=c(3,:);
bb1=c(3,60)
b2=c(4,:);
bb2=c(4,60)
ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); 
figure(2);
i=1:60;
plot(i,a1,'k',i,a2,'b',i,b1,'r',i,b2,'g') %画出辨识结果
grid on
legend('a1','a2','b1','b2');
title('递推阻尼最小二乘参数辨识')
figure(3); 
i=1:60; 
plot(i,ea1,'k',i,ea2,'b',i,eb1,'r',i,eb2,'g') %画出辨识结果的收敛情况
grid on
legend('a1','a2','b1','b2');
title('辨识精度')