DSP

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

2019-07-13 18:44发布

站内文章 / DSP
12091 0 1264

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');           

热门文章