最小二乘递推算法实例应用

1.1 问题描述

1.2 方法思路

模型类选CAR模型

递推算法

1.3 实验结果

1.4 实验代码（Matlab）

main.m

%%%%采用递推最小二乘法对四个真值[a1 a2 b1 b2]=[-1.5 0.7 1 0.5]进行估计
%%%第一步：选用的模型类为CAR模型，由已给的条件可知模型阶数为2
%产生M序列做为输入，周期为1023，幅值为1
u=F2(1023,1);%%利用作业2中代码封装成的函数产生周期为1023的M序列
%产生均值为0，方差为1的高斯白噪声
V=F1(0.1,0,550);%%利用作业1中代码封装成的函数产生白噪声
%四个参数估计值
Y=zeros(4,501);
%任意给定Y初值
Y(:,1:2)=[1 2 -1 3;0.8 -2 1 3]';
%给定P（0）满足课本式5-2-19，使P（0）尽量大
P=5000*eye(4); 
%给定增益矩阵初值
K=[1 1 1 1]';
%任意给定Z初值
Z(1:2)=[0 1];
%根据课本式5-1-1计算z(k)
for k=3:502
Z(k)=1.5*Z(k-1)-0.7*Z(k-2)+U(k-1)+0.5*U(k-2)+V(k-2);
end
%%%第二步：由课本式5-2-16，另加权因子为1，得到普通最小二乘递推算法RLS，其递推公式如下：
for k=3:502   %递推500次    
%%由课本式5-1-1可以得到h(k)的表达式，模型阶数为2，na=nb=2
h=[-Z(k-1) -Z(k-2) U(k-1) U(k-2)]';%此时的h为第k次的值
%%由课本式5-2-16可得RLS表达式如下：
K=P*h*inv(h'*P*h+1);%此时的P为第k-1次的值,h为第k次的值
Y(:,k)=Y(:,k-1)+K*(Z(k)-h'*Y(:,k-1));%此时的K为第k次的值,h为第k次的值
P=(eye(4)-K*h')*P;%此时的P为第k次的值
end
%%%第三步：画图，并比较递推500次后的估计值与真值
figure
t=[1:1:502];
%%画出a1的估计值
plot(t,Y(1,:),'r')
e1=-1.5-Y(1,502);%a1与真值的误差
text(502,Y(1,502),'a1');
hold on
%%画出a2的估计值
plot(t,Y(2,:),'g')
e2=0.7-Y(2,502);%a2与真值的误差
text(502,Y(2,502),'a2');
hold on
%%画出b1的估计值
plot(t,Y(3,:),'b')
e3=1-Y(3,502);%b1与真值的误差
text(502,Y(3,502),'b1');
hold on
%%画出b2的估计值
plot(t,Y(4,:),'m')
e4=0.5-Y(4,502);%b1与真值的误差
text(502,Y(4,502),'b2');
title('递推500次后的参数估计值的跟踪曲线')
legend('a1','a2','b1','b2') 

F1.m

function Y3=F1(v,m,n)
%利用乘同余法和变换抽样法产生均值为m，方差为v的正态分布的白噪声
%公式：Xi=A*Xi*mod(M)
M1=2^11;%为2的方幂
M2=2^17;%为2的方幂
A1=119;%%不能太小
A2=279;%%不能太小
x1=1;%伪随机序列1的初值
x2=11;%伪随机序列2的初值
X1=x1;
X2=x2;
Y1=[];%伪随机序列1
Y2=[];%伪随机序列1
% n=700; %伪随机序列长度
%%%%%采用乘同余法产生伪随机序列%%%%%
for i=1:n
    %产生伪随机序列1
    X1=mod(A1*X1,M1);
   Y1=[Y1 X1/M1];
    %产生伪随机序列2
    X2=mod(A2*X2,M2);
    Y2=[Y2 X2/M2];
end
%正态分布的白噪声
Y3=m+v*sqrt(-2*log(Y1)).*cos(2*pi*Y2);%正态分布的白噪声
figure
plot(Y3,'b');
title('正态分布的白噪声');
end

F2.m

function X=F2(T,a)
%%%采用伪随机信号发生器生成周期为1023的M序列
y=[1 1 0 1 1 0 1 0 1 0];%%%初始化10级移位寄存器，保证其初值不全为0
for i = 1:1:T
    Y(i)=y(10);%%%第n级的输出即是M序列的两种状态
    M(i)=mod(y(10)+y(7),2);%第n级与第N1级模二加结果，此时分别为10和7
    %%%%移位寄存器工作原理%%%
    y(10)=y(9);
    y(9)=y(8);
    y(8)=y(7);
    y(7)=y(6);
    y(6)=y(5);
    y(5)=y(4);
    y(4)=y(3);
    y(3)=y(2);
    y(2)=y(1);
    y(1)= M(i);
end
for i=1:1:T
    if Y(i)==1
        X(i)=a;
    else
        X(i)=-a;
    end
end
stairs(X,'b'); %画出伪随机信号
title('周期为1023、幅值为1的伪随机信号');
ylim([-1.2 1.2])
end