RLS自适应滤波算法

一. RLS(递归最小二乘)算法介绍

RLS(Recursive Least Squares,递归最小二乘)算法是一种自适应滤波算法,用于通过递归方式求解最小二乘问题。它是传统最小二乘法的在线学习版本,能够实时更新系统参数估计。

1.1 基本概念

RLS算法的主要特点:
- 递归计算:新数据到来时,不需要重新计算全部数据,只需基于前一时刻的估计进行更新
- 快速收敛:通常比LMS(最小均方)算法收敛更快
- 计算复杂度较高:相比LMS算法需要更多的计算资源

1. 2最小二乘问题回顾
标准最小二乘问题的目标是最小化误差平方和:

其中:
- w是待估计的参数向量(滤波器系数)
- u(i) 是第 i 个输入向量(回归变量)
-  d(i) 是期望输出(观测数据)
- n 是数据点数

在批处理最小二乘中,解为:

但批处理方法计算量大,且无法适应在线数据流。RLS 采用递归更新方式,逐样本更新估计值。

二. RLS算法推导

1.问题背景
给定线性模型:

其中:
- y(n) :第 n 时刻的观测值
- x(n) :n 时刻的输入向量( N * 1 )
- h(n):待估计的参数向量( N * 1 )
- v(n) :噪声

2.目标:通过递归方式最小化加权误差平方和:

       

       λ:遗忘因子(0<λ≤1),用于降低旧数据的影响       

  3.定义自相关与互相关矩阵       

      自相关矩阵:              

                                        

     互相关矩阵:

                                       

4:最小二乘解
       最优参数估计:

                                        

5:递归更新
利用矩阵求逆引理(Sherman-Morrison-Woodbury)实现递归更新:

                     

三.. RLS算法流程

        

四. 关键参数与特性

  • 遗忘因子λλ

    • λ≈1:长记忆(适用于稳态系统)。

    • λ<1:快速跟踪时变系统(典型值0.95-0.99)。

  • 计算复杂度

    • O(N2)每次迭代(因需更新P(n))。

    • 比LMS算法(O(N))更高,但收敛更快。

  • 优点

    • 指数收敛,无梯度噪声。

    • 适用于非平稳信号(如通信信道均衡)。

  • 缺点

    • 数值稳定性需注意(可通过平方根RLS改进)。

    • 高计算成本。                                               

五.代码实现(与LMS算法对比)

                                                      

clear all;close all;
k=5;%表示参数的个数
L=2000;%表示输入长度
w0=randn(k,1);%待估计的真实参数
Uk=randn(k,L);%更新输入的值
Vc=randn(1,L)*0.1;%噪声,0.1表示噪声的强度
Dk=w0'*Uk+Vc;%表示观测值
wk_rls=zeros(k,1);%初始化参数估计
P_rls=eye(k);%初始化逆自相关矩阵
lmd=1;%表示遗忘因子
%RLS
for n=1:L
    dn=Dk(n);%表示第n个观测值
    un=Uk(:,n);%表示n时刻的输入值,N*1
    en=dn-wk_rls'*un;%误差
    kn=P_rls*un/(lmd+un'*P_rls*un);%增益向量
    P_rls=1/lmd*(P_rls-kn*un'*P_rls);%更新逆自相关矩阵
    Err_rls(n)=norm(w0-wk_rls);%计算参数误差,欧氏距离
    wk_rls=wk_rls+kn*en;%更新参数
end
%LMS
wk_lms=zeros(k,1);%参数初始化
mu=0.005;%步长,需要调试
for n=1:L
    dn=Dk(n);%表示第k个观测值
    un=Uk(:,n);%表示k时刻的输入值,n*1
    en=dn-wk_lms'*un;%误差
    Err_lms(n)=norm(w0-wk_lms);
    wk_lms=wk_lms+mu*en*un;
end
figure,hold on;
grid on;
plot(log(Err_rls));
plot(log(Err_lms));
legend('RLS','LMS');

很明显,RLS收敛的更快。

用于滤波,rls自适应滤波 #include #include #include #define N 110000 float x[N]; float x1[110050]; float xout[N]; float w[50]; float p[50][50]; float u[50]; float kn[50]; float kn1[50]; float dn; float e; float outmiddle; //阶数50// //****************************************初始化程序用于计算w[][].p[]**********************************// float initial () { char i1,j1; for(i1=0;i1<50;i1++) { w[i1]=0; } for(i1=0;i1<50;i1++) { for(j1=0;j1<50;j1++) { p[i1][j1]=0; } } for(i1=0;i1<50;i1++) { p[i1][i1]=9999999.9; } return(0); } //****************************************************************************************************************// //****************************************e[n]的求解**************************************************************// float en(unsigned long int datanumber) { /*float result=0; int i2,k2; int j2; j2=datanumber; k2=datanumber+5; e=0; for(i2=0;i2<50;i2++)//计算w()*u() { e=e+w[i2]*u[i2]; } dn=0; for(i2=j2;i2<k2;i2++) { dn=dn+x[i2]; } dn=dn/5;//给定dn为0.2 e=dn-e; return(0);*/ } //****************************************************************************************************************// //****************************************k[n]的求解**************************************************************// void knkn() { int i3,j3; float knmiddle=0; float middle[50]; float pluse=0; for(i3=0;i3<50;i3++) { for(j3=0;j3<50;j3++) { knmiddle=knmiddle+p[i3][j3]*u[j3]; } kn[i3]=knmiddle; knmiddle=0; }//计算p(n-1)*u(n) for(i3=0;i3<50;i3++) { middle[i3]=kn[i3]; } for(i3=0;i3<50;i3++) { pluse=pluse+middle[i3]*u[i3]; }//计算uh(n)*p(n-1)*u(n) pluse=pluse+0.7;//衰减因数0.7 for(i3=0;i3<50;i3++) { kn[i3]=kn[i3]/pluse; } pluse=0; } //****************************************************************************************************************// //*******************************************用于pn[]的更新*******************************************************// float newpn() { float newpn[50][50]; float newpn1[50][50]; float newpn3[50][50]; float newpn2=0; int i4,j4,k4; for(i4=0;i4<50;i4++) { for(j4=0;j4<50;j4++) { newpn[i4][j4]=kn[i4]*u[j4]; newpn3[i4][j4]=p[i4][j4]; } }//计算k(n)*u(n) for(i4=0;i4<50;i4++) { for(j4=0;j4<50;j4++) { for(k4=0;k4<50;k4++) { newpn2=newpn2+newpn[i4][k4]*newpn3[k4][j4]; } newpn1[i4][j4]=newpn2; newpn2=0; } }//计算k()*u()*p() for(i4=0;i4<50;i4++) { for(j4=0;j4<50;j4++) { p[i4][j4]=p[i4][j4]-newpn1[i4][j4]; p[i4][j4]=p[i4][j4]/0.7;//衰减因子是0.7 } } return(0); } //*****************************************************************************************************************// //*************************************用于wn[]的更新**************************************************************// float newwn() { char i5,j5; for(i5=0;i5<50;i5++) { kn1[i5]=kn[i5]; } for(i5=0;i5<50;i5++) { kn1[i5]=kn1[i5]*e; w[i5]=w[i5]+kn1[i5]; } return(0); } main() { //, f[M][N], g[M][N], a[M][M], P[N], K[M]; //float x1[200]; //float Knum = 0, Kden = 0, sum = 0, eP; // int i, j, k, m, n; float dn; long int i=0, j=0,k=0,m=0,n1,n2,n3; float result=0; int i2,k2; int j2; char controller=1; FILE * fpr, * fp_out; int max[2][2]={1,2,3,4}; int max1[2]={1,2}; int max2[2]={2,3}; int maxmiddle; //a[0][1] = 1; if((fpr = fopen("wavewavedata.bin","rb")) == NULL) {printf("cannot open file\n"); exit(0); } for(i=0;i<N;i++) { fread(&x[i], 4, 1, fpr); //printf("%f\n",x[i]); } /*for(i=0;i<49;i++) { x1[i]=x[110000-i-1]; } for(i=49;i<110049;i++) { x1[i]=x[i-49]; }//得到要处理的110049个数据*/ fclose(fpr); for(i=0;i<N;i++) { xout[i]=x[i]; } if((fp_out = fopen("bbboutdata.bin","w")) == NULL) { printf("cannot open file\n"); return; } else printf("success fp_out = %d\n",fp_out); for(j=0; j<N; j++) { fwrite(&xout[j],4, 1, fp_out); //fread(&x1[i][j],4, 1, fpr); //printf("%f\n",xout[j]); } fclose(fp_out); getchar(); initial (); for(m=0;m<N;m++) { //m=0; n1=m; n2=m+50; k=49; for(i=n1;i<n2;i++) //取出要处理的50个数据放于u[n] { u[k]=x1[i]; k--; } //printf("%f\n",u[0]); while(controller) { //en(n3); j2=n2-3; k2=n2+2; e=0; for(i2=0;i2<50;i2++)//计算w()*u() { e=e+w[i2]*u[i2]; } dn=0; for(i2=j2;i2<k2;i2++) { dn=dn+x1[i2]; } dn=dn/5;//给定dn为0.2 e=dn-e; if(e<0.000000001) { controller=0; /*printf("success!",m); for(i=0;i<50;i++) { printf("%f\n",w[i]); } printf("success!",m); for(i=0;i<50;i++) { printf("%f\n",kn[i]); }*/ } else { knkn(); newpn(); newwn(); } } /*for(j=0;j<50;j++) { printf("%f\n",w[j]); }*/ outmiddle=0; for(i=0;i<50;i++) { outmiddle=outmiddle+w[i]*u[i]; } xout[m]=outmiddle; //printf("output",m); //printf("%f\n",xout[m]); } /*for(i=0;i<N;i++) { //xout[i]=x[i]; printf("%f\n",xout[i]); }*/ /*if((fp_out = fopen("bboutdata.bin","w")) == NULL) { printf("cannot open file\n"); return; } else printf("success fp_out = %d\n",fp_out); for(j=0; j<N; j++) { fwrite(&xout[j],4, 1, fp_out); //fread(&x1[i][j],4, 1, fpr); //printf("%f\n",xout[j]); } fclose(fp_out); getchar();*/ }
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值