精华内容
下载资源
问答
  • RLS算法实现修改了 RLS 的代码。 最简单的代码。 有什么不对的就告诉我。 谢谢。
  • RLS算法

    万次阅读 热门讨论 2019-04-07 00:55:40
    两个算法一起说了,而且学习的话,比较轻量级。 https://zhuanlan.zhihu.com/p/67250500 一般最小二乘法 一般最小二乘法是给定若干观测值,计算一个最有可能的估计。 [ y ( 1 ) y ( 2 ) ⋮ y ( k ) ] = [ h ...

    贴一篇我在知乎的文章,这里面也阐述了最小二乘法与卡尔曼滤波的关系。两个算法一起说了,而且学习的话,比较轻量级。

    https://zhuanlan.zhihu.com/p/67250500

    一般最小二乘法

    一般最小二乘法是给定若干观测值,计算一个最有可能的估计。
    [ y ( 1 ) y ( 2 ) ⋮ y ( k ) ] = [ h ( 1 ) 1 h ( 1 ) 2 ⋯ h ( 1 ) n h ( 2 ) 1 h ( 2 ) 2 ⋯ h ( 2 ) n ⋮ ⋮ ⋱ ⋮ h ( k ) 1 h ( k ) 2 ⋯ h ( k ) n ] [ x 1 x 2 ⋮ x n ] \left [\begin{matrix}y(1)\\ y(2)\\ \vdots\\ y(k) \end{matrix}\right]=\left [\begin{matrix}h(1)_1&h(1)_2& \cdots&h(1)_n\\ h(2)_1&h(2)_2& \cdots&h(2)_n\\ \vdots&\vdots&\ddots&\vdots\\ h(k)_1&h(k)_2& \cdots&h(k)_n \end{matrix}\right]\left [\begin{matrix}x_1\\ x_2\\ \vdots\\ x_n \end{matrix}\right] y(1)y(2)y(k)=h(1)1h(2)1h(k)1h(1)2h(2)2h(k)2h(1)nh(2)nh(k)nx1x2xn
    计算上述最有可能的 x x x取值,是使得代价函数:
    J ( k ) = ( Y ( k ) − H ( k ) x ) T ( Y ( k ) − H ( k ) x ) = Y T ( k ) Y ( k ) − Y T ( k ) H ( k ) x − x T H T ( x ) Y ( k ) + x T H ( k ) T H ( k ) x J(k)=(Y(k)-H(k)x)^T(Y(k)-H(k)x)\\ =Y^T(k)Y(k)-Y^T(k)H(k)x-x^TH^T(x)Y(k)+x^TH(k)^TH(k)x J(k)=(Y(k)H(k)x)T(Y(k)H(k)x)=YT(k)Y(k)YT(k)H(k)xxTHT(x)Y(k)+xTH(k)TH(k)x
    最小的 x x x的取值。令其为 x ^ ( k ) \hat x(k) x^(k),意思是 k k k时刻的估计值
    上述中, H ( k ) H(k) H(k)代表 k k k时刻以及 k k k时刻以前的观测,即为:
    H ( k ) = [ h ( 1 ) T ,   h ( 2 ) T , ⋯ h ( k ) T ] T H(k)=[h(1)^T,\ h(2)^T,\cdots h(k)^T]^T H(k)=[h(1)T, h(2)T,h(k)T]T
    Y ( k ) Y(k) Y(k)类似:
    Y ( k ) = [ y ( 1 ) , y ( 2 ) , ⋯ y ( k ) ] T Y(k)=[y(1),y(2),\cdots y(k)]^T Y(k)=[y(1),y(2),y(k)]T
    J ( k ) J(k) J(k)最小时,有:
    ∂ J ( k ) ∂ x = − 2 Y T ( k ) H ( k ) + 2 x T H T ( k ) H ( k ) = 0 \frac{\partial J(k)}{\partial x}=-2Y^T(k)H(k)+2x^TH^T(k)H(k)=0 xJ(k)=2YT(k)H(k)+2xTHT(k)H(k)=0
    此时, x = x ^ ( k ) x = \hat x(k) x=x^(k)
    所以: x ^ ( k ) = ( H T ( k ) H ( k ) ) − H T ( k ) Y ( k ) \hat x(k)=(H^T(k)H(k))^-H^T(k)Y(k) x^(k)=(HT(k)H(k))HT(k)Y(k)
    只有 H T ( k ) H ( k ) H^T(k)H(k) HT(k)H(k) 满秩时,才存在逆矩阵。
    所以只有当 k ≥ n k\geq n kn时, H T ( k ) H ( k ) H^T(k)H(k) HT(k)H(k)才可能存在逆矩阵。

    带权最小二乘法

    有时候,测量有好有坏,带权的最小二乘法就比较有必要了。可以人为赋予每个数据的一个置信度 r ( k ) r(k) r(k),有时候也会令 r 2 ( k ) = 1 D ( v ( k ) ) = 1 σ k 2 r^2(k)=\frac{1}{D(v(k))}=\frac{1}{\sigma_k^2} r2(k)=D(v(k))1=σk21
    其中, v ( k ) v(k) v(k) 是第 k k k测量的误差。区别于估计误差。上面的式子也比较好理解。测量的不确定性会削弱这次测量的置信度。如果利用测量误差的话,则不能存在 0 0 0误差的情况。改写 J ( k ) J(k) J(k)的形式:
    J ( k ) = ∑ i = 1 k r 2 ( k ) ( y ( k ) − y ^ ( k ) ) 2 J(k)=\sum_{i=1}^kr^2(k)(y(k)-\hat y(k))^2 J(k)=i=1kr2(k)(y(k)y^(k))2
    令: R ( k ) = [ r 2 ( 1 ) 0 ⋯ 0 0 r 2 ( 2 ) ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ r 2 ( n ) ] R(k)=\left [\begin{matrix}r^2(1)&0& \cdots&0\\ 0&r^2(2)& \cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0& \cdots&r^2(n) \end{matrix}\right] R(k)=r2(1)000r2(2)000r2(n)
    那么: J ( k ) = ( Y ( k ) − H ( k ) x ^ ) T R ( Y ( k ) − H ( k ) x ^ ) = Y ( k ) T R ( k ) Y ( k ) − Y T ( k ) R ( k ) H ( k ) x ^ − x ^ T H ( k ) T R ( k ) Y ( k ) + x ^ T H T ( k ) R ( K ) H T ( k ) x ^ J(k)=\Big(Y(k)-H(k)\hat x\Big)^TR\Big(Y(k)-H(k)\hat x\Big)\\=Y(k)^TR(k)Y(k)-Y^T(k)R(k)H(k)\hat x-\hat x^TH(k)^TR(k)Y(k)+\hat x^TH^T(k)R(K)H^T(k)\hat x J(k)=(Y(k)H(k)x^)TR(Y(k)H(k)x^)=Y(k)TR(k)Y(k)YT(k)R(k)H(k)x^x^TH(k)TR(k)Y(k)+x^THT(k)R(K)HT(k)x^
    令: ∂ J ( k ) ∂ x ^ = − 2 Y T ( k ) R ( k ) + 2 x T H T ( k ) R ( k ) H ( k ) = 0 \frac{\partial J(k)}{\partial\hat x}=-2Y^T(k)R(k)+2x^TH^T(k)R(k)H(k)=0 x^J(k)=2YT(k)R(k)+2xTHT(k)R(k)H(k)=0
    则: x ^ ( k ) = ( H T ( k ) R ( k ) H ( k ) ) − H T ( k ) R ( k ) Y ( k ) \hat x(k)= (H^T(k)R(k)H(k))^-H^T(k)R(k)Y(k) x^(k)=(HT(k)R(k)H(k))HT(k)R(k)Y(k)
    这样设置权重,其实是对噪声的归一化。

    R L S RLS RLS 递推算法

    有时候,我们不可能一次获得所有的观测结果,有时候随着观测结果的增加。从新修改参数使用最小二乘法计算,会很耗费时间。 R L S RLS RLS 就是等价于一般最小二乘法的递推计算方法。

    这里介绍一个精巧的矩阵反演:
    在这只要求 A , D A,D A,D可逆的情况下给定分块矩阵:
    [ A B C D ] \left [\begin{matrix} A& B\\ C&D \end{matrix}\right] [ACBD]
    令: E = D − C A − B ,    F = A − B D − C E =D-CA^-B,\ \ F = A-BD^-C E=DCAB,  F=ABDC
    计算可得:
    [ A B C D ] [ A − + A − B E − C A − − A − B E − − E − C A − E − ] = [ I O C A − + C A − B E − C A − − D E − C A − − C A − B E − + D E − ] = [ I O C A − + ( C A − B − D ) E − C A − ( − C A − B + D ) E − ] = [ I O O I ] \left [\begin{matrix} A&B\\ C&D \end{matrix}\right]\left [\begin{matrix} A^-+A^-BE^-CA^-&-A^-BE^-\\ -E^-CA^-&E^- \end{matrix}\right]\\=\left [\begin{matrix} I&O\\ CA^-+CA^-BE^-CA^--DE^-CA^-&-CA^-BE^-+DE^- \end{matrix}\right]\\=\left [\begin{matrix} I&O\\ CA^-+(CA^-B-D)E^-CA^-&(-CA^-B+D)E^- \end{matrix}\right]=\left [\begin{matrix} I&O\\ O&I \end{matrix}\right] [ACBD][A+ABECAECAABEE]=[ICA+CABECADECAOCABE+DE]=[ICA+(CABD)ECAO(CAB+D)E]=[IOOI]
    另一个方向: [ A B C D ] [ F − − A − B E − − E − C A − E − ] = [ I O O I ] \left [\begin{matrix} A&B\\ C&D \end{matrix}\right]\left [\begin{matrix}F^-&-A^-BE^-\\-E^-CA^-&E^- \end{matrix}\right]=\left [\begin{matrix} I&O\\ O&I \end{matrix}\right] [ACBD][FECAABEE]=[IOOI]
    这个方向可以自行计算。
    这也就是说,当矩阵 A , D , E , F A,D,E,F A,D,E,F可逆时,: ( A + B D − C ) − = A − − A − B ( D + C A − B ) − C A − (A+BD^-C)^-=A^--A^-B(D+CA^-B)^-CA^- (A+BDC)=AAB(D+CAB)CA

    P ( k ) = ( H T ( k ) H ( k ) ) − P(k)=(H^T(k)H(k))^- P(k)=(HT(k)H(k))
    则有:
    P ( k ) = ( H T ( k ) H ( k ) ) − = ( H T ( k − 1 ) H ( k − 1 ) + h T ( k ) h ( k ) ) − = ( P − ( k − 1 ) + h T ( k ) h ( k ) ) − = P ( k − 1 ) − P ( k − 1 ) h T ( k ) h ( k ) P ( k − 1 ) 1 + h ( k ) P ( k − 1 ) h T ( k ) P(k)=(H^T(k)H(k))^-\\=(H^T(k-1)H(k-1)+h^T(k)h(k))^-\\=(P^-(k-1)+h^T(k)h(k))^-\\=P(k-1)-\frac{P(k-1)h^T(k)h(k)P(k-1)}{1+h(k)P(k-1)h^T(k)} P(k)=(HT(k)H(k))=(HT(k1)H(k1)+hT(k)h(k))=(P(k1)+hT(k)h(k))=P(k1)1+h(k)P(k1)hT(k)P(k1)hT(k)h(k)P(k1)
    又有: x ^ ( k ) = P ( k ) H T ( k ) Y ( k ) \hat x(k)=P(k)H^T(k)Y(k) x^(k)=P(k)HT(k)Y(k)
    可得: P − ( k ) x ^ ( k ) = H T ( k ) Y ( k ) P^-(k)\hat x(k) = H^T(k)Y(k) P(k)x^(k)=HT(k)Y(k)
    又因为: P ( k ) = ( P − ( k − 1 ) + h T ( k ) h ( k ) ) − P(k)=(P^-(k-1)+h^T(k)h(k))^- P(k)=(P(k1)+hT(k)h(k))
    得: P − ( k − 1 ) = P − ( k ) − h T ( k ) h ( k ) P^-(k-1)=P^-(k)-h^T(k)h(k) P(k1)=P(k)hT(k)h(k)
    所以: x ^ ( k ) = P ( k ) H T ( k ) Y ( k ) = P ( k ) ( P − ( k − 1 ) x ^ ( k − 1 ) + h T ( k ) y ( k ) ) = P ( k ) ( ( P − ( k ) − h T ( k ) h ( k ) ) x ^ ( k − 1 ) + h T ( k ) y ( k ) ) = x ^ ( k − 1 ) + P ( k ) h T ( k ) ( y ( k ) − h ( k ) x ^ ( k − 1 ) ) \hat x(k) = P(k)H^T(k)Y(k)\\=P(k)\Big(P^-(k-1)\hat x(k-1)+h^T(k)y(k)\Big)\\=P(k)\Big((P^-(k)-h^T(k)h(k))\hat x(k-1)+h^T(k)y(k)\Big)\\=\hat x(k-1)+P(k)h^T(k)\Big(y(k)-h(k)\hat x(k-1)\Big) x^(k)=P(k)HT(k)Y(k)=P(k)(P(k1)x^(k1)+hT(k)y(k))=P(k)((P(k)hT(k)h(k))x^(k1)+hT(k)y(k))=x^(k1)+P(k)hT(k)(y(k)h(k)x^(k1))
    令: K ( k ) = P ( k − 1 ) h T ( k ) 1 + h ( k ) P ( k − 1 ) h T ( k ) = P ( k ) h T ( k ) K(k)=\frac{P(k-1)h^T(k)}{1+h(k)P(k-1)h^T(k)}=P(k)h^T(k) K(k)=1+h(k)P(k1)hT(k)P(k1)hT(k)=P(k)hT(k)
    则:
    P ( k ) = ( I − K ( k ) h ( k ) ) P ( k − 1 ) x ^ ( k ) = x ^ ( k − 1 ) + K ( k ) ( y ( k ) − h ( k ) x ^ ( k − 1 ) ) P(k)=(I-K(k)h(k))P(k-1)\\\hat x(k)=\hat x(k-1)+K(k)(y(k)-h(k)\hat x(k-1)) P(k)=(IK(k)h(k))P(k1)x^(k)=x^(k1)+K(k)(y(k)h(k)x^(k1))
    如何确定 P ( 0 ) P(0) P(0) 呢?
    显然: P ( k ) = ( ∑ i = 1 k h T ( i ) h ( i ) ) − ≈ ( ∑ i = 1 k h T ( i ) h ( i ) + 1 ∞ I ) − P(k)=\Big(\sum_{i=1}^kh^T(i)h(i)\Big)^- ≈ \Big(\sum_{i=1}^kh^T(i)h(i)+\frac{1}{∞ }I\Big)^- P(k)=(i=1khT(i)h(i))(i=1khT(i)h(i)+1I)
    P ( 0 ) = ∞ I P(0) = ∞I P(0)=I即可。那么此时 y ( 0 ) = h T ( 0 ) = x ^ ( 0 ) = O n × 1 y(0)=h^T(0)=\hat x(0) = O_{n\times 1} y(0)=hT(0)=x^(0)=On×1
    这样定义对代价函数影响忽略不计。

    带权RLS 递推计算

    当考虑带权递推时 x ^ ( k ) = ( H T ( k ) R ( k ) H ( k ) ) − H T ( k ) R ( k ) Y ( k ) \hat x(k)= \Big(H^T(k)R(k)H(k)\Big)^-H^T(k)R(k)Y(k) x^(k)=(HT(k)R(k)H(k))HT(k)R(k)Y(k)
    一般 R L S RLS RLS递推的 P ( k ) P(k) P(k)重新定义为: P ( k ) = ( H T ( k ) R ( k ) H ( k ) ) − = ( H T ( k − 1 ) R ( k − 1 ) H ( k − 1 ) + h T ( k ) r ( k ) h ( k ) ) − P(k) = (H^T(k)R(k)H(k))^-=\Big(H^T(k-1)R(k-1)H(k-1)+h^T(k)r(k)h(k)\Big)^- P(k)=(HT(k)R(k)H(k))=(HT(k1)R(k1)H(k1)+hT(k)r(k)h(k))
    令: γ ( k ) = 1 r ( k ) \gamma(k) = \frac{1}{r(k)} γ(k)=r(k)1
    矩阵反演有:
    P ( k ) = P ( k − 1 ) − P ( k − 1 ) h T ( k ) h ( k ) P ( k − 1 ) γ ( k ) − h ( k ) P ( k − 1 ) h T ( k ) P(k)=P(k-1)-\frac{P(k-1)h^T(k)h(k)P(k-1)}{\gamma(k)-h(k)P(k-1)h^T(k)} P(k)=P(k1)γ(k)h(k)P(k1)hT(k)P(k1)hT(k)h(k)P(k1)
    由: x ^ ( k ) = P ( k ) H T ( k ) R ( k ) Y ( k ) \hat x(k)=P(k)H^T(k)R(k)Y(k) x^(k)=P(k)HT(k)R(k)Y(k)
    有: P − ( k ) x ^ ( k ) = H T ( k ) R ( k ) Y ( k ) P^-(k)\hat x(k)=H^T(k)R(k)Y(k) P(k)x^(k)=HT(k)R(k)Y(k)
    由: P ( k ) = ( P − ( k − 1 ) + h T ( k ) r ( k ) h ( k ) ) P(k)=\Big(P^-(k-1)+h^T(k)r(k)h(k)\Big) P(k)=(P(k1)+hT(k)r(k)h(k))
    有: P − ( k − 1 ) = P − ( k ) − h T ( k ) r ( k ) h ( k ) P^-(k-1)=P^-(k)-h^T(k)r(k)h(k) P(k1)=P(k)hT(k)r(k)h(k)
    则: x ^ ( k ) = P ( k ) H T ( k ) R ( k ) Y ( k ) = P ( k ) ( P − ( k − 1 ) x ^ ( k − 1 ) + h T ( k ) r ( k ) y ( k ) ) = P ( k ) ( ( P − ( k ) − h T ( k ) r ( k ) h ( k ) ) x ^ ( k − 1 ) + h T ( k ) r ( k ) y ( k ) ) = x ^ ( k − 1 ) + P ( k ) h T ( k ) r ( k ) ( y ( k ) − h ( k ) x ( k ) ) \hat x(k)=P(k)H^T(k)R(k)Y(k)\\=P(k)\Big(P^-(k-1)\hat x(k-1)+h^T(k)r(k)y(k)\Big)\\=P(k)\Big((P^-(k)-h^T(k)r(k)h(k))\hat x(k-1)+h^T(k)r(k)y(k)\Big)\\ =\hat x(k-1)+P(k)h^T(k)r(k)\Big(y(k)-h(k)x(k)\Big) x^(k)=P(k)HT(k)R(k)Y(k)=P(k)(P(k1)x^(k1)+hT(k)r(k)y(k))=P(k)((P(k)hT(k)r(k)h(k))x^(k1)+hT(k)r(k)y(k))=x^(k1)+P(k)hT(k)r(k)(y(k)h(k)x(k))
    K ( k ) = P ( k ) h T ( k ) r ( k ) = P ( k − 1 ) h T ( k ) γ ( k ) − h ( k ) P ( k − 1 ) h T ( k ) K(k)=P(k)h^T(k)r(k)=\frac{P(k-1)h^T(k)}{\gamma(k)-h(k)P(k-1)h^T(k)} K(k)=P(k)hT(k)r(k)=γ(k)h(k)P(k1)hT(k)P(k1)hT(k)
    则: P ( k ) = ( I − K ( k ) h ( k ) ) P ( k − 1 ) P(k)=\Big(I-{K(k)h(k)}\Big)P(k-1) P(k)=(IK(k)h(k))P(k1)
    x ^ ( k ) = x ^ ( k − 1 ) + K ( k ) ( y ( k ) − h ( k ) x ^ ( k − 1 ) ) \hat x(k) = \hat x(k-1)+K(k)(y(k)-h(k)\hat x(k-1)) x^(k)=x^(k1)+K(k)(y(k)h(k)x^(k1))
    同上,当没有很好的初始设置数值时,可以讲, P ( 0 ) = ∞ I P(0)=∞I P(0)=I x ^ ( 0 ) = h T ( 0 ) = O n × 1 \hat x(0) = h^T(0)=O_{n\times 1} x^(0)=hT(0)=On×1

    带有遗忘的 R L S RLS RLS

    对于带有遗忘版本的 R L S RLS RLS算法,将对近期数据更为敏感。
    设置遗忘因子 λ ∈ ( 0 , 1 ] \lambda \in (0,1] λ(0,1]
    λ = 1 \lambda = 1 λ=1时,不会遗忘。
    每次更新,历史数据的权重会被整体缩小 λ \lambda λ
    那么此时:
    J ( k ) = λ J ( k − 1 ) + ( y ( k ) − h ( k ) x ^ ( k ) ) 2 J(k)=\lambda J(k-1) +(y(k)-h(k)\hat x(k))^2 J(k)=λJ(k1)+(y(k)h(k)x^(k))2
    则:
    J ( k ) = λ ∣ ∣ H ( k − 1 ) x ^ ( k ) − Y ( k − 1 ) ∣ ∣ 2 + ( y ( k ) − h ( k ) x ^ ( k ) ) 2 J(k)=\lambda ||H(k-1) \hat x(k)-Y(k-1)||^2+(y(k)-h(k)\hat x(k))^2 J(k)=λH(k1)x^(k)Y(k1)2+(y(k)h(k)x^(k))2
    从新定义
    H ( k ) = [ β H ( k − 1 ) ,   h ( k ) ] H(k)=[\beta H(k-1),\ h(k)] H(k)=[βH(k1), h(k)]
    Y ( k ) = [ β Y ( k − 1 ) ,   y ( k ) ] Y(k)=[\beta Y(k-1),\ y(k)] Y(k)=[βY(k1), y(k)]
    其中, β 2 = λ \beta^2= \lambda β2=λ
    P ( k ) = ( H T ( k ) H ( k ) ) − = ( H T ( k − 1 ) H ( k − 1 ) λ + h T ( k ) h ( k ) ) − P(k)=(H^T(k)H(k))^-\\=(H^T(k-1)H(k-1)\lambda+h^T(k)h(k))^- P(k)=(HT(k)H(k))=(HT(k1)H(k1)λ+hT(k)h(k))
    = ( P − ( k − 1 ) λ + h T ( k ) h ( k ) ) − =(P^-(k-1)\lambda+h^T(k)h(k))^- =(P(k1)λ+hT(k)h(k))
    = P ( k − 1 ) 1 λ − 1 λ 2 P ( k − 1 ) h T ( k ) h ( k ) P ( k − 1 ) 1 + 1 λ h ( k ) P ( k − 1 ) h T ( k ) =P(k-1)\frac{1}{\lambda}-\frac{\frac{1}{\lambda^2}P(k-1)h^T(k)h(k)P(k-1)}{1+\frac{1}{\lambda}h(k)P(k-1)h^T(k)} =P(k1)λ11+λ1h(k)P(k1)hT(k)λ21P(k1)hT(k)h(k)P(k1)
    = ( I − P ( k − 1 ) h T ( k ) h ( k ) λ + h ( k ) P ( k − 1 ) h T ( k ) ) P ( k − 1 ) 1 λ =\Big(I-\frac{P(k-1)h^T(k)h(k)}{\lambda+h(k)P(k-1)h^T(k)}\Big)P(k-1)\frac{1}{\lambda} =(Iλ+h(k)P(k1)hT(k)P(k1)hT(k)h(k))P(k1)λ1
    另一边:
    x ^ ( k ) = ( H T ( k ) H ( k ) ) − H T ( k ) Y ( k ) = P ( k ) H T ( k ) Y ( k ) \hat x(k)=(H^T(k)H(k))^-H^T(k)Y(k)\\ =P(k)H^T(k)Y(k) x^(k)=(HT(k)H(k))HT(k)Y(k)=P(k)HT(k)Y(k)
    x ^ ( k ) = P ( k ) H T ( k ) Y ( k ) \hat x(k)=P(k)H^T(k)Y(k) x^(k)=P(k)HT(k)Y(k)
    可得: P − ( k ) x ^ ( k ) = H T ( k ) Y ( k ) P^-(k)\hat x(k)=H^T(k)Y(k) P(k)x^(k)=HT(k)Y(k)
    同时:
    P ( k ) = ( P − ( k − 1 ) λ + h T ( k ) h ( k ) ) − P(k)=(P^-(k-1)\lambda +h^T(k)h(k))^- P(k)=(P(k1)λ+hT(k)h(k))
    P − ( k − 1 ) = 1 λ ( P − ( k ) − h T ( k ) h ( k ) ) P^-(k-1)=\frac{1}{\lambda}(P^-(k)-h^T(k)h(k)) P(k1)=λ1(P(k)hT(k)h(k))
    x ^ ( k ) = P ( k ) ( H T ( k − 1 ) Y ( k − 1 ) λ + h T ( k ) y ( k ) ) = P ( k ) ( P − ( k − 1 ) x ^ ( k − 1 ) λ + h T ( k ) y ( k ) ) = P ( k ) ( ( P − ( k ) − h T ( k ) h ( k ) ) x ^ ( k − 1 ) + h T ( k ) y ( k ) ) = x ^ ( k − 1 ) + P ( k ) h T ( k ) ( y ( k ) − h ( k ) x ^ ( k − 1 ) ) \hat x(k)=P(k)(H^T(k-1)Y(k-1)\lambda+h^T(k)y(k))\\ =P(k)(P^-(k-1)\hat x(k-1)\lambda +h^T(k)y(k))\\ =P(k)\Big(\big(P^-(k)-h^T(k)h(k)\big)\hat x(k-1)+h^T(k)y(k)\Big)\\ =\hat x(k-1)+P(k)h^T(k)\Big(y(k)-h(k)\hat x(k-1)\Big) x^(k)=P(k)(HT(k1)Y(k1)λ+hT(k)y(k))=P(k)(P(k1)x^(k1)λ+hT(k)y(k))=P(k)((P(k)hT(k)h(k))x^(k1)+hT(k)y(k))=x^(k1)+P(k)hT(k)(y(k)h(k)x^(k1))
    令: K ( k ) = P ( k ) h T ( k ) K(k)=P(k)h^T(k) K(k)=P(k)hT(k)
    那么:
    x ^ ( k ) = x ^ ( k − 1 ) + K ( k ) ( y ( k ) − h ( k ) x ^ ( k − 1 ) ) \hat x(k)=\hat x(k-1)+K(k)\Big(y(k)-h(k)\hat x(k-1)\Big) x^(k)=x^(k1)+K(k)(y(k)h(k)x^(k1))
    K ( k ) = P ( k ) h ( k ) = ( I − P ( k − 1 ) h T ( k ) h ( k ) λ + h ( k ) P ( k − 1 ) h T ( k ) ) P ( k − 1 ) 1 λ h T ( k ) = P ( k − 1 ) h T ( k ) λ + h ( k ) P ( k − 1 ) h T ( k ) K(k)=P(k)h(k)=\Big(I-\frac{P(k-1)h^T(k)h(k)}{\lambda+h(k)P(k-1)h^T(k)}\Big)P(k-1)\frac{1}{\lambda}h^T(k)\\ =\frac{P(k-1)h^T(k)}{\lambda +h(k)P(k-1)h^T(k)} K(k)=P(k)h(k)=(Iλ+h(k)P(k1)hT(k)P(k1)hT(k)h(k))P(k1)λ1hT(k)=λ+h(k)P(k1)hT(k)P(k1)hT(k)
    P ( k ) = ( I − K ( k ) h ( k ) ) P ( k − 1 ) 1 λ P(k)=(I-K(k)h(k))P(k-1)\frac{1}{\lambda} P(k)=(IK(k)h(k))P(k1)λ1

    总结

    RLS 算法一个比较精巧的在线滤波算法。 还可以增加遗忘因子,让其对近期近期数据更为敏感,比如笔记平滑种可能需要这种形式的滤波器。

    展开全文
  • LMS与RLS算法

    2019-04-03 20:26:28
    LMS与RLS算法,现代数字信号处理中极其重要的算法,本PPT简介了LMS与RLS算法的由来与推导,最小均方(LMS)自适应滤波器、递推最小二乘(RLS)滤波器格型滤波器都是我们自适应滤波器中较为强大的滤波器,也提及了他...
  • 自适应信号处理的理论和技术已经成为人们常用滤波和去噪技术。文中讲述了自适应滤波的原理以及LMS算法和RLS算法两种基本自适应算法的原理及步骤。并用MATLAB分别对两种算法进行了自适应滤波仿真和实现。
  • 时域LMS与RLS算法自适应滤波算法-RLS_finish2.m 本帖最后由 dingkillerwhale 于 2013-5-20 09:56 编辑 数据由于太大无法上传 以下是LMS算法以及RLS算法,其中RLS针对恒模信号与非恒模信号进行区分
  • 时域LMS与RLS算法自适应滤波算法-RLS_finish1.m 本帖最后由 dingkillerwhale 于 2013-5-20 09:56 编辑 数据由于太大无法上传 以下是LMS算法以及RLS算法,其中RLS针对恒模信号与非恒模信号进行区分
  • 快速RLS算法

    千次阅读 2019-08-15 18:15:37
    RLS算法需要的运算量比较大,在实际应用中会受到一些限制,因此出现了许多快速算法。本篇博文介绍了阶递归RLS和快速横向RLS。

    版权声明:本文为作者原创文章,如转载请注明出处。

    前言

    RLS算法需要的运算量比较大,在实际应用中会受到一些限制,因此出现了许多快速算法。比较常见的两种算法为阶递归自适应滤波器(Lattice-Based RLS)和快速横向滤波器(Fast Transversal RLS)。所有的快速算法都探索了信息数据的一些结构特性,以实现较低的计算复杂度。快速算法的模式是相似的,因为前向和后向预测滤波器是快速算法的基本组成部分。

    1 阶递归自适应滤波器

    对于一个预测过程,前向后验预测误差为
    ε f ( k , i + 1 ) = x ( k ) − w f T ( k , i + 1 ) x ( k − 1 , i + 1 ) \varepsilon_f(k,i+1)=x(k)-\bm{w}^T_f(k,i+1)\bm{x}(k-1,i+1) εf(k,i+1)=x(k)wfT(k,i+1)x(k1,i+1)式中 w f T ( k , i + 1 ) = [ w f 0 ( k ) , w f 1 ( k ) , w f 2 ( k ) , . . . w f i ( k ) ] T \bm{w}^T_f(k,i+1)=[w_{f0}(k),w_{f1}(k),w_{f2}(k),...w_{fi}(k)]^T wfT(k,i+1)=[wf0(k),wf1(k),wf2(k),...wfi(k)]T x ( k − 1 , i + 1 ) = [ x ( k − 1 ) , x ( k − 2 ) , . . . x ( k − i − 1 ) ] T \bm{x}(k-1,i+1)=[x(k-1),x(k-2),...x(k-i-1)]^T x(k1,i+1)=[x(k1),x(k2),...x(ki1)]T
    定义加权误差向量
    ε f ( k , i + 1 ) = x ^ ( k ) − X T ( k − 1 , i + 1 ) w f ( k , i + 1 ) \bm{\varepsilon}_f(k,i+1)=\widehat{\bm{x}}(k)-\bm{X}^T(k-1,i+1)\bm{w}_f(k,i+1) εf(k,i+1)=x (k)XT(k1,i+1)wf(k,i+1)
    式中 x ^ ( k ) = [ x ( k ) , λ 1 / 2 x ( k − 1 ) , , λ x ( k − 2 ) , . . . , λ k / 2 x ( 0 ) ] \widehat{\bm{x}}(k)=[x(k),\lambda^{1/2}x(k-1), ,\lambda x(k-2),...,\lambda^{k/2} x(0)] x (k)=[x(k),λ1/2x(k1),,λx(k2),...,λk/2x(0)]
    写为矩阵形式有
    ε f ( k , i + 1 ) = X T ( k , i + 2 ) [ 1 − w f ( k , i + 1 ) ] \bm{\varepsilon}_f(k,i+1)=\bm{X}^T(k,i+2)\begin{bmatrix} 1 \\ -\bm{w}_f(k,i+1) \\ \end{bmatrix} εf(k,i+1)=XT(k,i+2)[1wf(k,i+1)]
    均方误差为
    ξ f d ( k , i + 1 ) = ε f T ( k , i + 1 ) ε f ( k , i + 1 ) \xi^d_f(k,i+1)=\bm{\varepsilon}_f^T(k,i+1)\bm{\varepsilon}_f(k,i+1) ξfd(k,i+1)=εfT(k,i+1)εf(k,i+1)
    = ∑ l = 0 k λ k − l [ x ( l ) − x T ( l − 1 , i + 1 ) w f ( k , i + 1 ) ] 2 =\sum_{l=0}^k \lambda^{k-l}[x(l)-\bm{x}^T(l-1,i+1)\bm{w}_f(k,i+1)]^2 =l=0kλkl[x(l)xT(l1,i+1)wf(k,i+1)]2
    令导数为0可以得到
    w f ( k , i + 1 ) = [ ∑ l = 0 k λ k − l x ( l − 1 , i + 1 ) x T ( l − 1 , i + 1 ) ] − 1 ∑ l = 0 k λ k − l x ( l − 1 , i + 1 ) x ( l ) \bm{w}_f(k,i+1)=\left[ \sum_{l=0}^k \lambda^{k-l}\bm{x}(l-1,i+1)\bm{x}^T(l-1,i+1) \right]^{-1}\sum_{l=0}^k \lambda^{k-l} \bm{x}(l-1,i+1)x(l) wf(k,i+1)=[l=0kλklx(l1,i+1)xT(l1,i+1)]1l=0kλklx(l1,i+1)x(l)
    = [ X ( k − 1 , i + 1 ) X T ( k − 1 , i + 1 ) ] − 1 X ( k − 1 , i + 1 ) x ^ ( k ) =[\bm{X}(k-1,i+1)\bm{X}^T(k-1,i+1)]^{-1}\bm{X}(k-1,i+1)\widehat{\bm{x}}(k) =[X(k1,i+1)XT(k1,i+1)]1X(k1,i+1)x (k)
    = R D f − 1 ( k − 1 , i + 1 ) p D f ( k , i + 1 ) =\bm{R}^{-1}_{Df}(k-1,i+1)\bm{p}_{Df}(k,i+1) =RDf1(k1,i+1)pDf(k,i+1)
    将最佳权系数代入可以得到最小均方误差
    ξ f m i n d ( k , i + 1 ) = ∑ l = 0 k λ k − l x ( l ) [ x ( l ) − x T ( l − 1 , i + 1 ) w f ( k , i + 1 ) ] \xi^d_{fmin}(k,i+1)=\sum_{l=0}^k\lambda^{k-l}x(l)[x(l)-\bm{x}^T(l-1,i+1) \bm{w}_f(k,i+1)] ξfmind(k,i+1)=l=0kλklx(l)[x(l)xT(l1,i+1)wf(k,i+1)]
    = ∑ l = 0 k λ k − l x 2 ( l ) − p D f T ( k , i + 1 ) w f ( k , i + 1 ) =\sum_{l=0}^k\lambda^{k-l}x^2(l)-\bm{p}_{Df}^T(k,i+1)\bm{w}_f(k,i+1) =l=0kλklx2(l)pDfT(k,i+1)wf(k,i+1)
    = σ f 2 ( k ) − w f ( k , i + 1 ) T p D f ( k , i + 1 ) =\sigma^2_f(k)-\bm{w}_f(k,i+1)^T\bm{p}_{Df}(k,i+1) =σf2(k)wf(k,i+1)TpDf(k,i+1)
    将上面两式组合可以得到
    R D ( k , i + 2 ) = [ σ f 2 ( k ) p D f T ( k , i + 1 ) p D f T ( k , i + 1 ) R D f ( k − 1 , i + 1 ) ] \bm{R}_D(k,i+2)=\begin{bmatrix} \sigma^2_f(k)&\bm{p}^T_{Df}(k,i+1) \\ \bm{p}^T_{Df}(k,i+1) & \bm{R}_{Df}(k-1,i+1) \end{bmatrix} RD(k,i+2)=[σf2(k)pDfT(k,i+1pDfT(k,i+1RDf(k1,i+1)]

    R D ( k , i + 2 ) [ 1 − w f ( k , i + 1 ) ] = [ ξ f m i n d ( k , i + 1 ) 0 ] \bm{R}_D(k,i+2)\begin{bmatrix} 1 \\ -\bm{w}_{f}(k,i+1) \end{bmatrix}=\begin{bmatrix} \xi^d_{fmin}(k,i+1) \\ 0 \end{bmatrix} RD(k,i+2)[1wf(k,i+1]=[ξfmind(k,i+1)0]
    同理可得,后向预测模型为
    R D ( k , i + 2 ) [ − w b ( k , i + 1 ) 1 ] = [ 0 ξ b m i n d ( k , i + 1 ) ] \bm{R}_D(k,i+2)\begin{bmatrix} -\bm{w}_{b}(k,i+1) \\ 1 \end{bmatrix}=\begin{bmatrix} 0 \\ \xi^d_{b_{min}}(k,i+1) \end{bmatrix} RD(k,i+2)[wb(k,i+11]=[0ξbmind(k,i+1)]
    引入新变量 δ \delta δ
    R D ( k , i + 2 ) [ 1 − w f ( k , i ) 0 ] = [ ξ f m i n d ( k , i ) 0 δ f ( k , i ) ] \bm{R}_D(k,i+2)\begin{bmatrix} 1\\ -\bm{w}_{f}(k,i) \\ 0 \end{bmatrix}=\begin{bmatrix} \xi^d_{fmin}(k,i)\\ 0\\ \delta_f(k,i)\\ \end{bmatrix} RD(k,i+2)1wf(k,i0=ξfmind(k,i)0δf(k,i)

    式中 δ f ( k , i ) = ∑ l = 0 k λ k − l ε f ( l , i ) x ( l − i − 1 ) \delta_f(k,i)=\sum_{l=0}^k \lambda^{k-l}\varepsilon_f(l,i)x(l-i-1) δf(k,i)=l=0kλklεf(l,i)x(li1)
    同样的有
    R D ( k , i + 2 ) [ 0 − w b ( k − 1 , i ) 1 ] = [ δ b ( k , i ) 0 ξ b m i n d ( k , i ) ] \bm{R}_D(k,i+2)\begin{bmatrix} 0\\ -\bm{w}_{b}(k-1,i) \\ 1 \end{bmatrix}=\begin{bmatrix} \delta_b(k,i)\\ 0\\ \xi^d_{bmin}(k,i)\\ \end{bmatrix} RD(k,i+2)0wb(k1,i1=δb(k,i)0ξbmind(k,i)
    式中 δ b ( k , i ) = ∑ l = 0 k λ k − l ε b ( l − 1 , i ) x ( l ) \delta_b(k,i)=\sum_{l=0}^k \lambda^{k-l}\varepsilon_b(l-1,i)x(l) δb(k,i)=l=0kλklεb(l1,i)x(l)
    后验误差迭代为
    ξ b m i n d ( k , i + 1 ) = ξ b m i n d ( k − 1 , i ) − δ 2 ( k , i ) ξ f m i n d ( k , i ) \xi_{b_{min}}^d(k,i+1)=\xi_{b_{min}}^d(k-1,i)-\frac{\delta^2(k,i)}{\xi_{f_{min}}^d(k,i)} ξbmind(k,i+1)=ξbmind(k1,i)ξfmind(k,i)δ2(k,i)
    后验估计权系数迭代为
    w b ( k , i + 1 ) = [ 0 w b ( k − 1 , i ) ] − δ ( k , i ) ξ f m i n d ( k , i ) [ − 1 w f ( k , i ) ] \bm{w}_b(k,i+1)=\begin{bmatrix}0\\\bm{w}_{b}(k-1,i) \end{bmatrix}-\frac{\delta(k,i)}{\xi_{f_{min}}^d(k,i)}\begin{bmatrix}-1\\\bm{w}_{f}(k,i) \end{bmatrix} wb(k,i+1)=[0wb(k1,i)]ξfmind(k,i)δ(k,i)[1wf(k,i)]
    同样的,前验误差迭代为
    ξ f m i n d ( k , i + 1 ) = ξ f m i n d ( k , i ) − δ 2 ( k , i ) ξ b m i n d ( k − 1 , i ) \xi_{f_{min}}^d(k,i+1)=\xi_{f_{min}}^d(k,i)-\frac{\delta^2(k,i)}{\xi_{b_{min}}^d(k-1,i)} ξfmind(k,i+1)=ξfmind(k,i)ξbmind(k1,i)δ2(k,i)
    前验估计权系数迭代为
    w f ( k , i + 1 ) = [ w f ( k , i ) 0 ] − δ ( k , i ) ξ b m i n d ( k − 1 , i ) [ w b ( k − 1 , i ) − 1 ] \bm{w}_f(k,i+1)=\begin{bmatrix} \bm{w}_{f}(k,i) \\ 0 \end{bmatrix}-\frac{\delta(k,i)}{\xi_{b_{min}}^d(k-1,i)}\begin{bmatrix}\bm{w}_{b}(k-1,i) \\-1\end{bmatrix} wf(k,i+1)=[wf(k,i)0]ξbmind(k1,i)δ(k,i)[wb(k1,i)1]
    后验前向误差为
    ε f ( k , i + 1 ) = ε f ( k , i ) − κ f ( k , i ) ε b ( k − 1 , i ) \varepsilon_f(k,i+1)=\varepsilon_f(k,i)-\kappa_f(k,i)\varepsilon_b(k-1,i) εf(k,i+1)=εf(k,i)κf(k,i)εb(k1,i)
    式中 κ f ( k , i ) \kappa_f(k,i) κf(k,i)称为前向反射系数
    κ f ( k , i ) = δ ( k , i ) ξ b m i n d ( k − 1 , i ) \kappa_f(k,i)=\frac{\delta(k,i)}{\xi_{b_{min}}^d(k-1,i)} κf(k,i)=ξbmind(k1,i)δ(k,i)
    后验后向预测误差为
    ε b ( k , i + 1 ) = ε b ( k − 1 , i ) − κ b ( k , i ) ε f ( k , i ) \varepsilon_b(k,i+1)=\varepsilon_b(k-1,i)-\kappa_b(k,i)\varepsilon_f(k,i) εb(k,i+1)=εb(k1,i)κb(k,i)εf(k,i)
    式中 κ b ( k , i ) \kappa_b(k,i) κb(k,i)称为后向反射系数
    κ b ( k , i ) = δ ( k , i ) ξ f m i n d ( k , i ) \kappa_b(k,i)=\frac{\delta(k,i)}{\xi_{f_{min}}^d(k,i)} κb(k,i)=ξfmind(k,i)δ(k,i)
    在初始阶段
    ε b ( k , 0 ) = ε f ( k , 0 ) = x ( k ) \varepsilon_b(k,0)=\varepsilon_f(k,0)=x(k) εb(k,0)=εf(k,0)=x(k)
    ξ f m i n d ( k , 0 ) = ξ b m i n d ( k , 0 ) = x 2 ( k ) + λ ξ f m i n d ( k − 1 , 0 ) \xi_{f_{min}}^d(k,0)=\xi_{b_{min}}^d(k,0)=x^2(k)+\lambda\xi_{f_{min}}^d(k-1,0) ξfmind(k,0)=ξbmind(k,0)=x2(k)+λξfmind(k1,0)
    δ \delta δ的迭代为
    δ ( k , i ) = λ δ ( k − 1 , i ) + ε b ( k − 1 , i ) ε f ( k , i ) γ ( k − 1 , i ) \delta(k,i)=\lambda\delta(k-1,i)+\frac{\varepsilon_b(k-1,i)\varepsilon_f(k,i)}{\gamma(k-1,i)} δ(k,i)=λδ(k1,i)+γ(k1,i)εb(k1,i)εf(k,i)
    式中
    γ ( k , i + 1 ) = γ ( k , i ) − ε b 2 ( k , i ) ξ b m i n d ( k , i ) \gamma(k,i+1)=\gamma(k,i)-\frac{\varepsilon^2_b(k,i)}{\xi^d_{b_{min}}(k,i)} γ(k,i+1)=γ(k,i)ξbmind(k,i)εb2(k,i)
    加权均方误差为
    ξ d ( k , i + 1 ) = ∑ l = 0 k λ k − l ε 2 ( l , i + 1 ) \xi^d(k,i+1)=\sum^k_{l=0}\lambda^{k-l}\varepsilon^2(l,i+1) ξd(k,i+1)=l=0kλklε2(l,i+1)
    对应的
    δ D ( k , i ) = λ δ D ( k − 1 , i ) + ε ( k , i ) ε b ( k , i + 1 ) γ ( k , i ) \delta_D(k,i)=\lambda\delta_D(k-1,i)+\frac{\varepsilon(k,i)\varepsilon_b(k,i+1)}{\gamma(k,i)} δD(k,i)=λδD(k1,i)+γ(k,i)ε(k,i)εb(k,i+1)
    前向乘子为
    v i ( k ) = δ D ( k , i ) ξ b m i n d ( k , i ) v_i(k)=\frac{\delta_D(k,i)}{\xi^d_{b_{min}}(k,i)} vi(k)=ξbmind(k,i)δD(k,i)
    后验输出误差为
    ε ( k , i + 1 ) = ε ( k , i ) − v i ( k ) ε b ( k , i ) \varepsilon(k,i+1)=\varepsilon(k,i)-v_i(k)\varepsilon_b(k,i) ε(k,i+1)=ε(k,i)vi(k)εb(k,i)
    由此获得阶递归RLS算法的全过程。

    2 快速横向滤波器

    后验前向预测误差为
    ε f ( k , N ) = x ( k ) − w f T ( k , N ) x ( k − 1 , N ) \varepsilon_f(k,N)=x(k)-\bm{w}^T_f(k,N)\bm{x}(k-1,N) εf(k,N)=x(k)wfT(k,N)x(k1,N)
    = x T ( k , N + 1 ) [ 1 − w f ( k , N ) ] =\bm{x}^T(k,N+1)\begin{bmatrix}1\\-\bm{w}_f(k,N)\\\end{bmatrix} =xT(k,N+1)[1wf(k,N)]
    后验与先验前向预测误差之间的关系为
    e f ( k , N ) = ε f ( k , N ) γ ( k − 1 , N ) e_f(k,N)=\frac{\varepsilon_f(k,N)}{\gamma(k-1,N)} ef(k,N)=γ(k1,N)εf(k,N)
    式中
    γ ( k , N + 1 ) = λ ξ f m i n d ( k − 1 , N ) ξ f m i n d ( k , N ) γ ( k − 1 , N ) \gamma(k,N+1)=\frac{\lambda\xi^d_{f{min}}(k-1,N)}{\xi^d_{f{min}}(k,N)}\gamma(k-1,N) γ(k,N+1)=ξfmind(k,N)λξfmind(k1,N)γ(k1,N)
    前向最小加权均方误差为
    ξ f m i n d ( k , N ) = λ ξ f m i n d ( k − 1 , N ) + e f ( k , N ) ε f ( k , N ) \xi^d_{f_{min}}(k,N)=\lambda\xi^d_{f_{min}}(k-1,N)+e_f(k,N)\varepsilon_f(k,N) ξfmind(k,N)=λξfmind(k1,N)+ef(k,N)εf(k,N)
    前向预测权系数迭代
    w f ( k , N ) = w f ( k − 1 , N ) + ϕ ( k − 1 , N ) e f ( k , N ) \bm{w}_f(k,N)=\bm{w}_f(k-1,N)+\phi(k-1,N)e_f(k,N) wf(k,N)=wf(k1,N)+ϕ(k1,N)ef(k,N)
    式中
    ϕ ( k , N + 1 ) = [ 0 ϕ ( k − 1 , N ) ] + 1 ξ f m i n d ( k , N ) [ 1 w f ( k , N ) ] ε f ( k , N ) \phi(k,N+1)=\begin{bmatrix}0\\\phi(k-1,N)\end{bmatrix}+\frac{1}{\xi^d_{f_{min}}(k,N)}\begin{bmatrix}1\\\bm{w}_f(k,N)\\\end{bmatrix}\varepsilon_f(k,N) ϕ(k,N+1)=[0ϕ(k1,N)]+ξfmind(k,N)1[1wf(k,N)]εf(k,N)
    利用 ϕ ^ ( k , N + 1 ) \widehat\phi(k,N+1) ϕ (k,N+1)代替 ϕ ( k , N + 1 ) \phi(k,N+1) ϕ(k,N+1)
    ϕ ^ ( k , N + 1 ) = [ 0 ϕ ^ ( k − 1 , N ) ] + 1 λ ξ f m i n d ( k − 1 , N ) [ 1 − w f ( k − 1 , N ) ] e f ( k , N ) \widehat\phi(k,N+1)=\begin{bmatrix}0\\\widehat\phi(k-1,N)\end{bmatrix}+\frac{1}{\lambda\xi^d_{f_{min}}(k-1,N)}\begin{bmatrix}1\\-\bm{w}_f(k-1,N)\\\end{bmatrix}e_f(k,N) ϕ (k,N+1)=[0ϕ (k1,N)]+λξfmind(k1,N)1[1wf(k1,N)]ef(k,N)
    则前向预测权系数迭代可以改为
    w f ( k , N ) = w f ( k − 1 , N ) + ϕ ^ ( k − 1 , N ) ε f ( k , N ) \bm{w}_f(k,N)=\bm{w}_f(k-1,N)+\widehat\phi(k-1,N)\varepsilon_f(k,N) wf(k,N)=wf(k1,N)+ϕ (k1,N)εf(k,N)
    后验与先验预测误差之间的关系为
    ε b ( k , N ) = e b ( k , N ) γ ( k , N ) \varepsilon_b(k,N)=e_b(k,N)\gamma(k,N) εb(k,N)=eb(k,N)γ(k,N)
    转换因子
    γ ( k , N + 1 ) γ ( k , N ) = λ ξ b m i n d ( k − 1 , N ) ξ b m i n d ( k , N ) \frac{\gamma(k,N+1)}{\gamma(k,N)}=\frac{\lambda\xi^d_{b_{min}}(k-1,N)}{\xi^d_{b_{min}}(k,N)} γ(k,N)γ(k,N+1)=ξbmind(k,N)λξbmind(k1,N)
    后向先验预测误差为
    e b ( k , N ) = λ ξ b m i n d ( k − 1 , N ) ϕ N + 1 ( k , N + 1 ) e_b(k,N)=\lambda\xi^d_{b_{min}}(k-1,N)\phi_{N+1}(k,N+1) eb(k,N)=λξbmind(k1,N)ϕN+1(k,N+1)
    FTF中转换因子的迭代可以表示为
    γ − 1 ( k , N ) = γ − 1 ( k , N + 1 ) − ϕ ^ N + 1 ( k , N + 1 ) e b ( k , N ) \gamma^{-1}(k,N)=\gamma^{-1}(k,N+1)-\widehat\phi_{N+1}(k,N+1)e_b(k,N) γ1(k,N)=γ1(k,N+1)ϕ N+1(k,N+1)eb(k,N)
    后向最小加权均方误差为
    ξ b m i n d ( k , N ) = λ ξ b m i n d ( k − 1 , N ) + e b ( k , N ) ε b ( k , N ) \xi^d_{b_{min}}(k,N)=\lambda\xi^d_{b_{min}}(k-1,N)+e_b(k,N)\varepsilon_b(k,N) ξbmind(k,N)=λξbmind(k1,N)+eb(k,N)εb(k,N)
    则后向预测权系数迭代为
    w b ( k , N ) = w b ( k − 1 , N ) + ϕ ^ ( k , N ) ε b ( k , N ) \bm{w}_b(k,N)=\bm{w}_b(k-1,N)+\widehat\phi(k,N)\varepsilon_b(k,N) wb(k,N)=wb(k1,N)+ϕ (k,N)εb(k,N)
    式中
    [ ϕ ^ ( k , N ) 0 ] = ϕ ^ ( k , N + 1 ) − ϕ N + 1 ( k , N + 1 ) [ − w b ( k − 1 , N ) 1 ] \begin{bmatrix} \widehat\phi(k,N)\\ 0\\ \end{bmatrix}=\widehat\phi(k,N+1)-\phi_{N+1}(k,N+1)\begin{bmatrix}-\bm{w}_b(k-1,N)\\1\\\end{bmatrix} [ϕ (k,N)0]=ϕ (k,N+1)ϕN+1(k,N+1)[wb(k1,N)1]
    滤波先验误差为
    e ( k , N ) = d ( k ) − w T ( k − 1 , N ) x ( k , N ) e(k,N)=d(k)-\bm{w}^T(k-1,N)\bm{x}(k,N) e(k,N)=d(k)wT(k1,N)x(k,N)
    则后验误差为
    ε ( k , N ) = e ( k , N ) γ ( k , N ) \varepsilon(k,N)=e(k,N)\gamma(k,N) ε(k,N)=e(k,N)γ(k,N)
    滤波权系数迭代为
    w ( k , N ) = w ( k − 1 , N ) + ϕ ^ ( k , N ) ε ( k , N ) \bm{w}(k,N)=\bm{w}(k-1,N)+\widehat\phi(k,N)\varepsilon(k,N) w(k,N)=w(k1,N)+ϕ (k,N)ε(k,N)
    由此获得快速横向RLS算法的全过程。

    展开全文
  • RLS递归最小二乘方自适应算法源程序-rls算法.rar RLS(递归最小二乘方自适应算法源程序)
  • RLS算法MAIN函数

    2013-01-15 23:22:38
    这是RLS算法的MATLAB的MAIN函数。通过与上一个文档的联合应用,可以得到RLS算法
  • 包含MVDR算法,RLS算法,奇异值分解MVDR算法用于功率谱估计的MATLAB仿真,里面附有详细的注释,有不理解的地方可以评论区回复,适合初学者参考
  • 文讲述了LMS算法和RLS算法的基本原理,基于两种算法的推导过程较为繁琐,不便理解,本文以两种算法的主要计算环节为突破口,选取合适的迭代公式,对两种算法进行公式推导,方便读者理解。同时,本文通过原理推导、...
  • RLS算法的matlab代码

    2017-04-17 16:39:13
    本例比较了在四种特征值扩散度不同的情况下RLS算法的学习曲线
  • 时域LMS与RLS算法自适应滤波算法-LMS_finish.m 本帖最后由 dingkillerwhale 于 2013-5-20 09:56 编辑 数据由于太大无法上传 以下是LMS算法以及RLS算法,其中RLS针对恒模信号与非恒模信号进行区分
  • RLS算法Matlab实现

    千次阅读 2020-04-20 22:20:00
    RLS算法Matlab实现RLS迭代Matlab实现 递归最小二乘算法matlab代码实现 RLS迭代 Matlab实现 // RLS clc;clear Delta=[1 0.1 0.01 0.001]; Lamda=[1 0.9 0.5 0.1]; cishu = 1; for cishu=1:4 dotnumber=10000; u=...

    RLS算法Matlab实现

    递归最小二乘算法matlab代码实现

    RLS迭代

    在这里插入图片描述

    Matlab实现

    // RLS
    clc;clear
    Delta=[1 0.1 0.01 0.001];
    Lamda=[1 0.9 0.5 0.1];
    cishu = 1;
    for cishu=1:4
        
        dotnumber=10000;
        u=wgn(dotnumber,1,0);
        b=1;a=[1 -0.9];
        u=filter(b,a,u);
        u(dotnumber+1:end)=[];
        noise=wgn(dotnumber,1,-60);
        h=[-0.1 0.2 0.7 0.4 -0.2 -0.1 0.12 -0.25].';
        wo=[h;zeros(4,1)];
        d=filter(h,1,u);
        d(dotnumber+1:end)=[];
        d=d+noise;
        len=12;
        mu = 0.001;
        w = zeros(len,1);
        w_eror = zeros(dotnumber,1);
        e = zeros(dotnumber,1);
        
        vector_u = zeros(len,1);
        vector_d = zeros(len,1);
    %	delta = 0.01;    %固定delta
        delta = Delta(cishu);   %四次迭代,不同delta
        P = 1/delta*eye(len);
        lamda = 1;      %固定delta
    %    lamda = Lamda(cishu);   %四次迭代,不同delta    
    %    Pi = zeros(1,len);
        
    for n=1:dotnumber
        vector_u = [u(n); vector_u(1:end-1)];
    %    vector_d = [d(n); vector_d(1:end-1)];
        Pi = P * vector_u;
        k = Pi/(lamda + vector_u'*Pi);
        kesai = d(n) - w'*vector_u;
        w = w + k*conj(kesai);
        P =  (1/lamda)*P - (1/lamda)*k*vector_u'*P;
        w_error(n)=norm(w-wo)^2;
    end
    
    hold on
    plot(10*log10(w_error))
    title('lamda = 1')
    %title('delta = 0.01')
    xlabel('采样点')
    ylabel('系数误差 (dB)');
    %legend('RLS','LMS','AP');
    legend('delta=1','delta=0.1','delta=0.01','delta=0.001');
    %legend('lamda=1','lamda=0.9','lamda=0.5','lamda=0.1');
    end
    
    
    展开全文
  • 加权正交约束的ICA修正RLS算法
  • 用matlab编写的一个用RLS算法对一维信号降噪代码
  • 介绍了自适应均衡器下的LMS和RLS算法的基本原理,并分析了2种算法中的忘却因子μ对LMS和RLS算法收敛性能的影响。通过仿真可知,在相同忘却因子下,RLS算法的收敛速度明显快于LMS算法,并且误差也比LMS算法小。
  • RLS算法及遗忘因子对RLS的影响,修改RLS的参数来查看遗忘因子对其的性能的影响
  • (1)阅读参考资料和文献,明晰算法的计算过程,理解RLS算法基本过程; (2)主麦克风录制的语音信号是RLSprimsp.wav,参考麦克风录制的参考噪声是RLSrefns.wav,用Matlab指令读取; (3)根据算法编写相应的Matlab...
  • LMS和RLS算法分析

    2011-05-11 21:06:31
    自适应 LMS 算法 RLS 算法 算法仿真分析
  • 基于MATLAB的LMS和RLS算法滤波完整程序,一种是LMS算法,一种是RLS算法,并在MATLAB环境下对其进行了编程.对自适应滤波器性能进行分析
  • 具有稀疏RLS算法的在线顺序回波状态网络
  • LMS与RLS算法学习笔记

    万次阅读 多人点赞 2018-10-21 20:13:22
    LMS与RLS算法学习笔记一、 研究目的1.1最陡下降法理论1.2$LMS算法$1.3$RLS算法$1.4研究目标二、代码解析三、结果 一、 研究目的 1.1最陡下降法理论 LMS算法总是与最陡下降法联合使用,顾名思义,最陡下降法就是...


    实现代码点击 这里下载

    一、 研究目的

    1.1最陡下降法理论

    LMS算法总是与最陡下降法联合使用,顾名思义,最陡下降法就是沿着性能曲面最陡放方向向下(曲面负梯度方向)搜索曲面的最低点。迭代过程首先从曲面上某个初始点(对应与初始权矢量w(0) )出发,沿着该点负梯度方向搜索至第一点(对应与初始权矢量w(1) ,且w(1)等于初始值w(0)加上一个正比于负梯度的增量 )。以此类推,直到搜索到最低点 w* 。所以最陡下降法迭代计算权矢量的公式为:
    w ( n + 1 ) = w ( n ) + μ ( − ∇ ( n ) ) \boldsymbol{w(n+1)}=\boldsymbol{w(n)}+\mu(-\boldsymbol{\nabla}(n)) w(n+1)=w(n)+μ((n))
    式中, μ \mu μ是控制搜索步长的参数称为自适应增益常数。那么如何选取合适的 μ \mu μ值?
    由于推到十分复杂,这里简单给出结论:
    0 &lt; μ &lt; λ m a x − 1 0&lt;\mu&lt;\lambda^{-1}_{max} 0<μ<λmax1
    式中, λ m a x \lambda_{max} λmax R ( x ( n ) 自 相 关 矩 阵 ) R(x(n)自相关矩阵) Rx(n)的最大特征值。当然有时为了免去计算 R R R的特征值的麻烦,因为 R R R是正定的,所以有:
    t r [ R ] = ∑ k = 0 L λ k &gt; λ m a x t_r[R]=\sum_{k=0}^L\lambda_{k}&gt;\lambda_{max} tr[R]=k=0Lλk>λmax
    这里, t r [ R ] t_r[R] tr[R] R R R的迹,它可以用输入信号的取样值进行估计,即:
    t r [ R ] = ∑ k = 0 L E [ x i 2 ( n ) ] t_r[R]=\sum_{k=0}^LE[x^2_i(n)] tr[R]=k=0LE[xi2(n)]
    所以进一步有:
    0 &lt; μ &lt; t r − 1 [ R ] 0&lt;\mu&lt;t_r^{-1}[R] 0<μ<tr1[R]

    1.2 L M S LMS LMS算法

    由于最陡下降法每次迭代都需要知道性能曲面上某点的梯度值,而实际上梯度值只能根据观察数据进行估计。而 L M S LMS LMS实质上是用平方误差代替均方误差,即:
    ∇ ( n ) ≈ ∇ ^ ( n ) ≜ ∂ e 2 ( n ) ∂ w = [ ∂ e 2 ( n ) ∂ w 0   ∂ e 2 ( n ) ∂ w 1 ⋯ ∂ e 2 ( n ) ∂ w L ] T \boldsymbol{\nabla}(n)\approx\hat\boldsymbol{\nabla}(n)\triangleq\frac{\partial e^2(n)}{\partial w}=[\frac{\partial e^2(n)}{\partial w_0}\ \frac{\partial e^2(n)}{\partial w_1}\cdots\frac{\partial e^2(n)}{\partial w_L}]^T (n)^(n)we2(n)=[w0e2(n) w1e2(n)wLe2(n)]T
    可以得到:
    ∇ ^ ( n ) = 2 e ( n ) ∂ e ( n ) ∂ w = 2 e ( n ) x ( n ) \hat\boldsymbol{\nabla}(n)=2\boldsymbol{e(n)}\frac{\partial \boldsymbol{e(n)}}{\partial w}=2\boldsymbol{e(n)x(n)} ^(n)=2e(n)we(n)=2e(n)x(n)
    得到 L M S LMS LMS算法的基本关系式:
    w ( n + 1 ) = w ( n ) − μ ( ∇ ( n ) ) = w ( n ) + 2 μ e ( n ) x ( n ) \boldsymbol{w(n+1)}=\boldsymbol{w(n)}-\mu(\boldsymbol{\nabla}(n))=\boldsymbol{w(n)}+2\mu \boldsymbol{e(n)x(n)} w(n+1)=w(n)μ((n))=w(n)+2μe(n)x(n)
    我们可以证明 ∇ ^ ( n ) \hat\boldsymbol{\nabla}(n) ^(n) ∇ ( n ) \boldsymbol{\nabla}(n) (n)的无偏估计,换句话说,如果在每次迭代调整权矢量前能够进行多次观测,获得多个 x ( n ) x(n) x(n),然后一句梯度的统计平均值 E [ ∇ ^ ( n ) ] E[\hat\boldsymbol{\nabla}(n)] E[^(n)]来调整权矢量,则迭代结果必定能得到理想的最佳权矢量。

    1.3 R L S RLS RLS算法

    R L S RLS RLS算法的关键是用二乘方的时间平均准则取代最小均方准则,并按照时间进行迭代计算,换句话说,对从起始时刻到当前时刻所有误差的平方进行平均并使之最小化,即:
    ϵ ( n ) = ∑ k = 0 n e 2 ( k ) = m i n \epsilon(n)=\sum_{k=0}^{n}e^2(k)=min ϵ(n)=k=0ne2(k)=min
    其中:
    e ( k ) = d ( k ) − y ( k ) e(k)=d(k)-y(k) e(k)=d(k)y(k)
    对于,非平稳随机信号,为了更好的跟踪,引入一个指数加权因子对上式进行修正:
    ϵ ( n ) = ∑ k = 0 n λ n − k e 2 ( k ) = m i n \epsilon(n)=\sum_{k=0}^{n}\lambda^{n-k}e^2(k)=min ϵ(n)=k=0nλnke2(k)=min
    式中, λ \lambda λ称之为遗忘因子,它是小于1的正数,展开上式可以知道,新到的数据比就到的数据更重要。
    这里定义两个新的矩阵(被加权后的自相关和互相关),即:
    R ( n ) = ∑ k = 0 n λ n − k x ( k ) x T ( k ) \boldsymbol{R}(n)=\sum_{k=0}^{n}\lambda^{n-k}\boldsymbol{x(k)x^T(k)} R(n)=k=0nλnkx(k)xT(k)
    P ( n ) = ∑ k = 0 n λ n − k d ( k ) x ( k ) \boldsymbol{P}(n)=\sum_{k=0}^{n}\lambda^{n-k}\boldsymbol{d(k)x(k)} P(n)=k=0nλnkd(k)x(k)
    可以推出 R L S RLS RLS算法基本关系:
    w ( n ) = w ( n − 1 ) + k ( n ) e ( n ∣ n − 1 ) \boldsymbol{w(n)}=\boldsymbol{w(n-1)}+ \boldsymbol{k(n)e(n|n-1)} w(n)=w(n1)+