维纳滤波器
本文主要总结下维纳滤波器的相关知识,属于现代数字信号处理的内容
维纳滤波问题
已知 y ( n ) y(n) y(n)是参考信号,即我们所期望的输出信号, x ( n ) x(n) x(n)是输入信号, e ( n ) e(n) e(n)是误差信号。 y ( n ) y(n) y(n)和 x ( n ) x(n) x(n)是均值为0的平稳的离散时间信号,二阶矩(即自相关、互相关)已知。该滤波器是线性的(FIR、IIR)。
采用的准则:最小均方误差(MMSE)
m
i
n
i
m
i
z
e
E
[
e
2
(
n
)
]
=
E
{
[
y
(
n
)
−
y
~
(
n
)
]
2
}
(1)
minimize~E[e^{2}(n)]=E\{[y(n)-\tilde{y}(n)]^{2}\} \tag{1}
minimize E[e2(n)]=E{[y(n)−y~(n)]2}(1)
Wenier-Hopf方程
y
~
\tilde{y}
y~是我们滤波器的输出信号,即估计信号。
y
~
(
n
)
=
∑
i
h
i
x
(
n
−
i
)
(2)
\tilde{y}(n)=\sum_{i}h_ix(n-i) \tag{2}
y~(n)=i∑hix(n−i)(2)
采用准则最小化均方误差:
m
i
n
i
m
i
z
e
J
=
E
[
e
2
(
n
)
]
=
E
{
[
y
(
n
)
−
y
~
(
n
)
]
2
}
(3)
minimize~J=E[e^{2}(n)]=E\{[y(n)-\tilde{y}(n)]^{2}\} \tag{3}
minimize J=E[e2(n)]=E{[y(n)−y~(n)]2}(3)
我们的目的是求出最优的滤波器系数
h
i
h_i
hi,所以要对
h
i
h_i
hi求偏导:
∂
J
∂
h
i
=
∂
E
[
e
2
(
n
)
]
∂
h
i
=
2
E
[
e
(
n
)
∂
e
(
n
)
∂
h
i
]
(4)
\frac{\partial{J}}{\partial{h_i}}=\frac{\partial{E[e^2(n)]}}{\partial{h_i}}=2E[e(n)\frac{\partial{e(n)}}{\partial{h_i}}] \tag{4}
∂hi∂J=∂hi∂E[e2(n)]=2E[e(n)∂hi∂e(n)](4)
又因为
e
(
n
)
=
y
(
n
)
−
y
~
(
n
)
e(n)=y(n)-\tilde{y}(n)
e(n)=y(n)−y~(n),再代入公式(2)后,我们令偏导数得0,得到:
∂
J
∂
h
i
=
−
2
E
[
e
(
n
)
x
(
n
−
i
)
]
=
0
⟹
E
[
e
(
n
)
x
(
n
−
i
)
]
=
0
(5)
\frac{\partial{J}}{\partial{h_i}}=-2E[e(n)x(n-i)]=0\\ \implies ~~E[e(n)x(n-i)]=0 \tag{5}
∂hi∂J=−2E[e(n)x(n−i)]=0⟹ E[e(n)x(n−i)]=0(5)
因为
e
(
n
)
e(n)
e(n)展开后有变量
i
i
i,所以我们将(5)式中的
i
i
i换成
j
j
j,得
E
[
y
(
n
)
x
(
n
−
j
)
−
∑
i
h
i
x
(
n
−
i
)
x
(
n
−
j
)
]
=
0
(6)
E[y(n)x(n-j)-\sum_i{h_ix(n-i)x(n-j)}]=0 \tag{6}
E[y(n)x(n−j)−i∑hix(n−i)x(n−j)]=0(6)
令
r
y
x
(
j
)
=
E
[
y
(
n
)
y
(
n
−
j
)
]
r
x
x
(
j
)
=
E
[
x
(
n
)
x
(
n
−
j
)
]
(7)
r_{yx}(j)=E[y(n)y(n-j)]\\ r_{xx}(j)=E[x(n)x(n-j)] \tag{7}
ryx(j)=E[y(n)y(n−j)]rxx(j)=E[x(n)x(n−j)](7)
则(6)式可以被重新写成:
r
y
x
(
j
)
=
∑
i
h
i
r
x
x
(
j
−
i
)
.
(8)
r_{yx}(j)=\sum_{i}{h_ir_{xx}(j-i)}. \tag{8}
ryx(j)=i∑hirxx(j−i).(8)
公式(8)被称为Weiner-Hopf方程
正交原理
在Weiner-Hopf方程的推导过程中我们得知两个比较重要的信息:
-
E [ e ( n ) x ( n − j ) ] = 0 E[e(n)x(n-j)]=0 E[e(n)x(n−j)]=0
-
E [ e ( n ) y ~ ( n − j ) ] = 0 E[e(n)\tilde{y}(n-j)]=0 E[e(n)y~(n−j)]=0
即线性最优滤波器(维纳滤波器)的充要条件是滤波器的输出与误差正交
N阶FIR维纳滤波器的求解
假设我们的输入信号序列为:
X
(
n
)
=
[
x
(
n
)
,
x
(
n
−
1
)
,
.
.
.
x
(
n
−
N
+
1
)
]
T
(9)
\mathbf{X}(n)=[x(n),x(n-1),...x(n-N+1)]^T \tag{9}
X(n)=[x(n),x(n−1),...x(n−N+1)]T(9)
N阶的滤波器系数为:
H
=
[
h
0
,
h
1
,
h
2
,
.
.
.
h
N
−
1
]
T
(10)
\mathbf{H}=[h_0,h_1,h_2,...h_{N-1}]^T \tag{10}
H=[h0,h1,h2,...hN−1]T(10)
滤波器的输出为:
y
~
(
n
)
=
X
T
(
n
)
H
=
H
T
X
(
n
)
(11)
\tilde{y}(n)=\mathbf{X}^{T}(n)\mathbf{H}=\mathbf{H}^{T}\mathbf{X}(n) \tag{11}
y~(n)=XT(n)H=HTX(n)(11)
求最优滤波:
E
[
e
(
n
)
x
(
n
−
j
)
]
=
0
⟹
E
[
e
(
n
)
X
(
n
)
]
(12)
E[e(n)x(n-j)]=0\implies E[e(n)\mathbf{X}(n)]\tag{12}
E[e(n)x(n−j)]=0⟹E[e(n)X(n)](12)
E { [ y ( n ) − X T ( n ) H ] X ( n ) } = 0 E { [ y ( n ) − H T X ( n ) ] X ( n ) } = 0 E [ y ( n ) X ( n ) ] − E [ X ( n ) X T ( n ) ] H = 0 (13) E\{[y(n)-\mathbf{X}^{T}(n)\mathbf{H}]\mathbf{X}(n)\}=0\\ E\{[y(n)-\mathbf{H}^{T}\mathbf{X}(n)]\mathbf{X}(n)\}=0\\ E[y(n)\mathbf{X}(n)]-E[\mathbf{X(n)\mathbf{X^{T}(n)}}]\mathbf{H}=0\tag{13} E{[y(n)−XT(n)H]X(n)}=0E{[y(n)−HTX(n)]X(n)}=0E[y(n)X(n)]−E[X(n)XT(n)]H=0(13)
所以有:
r
y
x
=
R
x
x
H
o
p
t
⟹
H
o
p
t
=
R
x
x
−
1
r
y
x
(14)
r_{yx}=R_{xx}\mathbf{H}_{opt} \implies \mathbf{H}_{opt}=R_{xx}^{-1}r_{yx}\tag{14}
ryx=RxxHopt⟹Hopt=Rxx−1ryx(14)
N阶FIR维纳滤波器的最小均方误差 J m i n J_{min} Jmin
公式(3)展示了均方误差的计算方法,而公式(11)是N阶FIR维纳滤波器的输出输入关系,将(11)代入(3)可以计算出:
J
=
E
[
y
2
(
n
)
−
2
y
(
n
)
y
~
(
n
)
+
y
~
2
(
n
)
]
=
E
[
y
2
(
n
)
]
−
2
E
[
y
(
n
)
y
~
(
n
)
]
+
E
[
y
~
2
(
n
)
]
=
E
[
y
2
(
n
)
]
−
2
E
[
y
(
n
)
X
T
(
n
)
H
]
+
E
[
H
T
X
(
n
)
X
T
(
n
)
H
]
=
E
[
y
2
(
n
)
]
−
2
r
y
x
T
H
+
H
T
R
x
x
H
(15)
\begin{aligned} J&=E[y^2(n)-2y(n)\tilde{y}(n)+\tilde{y}^2(n)]\\ &=E[y^2(n)]-2E[y(n)\tilde{y}(n)]+E[\tilde{y}^2(n)]\\ &=E[y^2(n)]-2E[y(n)\mathbf{X}^T(n)\mathbf{H}]+E[\mathbf{H}^T\mathbf{X}(n)\mathbf{X}^{T}(n)\mathbf{H}]\\ &=E[y^2(n)]-2r_{yx}^T\mathbf{H}+\mathbf{H}^TR_{xx}\mathbf{H}\tag{15} \end{aligned}
J=E[y2(n)−2y(n)y~(n)+y~2(n)]=E[y2(n)]−2E[y(n)y~(n)]+E[y~2(n)]=E[y2(n)]−2E[y(n)XT(n)H]+E[HTX(n)XT(n)H]=E[y2(n)]−2ryxTH+HTRxxH(15)
在最优的滤波情况下,均方根误差达到最小,所以:
J
m
i
n
=
E
[
y
2
(
n
)
]
−
2
r
y
x
T
H
o
p
t
+
H
o
p
t
T
R
x
x
H
o
p
t
=
E
[
y
2
(
n
)
]
−
2
H
o
p
t
T
R
x
x
T
H
o
p
t
+
H
o
p
t
T
R
x
x
H
o
p
t
=
E
[
y
2
(
n
)
]
−
2
H
o
p
t
T
R
x
x
H
o
p
t
+
H
o
p
t
T
R
x
x
H
o
p
t
=
E
[
y
2
(
n
)
]
−
H
o
p
t
T
R
x
x
H
o
p
t
=
E
[
y
2
(
n
)
]
−
H
o
p
t
T
r
y
x
(16)
\begin{aligned} J_{min}& =E[y^2(n)]-2r_{yx}^T\mathbf{H}_{opt}+\mathbf{H}_{opt}^TR_{xx}\mathbf{H}_{opt}\\ & =E[y^2(n)]-2\mathbf{H}_{opt}^{T}R_{xx}^T\mathbf{H}_{opt}+\mathbf{H}_{opt}^TR_{xx}\mathbf{H}_{opt}\\ & =E[y^2(n)]-2\mathbf{H}_{opt}^{T}R_{xx}\mathbf{H}_{opt}+\mathbf{H}_{opt}^TR_{xx}\mathbf{H}_{opt}\\ & =E[y^2(n)]-\mathbf{H}_{opt}^{T}R_{xx}\mathbf{H}_{opt}\\ & =E[y^2(n)]-\mathbf{H}_{opt}^{T}r_{yx}\tag{16} \end{aligned}
Jmin=E[y2(n)]−2ryxTHopt+HoptTRxxHopt=E[y2(n)]−2HoptTRxxTHopt+HoptTRxxHopt=E[y2(n)]−2HoptTRxxHopt+HoptTRxxHopt=E[y2(n)]−HoptTRxxHopt=E[y2(n)]−HoptTryx(16)
N阶FIR维纳滤波器的误差性能曲面
由公式(15)我们可以得知,维纳滤波器的误差性能曲面是一个二次曲面: J = E [ y 2 ( n ) ] − 2 r y x T H + H T R x x H J=E[y^2(n)]-2r_{yx}^T\mathbf{H}+\mathbf{H}^TR_{xx}\mathbf{H} J=E[y2(n)]−2ryxTH+HTRxxH
在最优的情况下: J m i n = E [ y 2 ( n ) ] − H o p t T r y x J_{min}=E[y^2(n)]-\mathbf{H}_{opt}^{T}r_{yx} Jmin=E[y2(n)]−HoptTryx,说明该误差性能曲面有最低点 J m i n J_{min} Jmin。维纳滤波自适应的过程是自动调整系数,使均方误差达到最小值 J m i n J_{min} Jmin的过程,最常用的方法是梯度下降法和LMS算法,这里不做介绍。
下面我们来证明N阶FIR维纳滤波器有解,并且有唯一解:
计算:
J
(
H
)
−
J
m
i
n
=
E
[
y
2
(
n
)
]
−
2
r
y
x
T
H
+
H
T
R
x
x
H
−
{
E
[
y
2
(
n
)
]
−
H
o
p
t
T
r
y
x
}
=
−
2
r
y
x
T
H
+
H
T
R
x
x
H
+
H
o
p
t
T
r
y
x
(17)
\begin{aligned} J(\mathbf{H})-J_{min}&=E[y^2(n)]-2r_{yx}^T\mathbf{H}+\mathbf{H}^TR_{xx}\mathbf{H}-\{E[y^2(n)]-\mathbf{H}_{opt}^{T}r_{yx}\}\\ &=-2r_{yx}^T\mathbf{H}+\mathbf{H}^TR_{xx}\mathbf{H}+\mathbf{H}_{opt}^{T}r_{yx}\tag{17} \end{aligned}
J(H)−Jmin=E[y2(n)]−2ryxTH+HTRxxH−{E[y2(n)]−HoptTryx}=−2ryxTH+HTRxxH+HoptTryx(17)
公式(17)中等式右边的第一项
−
2
r
y
x
T
H
-2r_{yx}^T\mathbf{H}
−2ryxTH可以进行拆分:
−
2
r
y
x
T
H
=
−
r
y
x
T
H
−
r
y
x
T
H
=
−
H
o
p
t
T
R
x
x
H
−
H
T
r
y
x
=
−
H
o
p
t
T
R
x
x
H
−
H
T
R
x
x
H
o
p
t
(18)
\begin{aligned} -2r_{yx}^T\mathbf{H}&=-r_{yx}^T\mathbf{H}-r_{yx}^T\mathbf{H}\\ &=-\mathbf{H}^{T}_{opt}R_{xx}\mathbf{H}-\mathbf{H}^{T}r_{yx}\\ &=-\mathbf{H}^{T}_{opt}R_{xx}\mathbf{H}-\mathbf{H}^{T}R_{xx}\mathbf{H}_{opt}\tag{18} \end{aligned}
−2ryxTH=−ryxTH−ryxTH=−HoptTRxxH−HTryx=−HoptTRxxH−HTRxxHopt(18)
公式(17)中的等式有低钠最后一项
H
o
p
t
T
r
y
x
\mathbf{H}_{opt}^{T}r_{yx}
HoptTryx可以写成
H
o
p
t
T
r
y
x
=
H
o
p
t
T
R
x
x
H
o
p
t
(19)
\mathbf{H}_{opt}^{T}r_{yx}=\mathbf{H}_{opt}^{T}R_{xx}\mathbf{H}_{opt}\tag{19}
HoptTryx=HoptTRxxHopt(19)
将公式(18)(19)代入到公式(17)中去,得:
J
(
H
)
−
J
m
i
n
=
−
H
o
p
t
T
R
x
x
H
−
H
T
R
x
x
H
o
p
t
+
H
T
R
x
x
H
+
H
o
p
t
T
R
x
x
H
o
p
t
=
H
o
p
t
T
R
x
x
(
H
o
p
t
−
H
)
+
H
T
R
x
x
(
H
−
H
o
p
t
)
=
(
H
T
−
H
o
p
t
T
)
R
x
x
(
H
−
H
o
p
t
)
=
(
H
−
H
o
p
t
)
T
R
x
x
(
H
−
H
o
p
t
)
(20)
\begin{aligned} J(\mathbf{H})-J_{min}&=-\mathbf{H}^{T}_{opt}R_{xx}\mathbf{H}-\mathbf{H}^{T}R_{xx}\mathbf{H}_{opt}+\mathbf{H}^TR_{xx}\mathbf{H}+\mathbf{H}_{opt}^{T}R_{xx}\mathbf{H}_{opt}\\ &=\mathbf{H}^{T}_{opt}R_{xx}(\mathbf{H}_{opt}-\mathbf{H})+\mathbf{H}^TR_{xx}(\mathbf{H}-\mathbf{H}_{opt})\\ &=(\mathbf{H}^T-\mathbf{H}^{T}_{opt})R_{xx}(\mathbf{H}-\mathbf{H}_{opt})\\ &=(\mathbf{H}-\mathbf{H}_{opt})^{T}R_{xx}(\mathbf{H}-\mathbf{H}_{opt})\tag{20} \end{aligned}
J(H)−Jmin=−HoptTRxxH−HTRxxHopt+HTRxxH+HoptTRxxHopt=HoptTRxx(Hopt−H)+HTRxx(H−Hopt)=(HT−HoptT)Rxx(H−Hopt)=(H−Hopt)TRxx(H−Hopt)(20)
所以有:
J
(
H
)
=
J
m
i
n
+
(
H
−
H
o
p
t
)
T
R
x
x
(
H
−
H
o
p
t
)
(21)
J(\mathbf{H})=J_{min}+(\mathbf{H}-\mathbf{H}_{opt})^{T}R_{xx}(\mathbf{H}-\mathbf{H}_{opt})\tag{21}
J(H)=Jmin+(H−Hopt)TRxx(H−Hopt)(21)
证明完毕,说明N阶FIR维纳滤波器有解并且有唯一解。
之后对问题就是如何通过最陡下降法和LMS算法来在误差性能曲面上来求最优解。