RLS实现求解最小二乘确定性正则方程

RLS算法：递归最小二乘算法，使用迭代的方法求解最小二乘的确定性正则方程（奇异值分解也可以）。在用RLS时，数据矩阵会进一步得到扩充，将每次观测数据都放进来，而不是像奇异值分解只将每个滤波器都有数据时作为开始。

这里就不多写RLS的推导过程；大致写一些思路：

• 根据扩展后的观测数据矩阵A，定义出随时间和样本两个变化的时间相关矩阵和时间互相关矩阵，代入确定性正则方程；考虑到离当前时刻近的观测值对相关矩阵和互相关向量的影响较大，因此对时间相关矩阵和时间互相关矩阵乘上遗忘因子，为了让时间互相关矩阵可逆，又对其进行对角加载（加上一个对角矩阵），然后将其抽离出n-1时刻对应的形式。同理，时间互相关矩阵也抽离出n-1时刻。
• 利用抽离出的时间互相关矩阵递推形式，和矩阵求逆引理进行对比，写出A逆的形式。
• 引入增益向量k(n)，它是互相关矩阵的逆对输入向量u(n)的线性变换。
• 最后写出权向量的递推公式：

                        
        这里的大括号内的一项又定义为先验估计误差。LMS是考虑的后验估计误差。

以一阶AR模型u(n)=-0.99u(n-1)+v(n)为例；使用M=2的FIR滤波器，用RLS算法实现u(n)的线性预测，即求解权系数。

p2=1;
clc,clear all,close all
for exNum=1:500
%% 产生样本序列
%由题目得差分方程为 u(n)=-0.99u(n-1)+v(n)
% 写出差分方程初值条件以及 num 和 den
den=[1 0.99]; %差分左边
num=[1]; %差分右边
u0=zeros(length(den)-1,1); %初始值
vsigma=0.995;N=1000;
v=normrnd(0,vsigma,N,1);
% v=sqrt(vsigma)*rand(N,1);
% 生成初始条件
Zi=filtic(num,den,u0); %v(n)每一次已经固定了，所以初始条件没有
un=filter(num,den,v,Zi);%这里的 v 的信息需包括样本序列个数
% plot_n=1:N;
% figure(1),plot(plot_n,v,'R-',plot_n,un,'b--');
% xlabel('n');ylabel('y(n)')
% legend('加性白噪声 v ','样本序列 u(n)');
%% 产生期望相应信号和观测数据矩阵
n0=1;M=2;
b=un(n0+1:N);
L=length(b);
un1=[zeros(M-1,1).',un.'].';%扩展数据
A=zeros(M,L);
for k=1:L
A(:,k)=un1(M-1+k:-1:k);
end
%% RLS 求解
delta=0.004;
lambda=0.98;w=zeros(M,L+1);
epsilon=zeros(L,1);
P1=eye(M)/delta;
for k=1:L
PIn=P1*A(:,k);
denok=lambda+A(:,k)'*PIn;
kn=PIn/denok;
epsilon(k)=b(k)-w(:,k)'*A(:,k);
w(:,k+1)=w(:,k)+kn*conj(epsilon(k));P1=P1/lambda-kn*A(:,k)'*P1/lambda;
end
% figure,plot(1:N,w(1,:),'b-',1:N,w(2,:),'r--')
% xlabel('n');ylabel('w')'
% legend('w1','w2');
MSE=abs(epsilon).^2;
MSEsum(:,:,exNum)=MSE;
wSum(:,:,exNum)=w;
end
wmean=mean(wSum,3);
MSEmean=mean(MSEsum,3);
figure,plot(1:N-1,MSEmean),axis([0,600,0,inf]),xlabel('n');ylabel('MSE')
figure,plot(1:N,wmean(1,:),'b--',1:N,wmean(2,:),'r-'),axis([0,600,-2.5,0.5])
xlabel('n');ylabel('w')
legend('w1','w2');