从维纳霍夫方程到自适应稳定快速RLS算法（一）
前言
大家好，我是小潘！本系列将试图从维纳霍夫方程的角度一步一步理解自适应稳定快速RLS算法（SFTRLS Algorithm）的由来。(ps:第一次用Markdown编辑器，第一次用KaTex，而且预览模式后才发现排版与预想的不对，改死人了！经历修改之后，若文章当中有错误之处，恳请大家批评指正，共同学习~~)
一、什么是维纳霍夫(Wiener-Hopf)？
①首先知道什么是维纳滤波器(Wiener Filter)？
图1 维纳滤波器
这是一个FIR最优线性滤波器模型,假设序列{ s ( n ) s(n) s(n)},{ w ( n ) w(n) w(n)},{ d ( n ) d(n) d(n)}是零均值广义平稳过程,在最小均方误差(MMSE)意义上的最优线性滤波器称为维纳滤波器。
从图1所示滤波器的估计误差为:
e
(
n
)
=
d
(
n
)
−
y
(
n
)
(1.1.1)
e(n) = d(n) - y(n) \tag{1.1.1}
e(n)=d(n)−y(n)(1.1.1)
其中
d
(
n
)
d(n)
d(n)为期望序列,滤波器输出
y
(
n
)
y(n)
y(n)是对
d
(
n
)
d(n)
d(n)的估计，且假定滤波器的长度为M，权系数为{
h
k
\mathop{h}_k
hk，0
⩽
\leqslant
⩽k
⩽
\leqslant
⩽M-1}，这样其输出
y
(
n
)
y(n)
y(n)就和输入序列
x
(
n
)
,
x
(
n
−
1
)
,
.
.
.
.
.
x
(
n
−
M
+
1
)
x(n),x(n-1),.....x(n-M+1)
x(n),x(n−1),.....x(n−M+1)有关,其关系可表示为:
y
(
n
)
=
∑
k
=
0
M
−
1
h
(
n
)
x
(
n
−
k
)
(1.1.2)
y(n)=\displaystyle\sum_{k=0}^{M-1}h(n)x(n-k) \tag{1.1.2}
y(n)=k=0∑M−1h(n)x(n−k)(1.1.2)
为了方便，亦可用矢量表示：
y
(
n
)
=
h
M
T
(
n
)
X
M
(
n
)
(1.1.3)
y(n)=\mathop{h}_M^{T}(n)\mathop{X}_M(n) \tag{1.1.3}
y(n)=hMT(n)XM(n)(1.1.3)
其中
h
M
T
\mathop{h}_M^{T}
hMT为第n时刻的1xM维的滤波器系数矢量，
X
M
(
n
)
\mathop{X}_M(n)
XM(n)为第n时刻的Mx1的数据输入矢量，它们分别为
h
M
(
n
)
=
[
h
0
(
n
)
h
1
(
n
)
.
.
.
h
M
−
1
(
n
)
]
(1.1.4)
\mathop{h}_M^{}(n)=[\mathop{h}_0(n) \quad \mathop{h}_1(n) \quad ... \quad \mathop{h}_{M-1}(n)] \tag{1.1.4}
hM(n)=[h0(n)h1(n)...hM−1(n)](1.1.4)
和
X
M
(
n
)
=
[
x
(
n
)
x
(
n
−
1
)
.
.
.
x
(
n
−
M
+
1
)
]
T
(1.1.5)
\mathop{X}_M^{}(n)=[x(n) \quad x(n-1) \quad ...\quad x(n-M+1)]\mathop{}_{}^{T} \tag{1.1.5}
XM(n)=[x(n)x(n−1)...x(n−M+1)]T(1.1.5)
于是估计误差定义为
e
(
n
)
=
d
(
n
)
−
h
M
T
(
n
)
X
M
(
n
)
(1.1.6)
e(n)=d(n)-\mathop{h}_M^{T}(n)\mathop{X}_M(n) \tag{1.1.6}
e(n)=d(n)−hMT(n)XM(n)(1.1.6)
②最小均方误差准则(MMSE)与维纳霍夫方程(Wiener-Hopf equations)
在自适应滤波应用中，存在两个性能较好的准则:最小二乘方准则(the least squares criterion)及其在统计意义上的对等公式—均方误差准则(MSE,mean-square-error)。最小二乘方和MSE准则均产生了一个二次性能指标，且其为滤波器权系数矢量的二次函数，有单个全局最小值。
定义估计误差
e
(
n
)
e(n)
e(n)的均方误差(也可表示为平均功率)为
J
(
h
M
)
=
E
[
∣
e
(
n
)
∣
2
]
=
E
[
e
(
n
)
e
(
n
)
]
(1.2.1)
J(\mathop{h}_M)=E[\vert{e(n)}\vert\mathop{}_{}^{2}]=E[e(n)e(n)] \tag{1.2.1}
J(hM)=E[∣e(n)∣2]=E[e(n)e(n)](1.2.1)
把式(1.1.6)代入式(1.2.1),得
J
(
h
M
)
=
E
{
[
d
(
n
)
−
h
M
T
(
n
)
X
M
(
n
)
]
[
d
(
n
)
−
h
M
T
(
n
)
X
M
(
n
)
]
}
=
E
[
d
(
n
)
2
]
−
2
h
M
T
(
n
)
E
[
d
(
n
)
X
M
(
n
)
]
+
h
M
T
(
n
)
E
[
X
M
(
n
)
X
M
T
(
n
)
]
h
M
(
n
)
(1.2.2)
\begin{aligned} J(\mathop{h}_M)&=E\{[d(n)-\mathop{h}_M^{T}(n)\mathop{X}_M(n)][d(n)-\mathop{h}_M^{T}(n)\mathop{X}_M(n)] \} \\ &=E[d(n)\mathop{}_{}^{2}]-2\mathop{h}_M^{T}(n)E[d(n)\mathop{X}_M(n)]+\mathop{h}_M^{T}(n)E[\mathop{X}_M(n)\mathop{X}_M^{T}(n)]\mathop{h}_M^{}(n) \tag{1.2.2} \end{aligned}
J(hM)=E{[d(n)−hMT(n)XM(n)][d(n)−hMT(n)XM(n)]}=E[d(n)2]−2hMT(n)E[d(n)XM(n)]+hMT(n)E[XM(n)XMT(n)]hM(n)(1.2.2)
因为滤波器权系数矢量 h M \mathop{h}_M^{} hM是一个独立的确定量，故放在数学期望 E E E[ ]之外。
观察上式(1.2.2)第一项，定义 σ d 2 = E [ d ( n ) 2 ] \mathop{\sigma}_{d}^{2}=E[d(n)\mathop{}_{}^{2}] σd2=E[d(n)2]，这是期望响应的平均功率，也是期望响应的方差。
再观察上式(1.2.2)第二项，定义其中 Mx1维 的互相关矢量的数学期望 p \textbf{p} p为
p
=
E
[
d
(
n
)
X
M
(
n
)
]
=
[
E
{
x
(
n
)
d
(
n
)
}
E
{
x
(
n
−
1
)
d
(
n
)
}
.
.
.
E
{
x
(
n
−
M
+
1
)
d
(
n
)
}
]
=
[
p
(
0
)
p
(
−
1
)
.
.
.
p
(
−
M
+
1
)
]
(1.2.3)
\textbf{p}=E[d(n)\mathop{X}_M(n)]=\begin{bmatrix} E\{ x(n)d(n)\}\\ E\{ x(n-1)d(n)\} \\ . \\ . \\ . \\ E\{ x(n-M+1)d(n)\} \end{bmatrix}=\begin{bmatrix} p(0) \\ p(-1) \\ . \\ . \\ . \\ p(-M+1) \end{bmatrix} \tag{1.2.3}
p=E[d(n)XM(n)]=
E{x(n)d(n)}E{x(n−1)d(n)}...E{x(n−M+1)d(n)}
=
p(0)p(−1)...p(−M+1)
(1.2.3)
其中，
p
(
−
m
)
p(-m)
p(−m)为输入数据
x
(
n
−
m
)
x(n-m)
x(n−m)与期望响应d(n)的互相关函数。
再再观察上式(1.2.2)的第三项，定义自相关矩阵 的数学期望
R
\textbf{R}
R为
R
=
E
[
X
M
(
n
)
X
M
T
(
n
)
]
=
[
r
(
0
)
r
(
1
)
.
.
.
r
(
M
−
1
)
r
(
−
1
)
r
(
0
)
.
.
.
r
(
M
−
2
)
.
.
.
.
.
.
.
.
.
.
.
.
r
(
−
M
+
1
)
.
.
.
.
.
.
r
(
0
)
]
(1.2.4)
\textbf{R}=E[\mathop{X}_M(n)\mathop{X}_M^{T}(n)]=\begin{bmatrix} r(0) & r(1) &...&r(M-1) \\ r(-1) & r(0) &...&r(M-2)\\ ... & ... &...&...\\ r(-M+1) & ... &...&r(0)\\ \end{bmatrix}\tag{1.2.4}
R=E[XM(n)XMT(n)]=
r(0)r(−1)...r(−M+1)r(1)r(0)..................r(M−1)r(M−2)...r(0)
(1.2.4)
其中，自相关函数
r
(
i
−
k
)
r(i-k)
r(i−k)的定义为
r
(
i
−
k
)
=
E
[
x
(
n
−
k
)
x
(
n
−
i
)
]
(1.2.5)
r(i-k)=E[x(n-k)x(n-i)] \tag{1.2.5}
r(i−k)=E[x(n−k)x(n−i)](1.2.5)
通过
σ
d
2
，
p
，
R
\mathop{\sigma}_{d}^{2}，\textbf{p}，\textbf{R}
σd2，p，R重新定义式(1.2.2)的均方误差
J
(
h
M
)
=
σ
d
2
−
2
p
T
h
M
+
h
M
T
R
h
M
(1.2.6)
J(\mathop{h}_M)=\mathop{\sigma}_{d}^{2} - 2 \textbf{p} \mathop{}_{}^{T}\mathop{h}_M^{}+\mathop{h}_M^{T} \textbf{R} \mathop{h}_M^{} \tag{1.2.6}
J(hM)=σd2−2pThM+hMTRhM(1.2.6)
对式(1.2.6)的均方误差求偏导数，将得出关于MSE准则最小化的M维线性方程组
∂
[
J
(
h
M
)
]
∂
h
M
=
−
2
p
+
2
R
h
M
(1.2.7)
\cfrac{\partial[J(\mathop{h}_M)]}{\partial\mathop{h}_M} = -2\textbf{p}+2\textbf{R}\mathop{h}_M \tag{1.2.7}
∂hM∂[J(hM)]=−2p+2RhM(1.2.7)
令偏导数为0，即
p
=
R
h
M
o
p
t
(1.2.8)
\textbf{p} = \textbf{R} \mathop{h}_{Mopt} \tag{1.2.8}
p=RhMopt(1.2.8)
上面的式子就是维纳霍夫方程，其中
h
M
o
p
t
\mathop{h}_{Mopt}
hMopt为最佳滤波器权系数矢量，它使得均方误差
J
(
h
M
)
J(\mathop{h}_M)
J(hM)最小，也可以说是使得估计误差的平均功率最小，表达为
h
M
o
p
t
=
R
−
1
p
(1.2.9)
\mathop{h}_{Mopt}=\textbf{R} \mathop{}_{}^{-1} \textbf{p} \tag{1.2.9}
hMopt=R−1p(1.2.9)
二、最小二乘方准则(Least-Square)
上述的维纳滤波器是建立于最小均方误差准则(MMSE)之上，这个准则是基于统计意义上寻求最优滤波的，即需要无限个数据来加权求和，来求得数学期望。然而在实际应用中，只能够获取有限个观测数据，这就需要使用最小二乘方准则(LS)来寻求最优滤波。
①最小二乘方准则与维纳滤波的关系
首先，在维纳滤波(MMSE)估计中，它的代价函数是均方误差信号
J
(
h
M
)
=
E
[
∣
e
(
n
)
∣
2
]
(2.1.1)
J(\mathop{h}_M)=E[\ \vert{e(n)}\vert\mathop{}_{}^{2}] \tag{2.1.1}
J(hM)=E[ ∣e(n)∣2](2.1.1)
其中“离散数学期望E[]”概念如下图(网上找的)
图2 离散数学期望的概念
可以看到在MMSE估计中，需要无限个估计误差数据来加权求和，使得趋于收敛某个数值，这个数值就称为数学期望，而在最小二乘方准则(LS估计)中，代价函数定义为误差数据的有限个的模的平方和，即
J
=
∑
n
=
0
M
−
1
∣
e
(
n
)
∣
2
(2.1.2)
J=\displaystyle\sum_{n=0}^{M-1}\vert{e(n)}\vert\mathop{}_{}^{2} \tag{2.1.2}
J=n=0∑M−1∣e(n)∣2(2.1.2)
通过人们的验证(此处省略！)，可以认为当观测样本数趋于无穷大的时候，最小二乘方法(LS估计)逼近维纳滤波(MMSE估计)；换句话说，最小二乘方法(LS估计)是维纳滤波(MMSE估计)在有限个观测值时的时间平均近似！
②递归最小二乘方法(RLS)
递归最小二乘法是利用迭代的方法，去求解最小二乘准则的最优解，其基本思路是已知 n − 1 n-1 n−1时刻的滤波器权系数矢量，利用当前n时刻新得到的观测数据，用时间迭代的方法计算出 n n n时刻的滤波器权系数矢量。
定义误差为期望序列
d
(
n
)
d(n)
d(n)与滤波器输出值
y
(
n
)
y(n)
y(n)之间的差值
e
M
(
l
,
n
)
=
d
(
l
)
−
h
M
T
(
n
)
X
M
(
l
)
l
=
0
,
1
,
2
,
.
.
.
,
n
(2.2.1)
\mathop{e}_M(l,n)=d(l)-\mathop{h}_M^{T}(n)\mathop{X}_M(l) \qquad \qquad \qquad l=0,1,2,...,n \tag{2.2.1}
eM(l,n)=d(l)−hMT(n)XM(l)l=0,1,2,...,n(2.2.1)
其中第n时刻的滤波器系数矢量定义为
h
M
(
n
)
=
[
h
0
(
n
)
h
1
(
n
)
.
.
.
h
M
−
1
(
n
)
]
T
(2.2.2)
\mathop{h}_M^{}(n)=[\mathop{h}_0^{}(n) \quad \mathop{h}_1^{}(n) \quad ...\quad \mathop{h}_{M-1}^{}(n)]\mathop{}_{}^{T} \tag{2.2.2}
hM(n)=[h0(n)h1(n)...hM−1(n)]T(2.2.2)
下表M表示为滤波器的长度。同样，在第n时刻输入到滤波器的数据矢量表示为
X
M
(
n
)
=
[
x
(
n
)
x
(
n
−
1
)
.
.
.
x
(
n
−
M
+
1
)
]
T
(2.2.3)
\mathop{X}_M^{}(n)=[\mathop{x}_{}^{}(n) \quad \mathop{x}_{}^{}(n-1) \quad ...\quad \mathop{x}_{}^{}(n-M+1)]\mathop{}_{}^{T} \tag{2.2.3}
XM(n)=[x(n)x(n−1)...x(n−M+1)]T(2.2.3)
定义代价函数为误差的模的平方的加权和
J
=
∑
l
=
0
n
λ
n
−
l
∣
e
M
(
l
,
n
)
∣
2
(2.2.4)
J=\displaystyle\sum_{l=0}^{n}\mathop{\lambda}_{}^{n-l}\vert{\mathop{e}_M(l,n)}\vert\mathop{}_{}^{2} \tag{2.2.4}
J=l=0∑nλn−l∣eM(l,n)∣2(2.2.4)
λ
\lambda
λ是范围0
<
<
<
λ
\lambda
λ
≤
\le
≤ 1的加权因子(指数加权)，又叫遗忘因子(Forgetting Factor)。当
λ
\lambda
λ = 1时，表示该代价函数为过去所有估计误差的平均加权(人人有份)。而当0
<
<
<
λ
\lambda
λ
<
<
< 1时，离当前时刻最近的估计误差对代价函数的贡献最大，反之，离当前时刻越远的估计误差对代价函数的贡献越小，这样滤波器系数能够适应时变的数据。
同样地，对式(2.2.4)进行求偏导，并令它为0可得
∑
l
=
0
n
λ
n
−
l
X
M
(
l
)
X
M
T
(
l
)
h
M
(
n
)
=
∑
l
=
0
n
λ
n
−
l
X
M
(
l
)
d
(
l
)
(2.2.5)
\displaystyle\sum_{l=0}^{n}\mathop{\lambda}_{}^{n-l}\mathop{X}_{M}^{}(l)\mathop{X}_{M}^{T}(l)\mathop{h}_{M}^{}(n)=\displaystyle\sum_{l=0}^{n}\mathop{\lambda}_{}^{n-l}\mathop{X}_{M}^{}(l)d(l) \tag{2.2.5}
l=0∑nλn−lXM(l)XMT(l)hM(n)=l=0∑nλn−lXM(l)d(l)(2.2.5)
为了方便，用矢量表示(往后都用矢量表示)。定义信号相关矩阵
R
M
(
n
)
\mathop{R}_{M}^{}(n)
RM(n)
R
M
(
n
)
=
∑
l
=
0
n
λ
n
−
l
X
M
(
l
)
X
M
T
(
l
)
(2.2.6)
\mathop{R}_{M}^{}(n)=\displaystyle\sum_{l=0}^{n}\mathop{\lambda}_{}^{n-l}\mathop{X}_{M}^{}(l)\mathop{X}_{M}^{T}(l) \tag{2.2.6}
RM(n)=l=0∑nλn−lXM(l)XMT(l)(2.2.6)
定义信号互相关矢量
D
M
(
n
)
\mathop{D}_{M}^{}(n)
DM(n)
D
M
(
n
)
=
∑
l
=
0
n
λ
n
−
l
X
M
(
l
)
d
(
l
)
(2.2.7)
\mathop{D}_{M}^{}(n)=\displaystyle\sum_{l=0}^{n}\mathop{\lambda}_{}^{n-l}\mathop{X}_{M}^{}(l)d(l) \tag{2.2.7}
DM(n)=l=0∑nλn−lXM(l)d(l)(2.2.7)
则式(2.2.5)可改为
R
M
(
n
)
h
M
(
n
)
=
D
M
(
n
)
(2.2.8)
\mathop{R}_{M}^{}(n)\mathop{h}_{M}^{}(n)=\mathop{D}_{M}^{}(n) \tag{2.2.8}
RM(n)hM(n)=DM(n)(2.2.8)
所以其最小二乘方最优解为
h
M
(
n
)
=
R
M
−
1
(
n
)
D
M
(
n
)
(2.2.9)
\mathop{h}_{M}^{}(n)=\mathop{R}_{M}^{-1}(n)\mathop{D}_{M}^{}(n) \tag{2.2.9}
hM(n)=RM−1(n)DM(n)(2.2.9)
从式(2.2.9)可以看到，最小二乘方的最优解与信号相关矩阵的逆和信号互相关矢量有关。一般情况下，习惯上在初始化时就要将 R M ( n ) \mathop{R}_{M}^{}(n) RM(n)加上矩阵 δ I M \delta\mathop{I}_{M}^{} δIM，其中 I M \mathop{I}_{M}^{} IM是单位矩阵， δ \delta δ是一个很小的常数，以保证 R M ( n ) \mathop{R}_{M}^{}(n) RM(n)矩阵在迭代的每一步都是非奇异的(为了保证矩阵可逆)。引入 λ n − l \mathop{\lambda}_{}^{n-l} λn−l指数权，可使 δ I M \delta\mathop{I}_{M}^{} δIM增加的影响随时间流逝而消失。
把
R
M
(
n
)
\mathop{R}_{M}^{}(n)
RM(n)按照时间递归拆分为 (把式(2.2.6)拆开就知道了，类似于高中的数列)
R
M
(
n
)
=
λ
R
M
(
n
−
1
)
+
X
M
(
n
)
X
M
T
(
n
)
(2.2.10)
\mathop{R}_{M}^{}(n)=\lambda\mathop{R}_{M}^{}(n-1)+\mathop{X}_{M}^{}(n)\mathop{X}_{M}^{T}(n) \tag{2.2.10}
RM(n)=λRM(n−1)+XM(n)XMT(n)(2.2.10)
因为需要的是
R
M
(
n
)
\mathop{R}_{M}^{}(n)
RM(n)的逆，所以研究人员们动用了黑科技- - -矩阵转化定理：
设 M x M 矩阵 A，B，N x N 矩阵D为正定矩阵，且满足关系
A
=
B
−
1
+
C
D
−
1
C
T
(2.2.11)
A=\mathop{B}_{}^{-1}+\mathop{C}_{}^{}\mathop{D}_{}^{-1}\mathop{C}_{}^{T} \tag{2.2.11}
A=B−1+CD−1CT(2.2.11)
其中矩阵C为 M x N 矩阵，由正定性，矩阵A、B和D是非奇异的。
根据矩阵转化定理，矩阵A的逆为
A
−
1
=
B
−
B
C
(
D
+
C
T
B
C
)
−
1
C
T
B
(2.2.12)
\mathop{A}_{}^{-1}=\mathop{B}_{}^{}-\mathop{B}_{}^{}\mathop{C}_{}^{}(\mathop{D}_{}^{}+\mathop{C}_{}^{T}\mathop{B}_{}^{}\mathop{C}_{}^{})\mathop{}_{}^{-1}\mathop{C}_{}^{T}\mathop{B}_{}^{} \tag{2.2.12}
A−1=B−BC(D+CTBC)−1CTB(2.2.12)
令 A = R M ( n ) B − 1 = λ R M ( n − 1 ) C = X M ( n ) D = 1 (2.2.13) \begin{aligned} A&= \mathop{R}_{M}^{}(n) \\ \mathop{B}_{}^{-1}&= \lambda\mathop{R}_{M}^{}(n-1) \\ C&= \mathop{X}_{M}^{}(n) \\ D&=1 \end{aligned} \tag{2.2.13} AB−1CD=RM(n)=λRM(n−1)=XM(n)=1(2.2.13)
将上式代入式(2.2.12)得到了相关矩阵的逆，为
R
M
−
1
(
n
)
=
1
λ
[
R
M
−
1
(
n
−
1
)
−
R
M
−
1
(
n
−
1
)
X
M
(
n
)
X
M
T
(
n
)
R
M
−
1
(
n
−
1
)
λ
+
X
M
T
(
n
)
R
M
−
1
(
n
−
1
)
X
M
(
n
)
]
(2.2.14)
\mathop{R}_{M}^{-1}(n)=\cfrac{1}{\lambda}[\mathop{R}_{M}^{-1}(n-1)-\cfrac{\mathop{R}_{M}^{-1}(n-1)\mathop{X}_{M}^{}(n)\mathop{X}_{M}^{T}(n)\mathop{R}_{M}^{-1}(n-1)}{\lambda+\mathop{X}_{M}^{T}(n)\mathop{R}_{M}^{-1}(n-1)\mathop{X}_{M}^{}(n)}] \tag{2.2.14}
RM−1(n)=λ1[RM−1(n−1)−λ+XMT(n)RM−1(n−1)XM(n)RM−1(n−1)XM(n)XMT(n)RM−1(n−1)](2.2.14)
为了方便，把
R
M
−
1
(
n
)
\mathop{R}_{M}^{-1}(n)
RM−1(n)定义为
P
M
(
n
)
\mathop{P}_{M}^{}(n)
PM(n)，且再定义一个M维矢量
K
M
(
n
)
\mathop{K}_{M}^{}(n)
KM(n)，称为卡尔曼增益矢量(Kalman gain vector)，其形式为
K
M
(
n
)
=
P
M
(
n
−
1
)
X
M
(
n
)
λ
+
X
M
T
(
n
)
P
M
(
n
−
1
)
X
M
(
n
)
(2.2.15)
\mathop{K}_{M}^{}(n)=\cfrac{\mathop{P}_{M}^{}(n-1)\mathop{X}_{M}^{}(n)}{\lambda+\mathop{X}_{M}^{T}(n)\mathop{P}_{M}^{}(n-1)\mathop{X}_{M}^{}(n)} \tag{2.2.15}
KM(n)=λ+XMT(n)PM(n−1)XM(n)PM(n−1)XM(n)(2.2.15)
利用上面两个定义，把式(2.2.14)改写为
P
M
(
n
)
=
1
λ
[
P
M
(
n
−
1
)
−
K
M
(
n
)
X
M
T
(
n
)
P
M
(
n
−
1
)
]
(2.2.16)
\mathop{P}_{M}^{}(n)=\cfrac{1}{\lambda}[\mathop{P}_{M}^{}(n-1)-\mathop{K}_{M}^{}(n)\mathop{X}_{M}^{T}(n)\mathop{P}_{M}^{}(n-1)] \tag{2.2.16}
PM(n)=λ1[PM(n−1)−KM(n)XMT(n)PM(n−1)](2.2.16)
把上式乘以
X
M
(
n
)
\mathop{X}_{M}^{}(n)
XM(n)，可整理得
P
M
(
n
)
X
M
(
n
)
=
1
λ
[
P
M
(
n
−
1
)
X
M
(
n
)
−
K
M
(
n
)
X
M
T
(
n
)
P
M
(
n
−
1
)
X
M
(
n
)
]
=
K
M
(
n
)
(2.2.17)
\begin{aligned} \mathop{P}_{M}^{}(n)\mathop{X}_{M}^{}(n)&=\cfrac{1}{\lambda}[\mathop{P}_{M}^{}(n-1)\mathop{X}_{M}^{}(n)-\mathop{K}_{M}^{}(n)\mathop{X}_{M}^{T}(n)\mathop{P}_{M}^{}(n-1)\mathop{X}_{M}^{}(n)] \\ &=\mathop{K}_{M}^{}(n) \end{aligned} \tag{2.2.17}
PM(n)XM(n)=λ1[PM(n−1)XM(n)−KM(n)XMT(n)PM(n−1)XM(n)]=KM(n)(2.2.17)
这个式子在后续推导快速RLS算法过程中起到很重要的作用。
回到式(2.2.9)，改写为
h
M
(
n
)
=
P
M
(
n
)
D
M
(
n
)
(2.2.18)
\mathop{h}_{M}^{}(n)=\mathop{P}_{M}^{}(n)\mathop{D}_{M}^{}(n) \tag{2.2.18}
hM(n)=PM(n)DM(n)(2.2.18)
同样地，类似于式(2.2.10)，把
D
M
(
n
)
\mathop{D}_{M}^{}(n)
DM(n)分解成时间递归形式
D
M
(
n
)
=
λ
D
M
(
n
−
1
)
+
d
(
n
)
X
M
(
n
)
(2.2.19)
\mathop{D}_{M}^{}(n)=\lambda\mathop{D}_{M}^{}(n-1)+d(n)\mathop{X}_{M}^{}(n) \tag{2.2.19}
DM(n)=λDM(n−1)+d(n)XM(n)(2.2.19)
将式(2.2.16)和式(2.2.19)代入式(2.2.18)，得出
h
M
(
n
)
=
1
λ
[
P
M
(
n
−
1
)
−
K
M
(
n
)
X
M
T
(
n
)
P
M
(
n
−
1
)
]
∗
[
λ
D
M
(
n
−
1
)
+
d
(
n
)
X
M
(
n
)
]
=
P
M
(
n
−
1
)
D
M
(
n
−
1
)
+
1
λ
d
(
n
)
P
M
(
n
−
1
)
X
M
(
n
)
−
K
M
(
n
)
X
M
T
(
n
)
P
M
(
n
−
1
)
D
M
(
n
−
1
)
−
1
λ
d
(
n
)
K
M
(
n
)
X
M
T
(
n
)
P
M
(
n
−
1
)
X
M
(
n
)
=
h
M
(
n
−
1
)
+
K
M
(
n
)
[
d
(
n
)
−
X
M
T
(
n
)
h
M
(
n
−
1
)
]
(
2.2.20
)
\begin{aligned} \mathop{h}_{M}^{}(n)&=\cfrac{1}{\lambda}[\mathop{P}_{M}^{}(n-1)-\mathop{K}_{M}^{}(n)\mathop{X}_{M}^{T}(n)\mathop{P}_{M}^{}(n-1)]*[\lambda\mathop{D}_{M}^{}(n-1)+d(n)\mathop{X}_{M}^{}(n)] \\ &=\mathop{P}_{M}^{}(n-1)\mathop{D}_{M}^{}(n-1)+\cfrac{1}{\lambda}d(n)\mathop{P}_{M}^{}(n-1)\mathop{X}_{M}^{}(n)-\mathop{K}_{M}^{}(n)\mathop{X}_{M}^{T}(n)\mathop{P}_{M}^{}(n-1)\mathop{D}_{M}^{}(n-1) -\cfrac{1}{\lambda}d(n)\mathop{K}_{M}^{}(n)\mathop{X}_{M}^{T}(n)\mathop{P}_{M}^{}(n-1)\mathop{X}_{M}^{}(n) \\ &=\mathop{h}_{M}^{}(n-1)+\mathop{K}_{M}^{}(n)[d(n)-\mathop{X}_{M}^{T}(n)\mathop{h}_{M}^{}(n-1)] &(2.2.20) \end{aligned}
hM(n)=λ1[PM(n−1)−KM(n)XMT(n)PM(n−1)]∗[λDM(n−1)+d(n)XM(n)]=PM(n−1)DM(n−1)+λ1d(n)PM(n−1)XM(n)−KM(n)XMT(n)PM(n−1)DM(n−1)−λ1d(n)KM(n)XMT(n)PM(n−1)XM(n)=hM(n−1)+KM(n)[d(n)−XMT(n)hM(n−1)](2.2.20)
注意看，这个男人叫。。。，不对，注意到式(2.2.20)右边那项
d
(
n
)
−
X
M
T
(
n
)
h
M
(
n
−
1
)
(2.2.21)
d(n)-\mathop{X}_{M}^{T}(n)\mathop{h}_{M}^{}(n-1) \tag{2.2.21}
d(n)−XMT(n)hM(n−1)(2.2.21)
根据式(2.2.1)中误差的定义，可以把式(2.2.21)定义为
ξ
M
(
n
)
≡
e
M
(
n
,
n
−
1
)
=
d
(
n
)
−
h
M
T
(
n
−
1
)
X
M
(
n
)
(2.2.22)
\mathop{\xi}_{M}^{}(n)≡\mathop{e}_M(n,n-1)=d(n)-\mathop{h}_M^{T}(n-1)\mathop{X}_M(n) \tag{2.2.22}
ξM(n)≡eM(n,n−1)=d(n)−hMT(n−1)XM(n)(2.2.22)
从式子可以看出 h M T ( n − 1 ) X M ( n ) \mathop{h}_M^{T}(n-1)\mathop{X}_M(n) hMT(n−1)XM(n)为自适应滤波器利用第n-1时刻的滤波器权系数去得到第n时刻的输出，所以式(2.2.22)定义为先验估计误差 ξ M ( n ) \mathop{\xi}_{M}^{}(n) ξM(n) （类似于利用过去的经验来验证现在发生的事）
由此得出
h
M
(
n
)
\mathop{h}_M^{}(n)
hM(n)的时间递归公式可以表示为
h
M
(
n
)
=
h
M
(
n
−
1
)
+
K
M
(
n
)
ξ
M
(
n
)
(2.2.23)
\mathop{h}_{M}^{}(n)=\mathop{h}_{M}^{}(n-1)+\mathop{K}_{M}^{}(n)\mathop{\xi}_{M}^{}(n) \tag{2.2.23}
hM(n)=hM(n−1)+KM(n)ξM(n)(2.2.23)
而第n时刻的估计误差为
e
M
(
n
)
≡
e
M
(
n
,
n
)
=
d
(
n
)
−
h
M
T
(
n
)
X
M
(
n
)
(2.2.24)
\mathop{e}_M(n)≡\mathop{e}_M(n,n)=d(n)-\mathop{h}_M^{T}(n)\mathop{X}_M(n) \tag{2.2.24}
eM(n)≡eM(n,n)=d(n)−hMT(n)XM(n)(2.2.24)
比较式(2.2.22)和上式(2.2.24)可以看出，计算 ξ M ( n ) \mathop{\xi}_{M}^{}(n) ξM(n)的时候，第n时刻的滤波器输出是用第n-1时刻的滤波器权系数矢量计算的，而计算 e M ( n ) \mathop{e}_M(n) eM(n)的时候，第n时刻的滤波器输出是用第n时刻的滤波器权系数矢量计算的，因此 e M ( n ) \mathop{e}_M(n) eM(n)也称为后验估计误差。 (先验估计误差与后验估计误差在推导快速RLS算法的过程上也起到很重要的作用)
结合上面推导的式子，整理出RLS算法的步骤：
初始化：
P
M
(
0
)
=
1
/
δ
I
M
\mathop{P}_M(0)=1/\delta\mathop{I}_{M}^{}
PM(0)=1/δIM，其中
δ
\delta
δ是小的正数。
h
M
(
0
)
=
0
\mathop{h}_M^{}(0)=0
hM(0)=0
0
<
λ
≤
1
，
λ
0<\lambda\le1，\lambda
0<λ≤1，λ一般取0.95以上
算法主体：
当
n
=
1
,
2
,
.
.
.
N
n=1,2,...N
n=1,2,...N时，按照步骤完成以下时间迭代运算
①计算滤波器输出
y
(
n
)
=
X
M
T
(
n
)
h
M
(
n
−
1
)
(2.2.25)
y(n)=\mathop{X}_{M}^{T}(n)\mathop{h}_{M}^{}(n-1) \tag{2.2.25}
y(n)=XMT(n)hM(n−1)(2.2.25)
②计算先验误差
ξ
M
(
n
)
=
d
(
n
)
−
y
(
n
)
(2.2.26)
\mathop{\xi}_{M}^{}(n)=d(n)-y(n) \tag{2.2.26}
ξM(n)=d(n)−y(n)(2.2.26)
③计算卡尔曼增益矢量
K
M
(
n
)
=
P
M
(
n
−
1
)
X
M
(
n
)
λ
+
X
M
T
(
n
)
P
M
(
n
−
1
)
X
M
(
n
)
(2.2.27)
\mathop{K}_{M}^{}(n)=\cfrac{\mathop{P}_{M}^{}(n-1)\mathop{X}_{M}^{}(n)}{\lambda+\mathop{X}_{M}^{T}(n)\mathop{P}_{M}^{}(n-1)\mathop{X}_{M}^{}(n)} \tag{2.2.27}
KM(n)=λ+XMT(n)PM(n−1)XM(n)PM(n−1)XM(n)(2.2.27)
④更新相关矩阵的逆
P
M
(
n
)
=
1
λ
[
P
M
(
n
−
1
)
−
K
M
(
n
)
X
M
T
(
n
)
P
M
(
n
−
1
)
]
(2.2.28)
\mathop{P}_{M}^{}(n)=\cfrac{1}{\lambda}[\mathop{P}_{M}^{}(n-1)-\mathop{K}_{M}^{}(n)\mathop{X}_{M}^{T}(n)\mathop{P}_{M}^{}(n-1)] \tag{2.2.28}
PM(n)=λ1[PM(n−1)−KM(n)XMT(n)PM(n−1)](2.2.28)
⑤更新滤波器权系数矢量
h
M
(
n
)
=
h
M
(
n
−
1
)
+
K
M
(
n
)
ξ
M
(
n
)
(2.2.29)
\mathop{h}_{M}^{}(n)=\mathop{h}_{M}^{}(n-1)+\mathop{K}_{M}^{}(n)\mathop{\xi}_{M}^{}(n) \tag{2.2.29}
hM(n)=hM(n−1)+KM(n)ξM(n)(2.2.29)
从式(2.2.23)可以看到，RLS算法(往后称为RLS直接型算法)中的滤波器权系数矢量的时间更新是由卡尔曼增益矢量
K
M
(
n
)
\mathop{K}_{M}^{}(n)
KM(n)的M个元素共同控制的，理论上比起传统LMS算法(基于维纳霍夫方程的近似值)会更快获得收敛，因为后者只有单一的参数- - -步长
μ
\mu
μ用于控制系数调整速率。虽然RLS直接型算法的拥有收敛速度快的优点，但也有它的致命缺点。
RLS直接型算法大约需要20~24位精度才能很好的工作，低于该精度下，RLS直接型算法会由于舍入误差的积聚变得不稳定，滤波器权系数时不时就会发散，导致整个模型崩溃。而后面要讲的快速RLS算法低至11位也能很好地工作，但随着更多的迭代次数，由于舍入误差的积聚，也会像RLS直接型算法一样变得不稳定。对于传统LMS算法来说，它虽然收敛速度远不及前面的两个快，但它对舍入误差具有很强的健壮性。当滤波器权系数精度降低时，虽说整个性能上的退化是非常明显的，但不至于发生不稳定的情况。
由于RLS直接型算法里涉及到矩阵的乘法运算，它的计算复杂度正比于
M
2
\mathop{M}_{}^{2}
M2，当滤波器的长度
M
M
M越长的时候，整个算法的计算复杂度将以平方的级数上升。在后面将要推导的快速RLS算法中，通过利用前向预测和后向预测公式，就可能得出卡尔曼增益矢量的时间迭代方程，从而避免了矩阵乘法运算，使计算复杂度正比于
M
M
M。
总结
维纳霍夫方程是维纳滤波器在基于MMSE准则上求得的最优解的线性方程，而MMSE是基于统计意义上，需要无穷多的数据去确立维纳霍夫方程的解，这显然是不实际的，而最小二乘方(LS)使用的是确定性思想，通过有限个观测数据去求维纳霍夫方程的近似解(当观测数据足够多的时候，两者结果趋于一致)。
在最小二乘方(LS)求最优解的问题上，使用了时间迭代的方法- - -RLS直接型算法，但RLS直接型算法的复杂度为
O
(
M
2
)
O(\mathop{M}_{}^{2})
O(M2),且对精度较低带来的舍入误差的积聚很敏感。于是提出了快速RLS算法，其复杂度为
O
(
7
M
)
O(7M)
O(7M)，虽然复杂度相对于RLS直接型算法降低了很多，但仍存在对舍入误差敏感，并且显示出不稳定的问题。针对上述问题，又提出了使快速RLS算法稳定的修改算法，称为稳定快速RLS算法(SFTRLS),其复杂度为
O
(
8
M
)
O(8M)
O(8M)。
下个系列将重点放在对快速RLS算法以及稳定快速RLS算法的推导过程的剖析 上，并附加上Matlab代码。
5658
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
