rgls.m

广义最小二乘辨识（噪信比大的时候会收敛到局部极值

clear all;
clc;
display('广义最小二乘在线辨识 20081115');

y1=1;y2=0;y3=1;y4=0;
u=zeros(1,500);
for k=1:50
    if y4==0
        u(k)=1;
    else
        u(k)=-1;
    end
    x1=xor(y3,y4); %计算下一时刻的状态，y是当前输出
    x2=y1;         %计算下一时刻的状态
    x3=y2;         %计算下一时刻的状态
    x4=y3;         %计算下一时刻的状态 
    y1=x1;         %计算下一时刻的当前输出
    y2=x2;
    y3=x3;
    y4=x4;
end   
%生成M序列

v=0.03*randn(1,500);
%白噪声序列

y=zeros(1,500);
y(1)=0;y(2)=0;
for k=3:500
y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-1)+0.5*u(k-2)+v(k);  %实际系统输出,其中C(q-1)=1
end
%得到系统实际输出序列（加了噪声）

N=0.000000005;
p0=10^6*eye(4,4);  %计算模型theta的P矩阵初值
pe0=10^6*eye(2,2); %计算滤波器thetae的P矩阵初值
yf=zeros(1,500);
yf(1)=0;yf(2)=0;
uf=zeros(1,500);
uf(1)=0;uf(2)=0;
e=zeros(1,500);   %残差序列
e(1)=0;e(2)=0;   
theta=zeros(4,500);  %模型参数每次辨识的结果
thetae=zeros(2,500); %滤波器参数每次辨识的结果
theta0=[0.001;0.001;0.001;0.001]; %系统参数初值
thetae0=[0.001;0.001]; %滤波器C(q-1)的参数初值
for k=3:500
    yf(k)=y(k)+thetae0(1)*y(k-1)+thetae0(2)*y(k-2);  %滤波后的输入序列
    uf(k)=u(k)+thetae0(1)*u(k-1)+thetae0(2)*u(k-2);  %滤波后的输出序列
    h1=[-yf(k-1);-yf(k-2);uf(k-1);uf(k-2)]; %滤波后的输出输出组成的辨识数据
    x=1+h1'*p0*h1;
    x1=inv(x);
    k1=p0*h1*x1;  %求当前时刻的K矩阵
    d1=yf(k)-h1'*theta0;
    theta1=theta0+k1*d1;
    theta(:,k)=theta1;
    e1=theta1-theta0;  
    e2=e1./theta0;   
    E(:,k)=e1;   %保存模型参数的误差
    theta0=theta1; %为下时刻准备
    p1=p0-k1*k1'*[1+h1'*p0*h1]; %为下时刻准备
    p0=p1;  %为下时刻准备
    %%%递推算法的实质就是theta和P在原来的基础上计算，因此要为下时刻的这两个值做准备
    
    e(k)=y(k)-[-y(k-1),-y(k-2),u(k-1),u(k-2)]*theta1;  %计算残差e(k),注意这里的h的组成
    
    he1=[-e(k-1);-e(k-2)];
    xe=1+he1'*pe0*he1;
    xe1=inv(xe);
    ke1=pe0*he1*xe1;  %求当前的K矩阵
    de1=e(k)-he1'*thetae0;
    thetae1=thetae0+ke1*de1;
    thetae(:,k)=thetae1;  %保存滤波器参数
    ee1=thetae1-thetae0;
    ee2=ee1./thetae0;      %模型参数的相对误差
    Ee(:,k)=ee1;
    thetae0=thetae1;   %为下一时刻做准备
    pe1=pe0-ke1*ke1'*[he1'*pe0*he1+1]; %为下一时刻赋值
    pe0=pe1; %为下一时刻P矩阵赋值
    %%%递推算法的实质就是theta和P在原来的基础上计算，因此要为下时刻的这两个值做准备
end

display('->Done!');