第五章 维纳滤波
(Wiener Filtering)
随机信号或随机过程(random process)是普遍存在的。一方面,任何确定性信号经过测
量后往往就会引入随机性误差而使该信号随机化;另一方面,任何信号本身都存在随机干扰,
通常把对信号或系统功能起干扰作用的随机信号称之为噪声。噪声按功率谱密度划分可以分
为白噪声(white noise)和色噪声(color noise),我们把均值为 0 的白噪声叫纯随机信
号(pure random signal)。因此,任何其它随机信号都可看成是纯随机信号与确定性信号
并存的混合随机信号或简称为随机信号。要区别干扰(interference)和噪声( noise)两种
事实和两个概念。非目标信号(nonobjective signal)都可叫干扰。干扰可以是确定信号,
如国内的 50Hz 工频干扰。干扰也可以是噪声,纯随机信号(白噪声)加上一个直流成分(确
定性信号),就成了最简单的混合随机信号。医学数字信号处理的目的是要提取包含在随机
信号中的确定成分,并探求它与生理、病理过程的关系,为医学决策提供一定的依据。例如
从自发脑电中提取诱发脑电信号,就是把自发脑电看成是干扰信号,从中提取出需要的信息
成分。因此我们需要寻找一种最佳线性滤波器,当信号和干扰以及随机噪声同时输入该滤波
器时,在输出端能将信号尽可能精确地表现出来。
维纳滤波和卡尔曼滤波就是用来解决这样一类问题的方法:从噪声中提取出有用的信
号。实际上,这种线性滤波方法也被看成是一种估计问题或者线性预测问题。
设有一个线性系统,它的单位脉冲响应是
)(nh
,当输入一个观测到的随机信号
)(nx
,
简称观测值,且该信号包含噪声
)(nw
和有用信号
)(ns
,简称信号,也即
nx
)(
=
ns
)(
+
nw
)(
(5-1)
则输出
)(ny
为
ny
)(
=
nx
)(
∗
nh
)(
=
∑+∞
mnxmh
(
()
−
)
(5-2)
m
−∞=
我们希望输出得到的
)(ny
与有用信号
)(ns
尽量接近,因此称
)(ny
为
)(ns
的估计值,用
)(ˆ ns
来表示
)(ny
,我们就有了维纳滤波器的系统框图,如图 5.1。这个系统地单位脉冲响
应也称为对于
)(ns
的一种估计器。
nx
)(
=
ns
)(
+
nw
)(
)(nh
ny
)(
=
)(ˆ
ns
图 5.1 维纳滤波器的输入输出关系
如果该系统是因果系统,式(5-2)的 m=0,1,2,…,则输出的
)(ˆ ns
可以看成是由当前
时刻的观测值
)(nx
和过去时刻的观测值
( −nx
)1
、
( −nx
)2
、
( −nx
)3
…的估计值。用当
前的和过去的观测值来估计当前的信号
ny
)(
=
)(ˆ
ns
称为滤波;用过去的观测值来估计当前
的或将来的信号
ny
)(
=
(ˆ
Nns
+
)
,N ,称为预测;用过去的观测值来估计过去的信号
0≥
ny
)(
=
(ˆ
Nns
−
)
,N ,称为平滑或者内插。本章将讨论滤波和预测问题。
1≥
从图 5.1 的系统框图中估计到的
)(ˆ ns
信号和我们期望得到的有用信号
)(ns
不可能完全
相同,这里用
)(ne
来表示真值和估计值之间的误差
ne
)(
=
ns
)(
−
)(ˆ
ns
(5-3)
显然
)(ne
是随机变量,维纳滤波和卡尔曼滤波的误差准则就是最小均方误差准则
[
neE
)(
2
]
=
[
nsE
)((
−
(ˆ
ns
))
]2
(5-4)
维纳滤波和卡尔曼滤波都是解决线性滤波和预测问题的方法,并且都是以均方误差最小
为准则的,在平稳条件下两者的稳态结果是一致的。但是它们解决问题的方法有很大区别。
维纳滤波是根据全部过去观测值和当前观测值来估计信号的当前值,因此它的解形式是系统
的传递函数
)(zH
或单位脉冲响应
)(nh
;卡尔曼滤波是用当前一个估计值和最近一个观测
值来估计信号的当前值,它的解形式是状态变量值。
维纳滤波只适用于平稳随机过程,卡尔曼滤波就没有这个限制。设计维纳滤波器要求已
知信号与噪声的相关函数,设计卡尔曼滤波要求已知状态方程和量测方程,当然两者之间也
有联系。
第一节 维纳滤波器的时域解(Time domain solution of
the Wiener filter)
设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应
)(nh
或传递
函数
)(zH
的表达式,其实质就是解维纳-霍夫(Wiener-Hopf)方程。我们从时域入手求
最小均方误差下的
)(nh
,用
hopt
)(n
表示最佳线性滤波器。这里只讨论因果可实现滤波器的
设计。
5.1.1 因果的维纳滤波器
设
)(nh
是物理可实现的,也即是因果序列:
nh
)(
=
,0
n
当
<
0
因此,从式(5-1)、(5-2)、(5-3)、(5-4)推导:
∑+∞
ny
)(
)(ˆ
ns
=
=
m
=
0
mnxmh
(
()
−
)
(5-5)
[
neE
)(
2
]
=
nsE
)((
⎡
⎢
⎣
−
∑+∞
m
=
0
mnxmh
(
()
−
2
))
⎤
⎥
⎦
(5-6)
要使得均方误差最小,则将上式对各
(mh
)
,m=0,1,…,求偏导,并且等于零,得:
2
nsE
)((
⎡
⎢
⎣
即
− ∑+∞
m
=
0
nxmnxmh
())
opt
()
−
(
−
j
)
⎤
=⎥
⎦
0
j
=
L2,1,0
(5-7)
[
nxnsE
()(
−
j
)
]
=
∑+∞
m
=
0
nxmnxEmh
()
opt
−
(
)
(
[
−
j
)
]
j
≥
0
(5-8)
用相关函数 R 来表达上式,则得到维纳-霍夫方程的离散形式:
R
xs
j
)(
= ∑+∞
m
=
0
Rmh
opt
)
(
(
mj
−
)
xx
j
≥
0
(5-9)
从维纳-霍夫方程中解出的 h 就是最小均方误差下的最佳 h,
hopt
)(n
。求到
hopt
)(n
,这时
的均方误差为最小:
[
neE
)(
2
]
min
=
nsE
)((
⎡
⎢
⎣
−
∑+∞
m
=
0
mnxmh
opt
()
−
(
2
))
⎤
⎥
⎦
=
2
nsE
[
ns
)(2)(
−
+∞
∑
m
=
0
mnxmh
(
()
−
)
+
=
R
ss
2)0(
−
+∞
∑
m
=
0
mRmh
opt
)
(
(
xs
)
+
+∞
+∞
∑∑
m
=
0
r
=
0
hmnxmh
)
opt
()
−
(
opt
rnxr
()(
−
])
+∞
∑
m
=
0
mh
opt
(
)
⎡
⎢
⎣
+∞
∑
r
=
0
h
opt
rmRr
)(
−
(
xx
)
⎤
⎥
⎦
由式(5-9)进一步化简得:
2
neE
[
(
)]
min
=
R
ss
)0(
−
∑+∞
m
=
0
mRmh
opt
(
)
(
xs
)
(5-10)
5.1.2 有限脉冲响应法求解维纳-霍夫方程
如何去求解维纳-霍夫方程,即式(5-9)中解
hopt
)(n
的问题,设
)(nh
是一个因果序
列且可以用有限长(N 点长)的序列去逼进它,则式(5-5) -(5-10)分别发生变化:
ny
)(
=
)(ˆ
ns
=
1
N
∑−
m
=
0
mnxmh
(
()
−
)
(5-11)
[
neE
)(
2
]
=
nsE
)((
⎡
⎢
⎣
−
1
N
∑−
m
=
0
mnxmh
(
()
−
2
))
⎤
⎥
⎦
(5-12)
2
nsE
)((
⎡
⎢
⎣
1
N
− ∑−
m
=
0
nxmnxmh
())
opt
()
−
(
−
j
)
⎤
=⎥
⎦
0
j
=
2,1,0
L
N
−
1
(5-13)
[
nxnsE
()(
−
j
)
]
=
1
N
∑−
m
=
0
nxmnxEmh
()
opt
−
(
(
)
[
−
j
)
]
j
=
,1,0
L
N
−
1
(5-14)
R
xs
j
)(
1
N
= ∑−
m
=
0
于是得到 N 个线性方程:
Rmh
opt
)
(
(
mj
−
)
xx
j
=
,2,1,0
,
N
−
1
L
(5-15)
⎧
⎪
⎪
⎨
⎪
⎪
⎩
j
j
=
=
0
1
R
R
xs
xs
)0(
)1(
=
=
h
R
)0(
R
h
)0(
Rh
)0(
)1(
+
xx
Rh
)1(
)1(
+
)1(
xx
)0(
xx
xx
+
+
+
L
+
L
Nh
(
Nh
(
NR
(
)1
−
xx
NR
)1
(
−
xx
−
−
)1
)2
M
Nj
=
−
1
NR
xs
(
−
)1
=
h
)0(
NR
xx
(
)1
+−
M
NRh
)1(
(
xx
−
)2
+
L
+
Nh
(
−
)1
R
)0(
xx
写成矩阵形式有:
Rxx=
⎡
⎢
⎢
⎢
⎢
⎣
R
xx
R
xx
)0(
)1(
R
R
xx
xx
)1(
)0(
M
NR
xx
(
M
NR
xx
(
)2
只要Rxx是非奇异的,就可以求到H:
)1
−
−
L
L
L
L
NR
(
xx
NR
(
xx
−
−
)1
)2
M
)0(
xx
R
⎤
⎥
⎥
⎥
⎥
⎦
,是自相关矩阵。
⎡
⎢
⎢
⎢
⎢
⎣
R
xx
R
xx
)0(
)1(
R
R
xx
xx
)1(
)0(
M
NR
xx
简化形式:
(
−
)1
M
NR
xx
(
−
)2
L
L
L
L
NR
(
xx
NR
(
xx
−
−
)1
)2
M
)0(
xx
R
⎤
⎥
⎥
⎥
⎥
⎦
⎡
⎢
⎢
⎢
⎢
⎣
h
)0(
h
)1(
M
Nh
(
−
)1
⎤
⎥
⎥
⎥
⎥
⎦
=
⎡
⎢
⎢
⎢
⎢
⎣
R
xs
R
xs
)0(
)1(
M
NR
xs
(
−
)1
⎤
⎥
⎥
⎥
⎥
⎦
(5-16)
式中,H=[h(0) h(1) …h(N-1)]′,是待求的单位脉冲响应;
RxxH=Rxs (5-17)
Rxs=[
R
),0(
R
xs
),1(
xs
L
−NR
xs
(
])1
′,是互相关序列;
H=Rxx
-1Rxs (5-18)
求得
hopt
)(n
后,这时的均方误差为最小:
[
neE
)(
2
]
min
=
nsE
)((
⎡
⎢
⎣
−
1
N
∑−
m
=
0
mnxmh
opt
()
−
(
2
))
⎤
⎥
⎦
=
2
nsE
[
ns
)(2)(
−
N
1
−
∑
m
=
0
mnxmh
(
()
−
)
+
N
1
−
N
1
−
∑∑
m
=
0
r
=
0
hmnxmh
)
opt
()
−
(
opt
rnxr
()(
−
])
=
R
ss
2)0(
−
N
1
−
∑
m
=
0
mRmh
opt
)
(
(
xs
)
+
N
1
−
∑
m
=
0
mh
opt
(
)
⎡
⎢
⎣
N
1
−
∑
r
=
0
h
opt
rmRr
)(
−
(
xx
)
⎤
⎥
⎦
由式(5-15)进一步化简得:
2
neE
[
(
)]
min
=
R
ss
)0(
−
1
N
∑−
m
=
0
mRmh
opt
(
)
(
xs
)
(5-19)
用有限长的
)(nh
来实现维纳滤波时,当已知观测值的自相关和观测值与信号的互相关
时就可以按照式(5-15)在时域里求解
hopt
)(n
。但是当 N 比较大时,计算量很大,并且涉
及到求自相关矩阵的逆矩阵问题。注意到式(5-15)的表现形式和第三章的 AR 模型参数估
计的矩阵形式类似,因而也可以用前面介绍的 L-D 快速算法实现求解。
若信号
)(ns
与噪声
)(nw
互不相关,即,
则有
mR
sw
(
)
=
mR
ws
(
)
=
0
mnsnxEmR
xs
()([
+
=
)
(
)]
=
mnsnwmnsnsE
()([
()(
+
+
+
)
)]
=
mR
ss
(
)
mR
xx
(
)
=
nwnsE
[(
(
)(
+
))(
则式(5-15)和式(5-19)化为:
mnwmns
(
+
+
+
(
)
))]
=
mRmR
ss
+
ww
(
)
(
)
jR
)(
ss
1
N
= ∑−
m
=
0
RmjRmh
opt
)[
+
−
(
(
)
ss
(
mj
−
)]
ww
j
=
,2,1,0
L
,
N
−
1
(5-20)
2
neE
[
(
)]
min
=
R
ss
)0(
−
1
N
∑−
m
=
0
mRmh
opt
(
)
(
ss
)
(5-21)
【例 5-1】已知图 5.1 中
nx
)(
=
ns
)(
+
nw
)(
,且
)(ns
与
)(nw
统计独立,其中
)(ns
的
自相关序列为
ss mR
(
)
=
6.0
m
,
)(nw
是方差为 1 的单位白噪声,试设计一个 N=2 的维纳
滤波器来估计
)(ns
,并求最小均方误差。
解:依题意,已知信号的自相关和噪声的自相关为:
ss mR
(
)
=
6.0
m
,
mRww
(
)
(
δ=
m
)
,
代入式(5-20)得
j
j
=
=
⎧
⎨
⎩
0
=
6.01
=
h
h
)1(6.0)0(21
h
)1(2)0(6.0
+
+
h
解得:
)0(h
=0.451,
)1(h
=0.165。
将上述结果代入式(5-21),求得最小均方误差:
2
neE
[
(
)]
min
=
R
ss
)0(
−
1
∑
m
=
0
mRmh
(
(
)
ss
1)
−=
h
)1(6.0)0(
−
h
=
45.0
若要进一步减小误差可以适当增加维纳滤波的阶数,但相应的计算量也会增加。
5.1.3 预白化法求解维纳-霍夫方程
从上面分析知求解维纳-霍夫方程比较复杂,本节用波德(Bode)和香农(Shannon)
提出的白化的方法求解维纳-霍夫方程,得到系统函数
)(zH
。
由第三章的知识,我们知道随机信号都可以看成是由一白色噪声
)(1 nw
激励一个物理可
实现的系统或模型的响应,如图 5.2 所示,其中 A(z)表示系统的传递函数。由于
nx
)(
=
ns
)(
+
nw
)(
,在图 5.2 的基础上给出
)(nx
的信号模型,图 5.3 所示。把这两个模型
合并最后得到维纳滤波器的信号模型,图 5.4 所示,其中传递函数用 B(z)表示。
)(1 nw
)(zA
)(ns
)(1 nw
)(zA
)(ns
)(nw
⊕
)(nx
图 5.2
)(ns
的信号模型 图 5.3
)(nx
的信号模型
)(1 nw
)(zB
)(nx
白噪声的自相关函数为
,它的 z 变换就等于 。图 5.2 中输出
2
1wσ
图 5.4 维纳滤波器的输入信号模型
R
(
δσ=
m
m
)
)
(
ww
11
2
w
1
信号的自相关函数为
(mRss
)
,根据卷积性质有
mR
ss
(
)
=
mnsnsE
()([
+
)]
=
E
[
+∞
∑
k
−∞=
knwka
)(
−
(
)
⋅
1
+∞
∑
r
−∞=
rmnwra
)(
−
+
(
])
1
+∞
= ∑
k
−∞=
ka
)(
+∞
∑
r
−∞=
Rra
)(
(
km
−+
r
)
ww
11
令
l
−=
r
k
,
上式
=
+∞
∑
k
−∞=
ka
)(
+∞
∑
−∞=
l
ka
(
+
Rl
)
ww
11
(
lm
−
)
=
+∞
∑
−∞=
l
R
ww
11
(
lm
−
)
+∞
∑
k
−∞=
kaka
()(
+
l
)
令
f
l
)(
= ∑+∞
l
−∞=
kaka
()(
+
l
)
=
la
)(
la
(
−∗
)
,代入上式得
mR
ss
(
)
= ∑+∞
l
−∞=
R
ww
11
(
lm
−
)
f
l
)(
=
R
(
ww
11
mfm
∗
(
)
)
=
R
(
ww
11
mamam
(
−∗
∗
)
(
)
)
(5-22)
对式(5-22)进行 Z 变换得到系统函数和相关函数的 z 变换之间的关系:
同样,对图 5.4 进行 z 变换得
zR
)(
ss
=
2
σ
w
1
zAzA
(
)(
1
−
)
(5-23)
zR
)(
xx
=
2
σ
w
1
zBzB
(
)(
1
−
)
(5-24)
图 5.4 中利用卷积性质还可以找到互相关函数之间的关系:
mR
xs
(
)
=
mnsnxE
()([
+
)]
=
E
[
∑+∞
k
−∞=
knwkb
)(
−
(
)
⋅
mns
(
+
)]
1
= ∑+∞
k
−∞=
kmRkb
)(
+
(
sw
1
)
=
mbmR
sw
1
(
−∗
)
(
)
两边 z 变换得到
zR
)(
xs
=
R
sw
1
zBz
(
)(
1
−
)
(5-25)
如果已知观测信号的自相关函数,求它的 z 变换,然后找到该函数的成对零点、极点,取其
中在单位圆内的那一半零点、极点构成
)(zB
,另外在单位圆外的零、极点构成
1−zB
(
)
,这
样就保证了
)(zB
是因果的,并且是最小相位系统。
从图 5.4 可得
zW
)(1
=
1
zB
)(
zX
)(
(5-26)
由于系统函数
)(zB
的零点和极点都在单位圆内,即是一个物理可实现的最小相位系统,则
1
zB 也是一个物理可实现的最小相移网络函数。我们就可以利用式(5-25)对
)(
)(nx
进行
白化,即把
)(nx
当作输入,
)(1 nw
当作输出,
1
zB 是系统传递函数。
)(
将图 5.1 重新给出,待求的问题就是最小均方误差下的最佳
)(zH
,如图 5.5(a)所示,
为了便于求这个
H opt
)(z
,将图 5.5(a)的滤波器分解成两个级联的滤波器:
1
zB 和 G(z),
)(
如图 5.5(b)所示,则
zGzH
)(
)(
zB
)(
=
(5-27)
nx
)(
=
ns
)(
+
nw
)(
)(nx
1
zB
)(
)(nh
(a)
)(1 nw
(b)
ny
)(
=
)(ˆ
ns
)(zG
ny
)(
=
)(ˆ
ns
图 5.5 利用白化方法求解模型
有了上述的模型后,白化法求解维纳-霍夫方程步骤如下:
1) 对观测信号
)(nx
的自相关函数
(mRxx
)
求 z 变换得到
)(zRxx
2) 利用等式
zR
)(
xx
=
2
σ
w
1
zBzB
(
)(
1
−
)
,找到最小相位系统
)(zB
3) 利用均方误差最小原则求解因果的 G(z)