Signals, Systems and Inference, Chapter 11: Wiener Filtering (mit.edu)
基本
在图像处理的时候, 遇到了这个维纳滤波, 其推导的公式不是很理解, 于是上网查了查, 并做个简单的总结.
| 符号 | 说明 |
|---|---|
| x [ k ] ] x[k]] x[k]] | 观测信号 x x x的第k个元素 |
| y ^ \hat{y} y^ | 为 y y y的一个估计 |
| v v v | 噪声信号 |
| e [ k ] e[k] e[k] | 误差, 为 e [ k ] = y ^ [ k ] − y [ k ] e[k]=\hat{y}[k] - y[k] e[k]=y^[k]−y[k] |
| R x y ( i , j ) R_{xy}(i, j) Rxy(i,j) | 相关系数: E { x [ n − i ] y [ n − j ] } E\{x[n-i]y[n-j]\} E{x[n−i]y[n−j]} |
| S x y S_{xy} Sxy | 傅里叶变换: F { R x y } \mathcal{F} \{R_{xy}\} F{Rxy} |
基本的假设:
y [ k ] y[k] y[k]服从wide-sense stationary (WSS), 即
- E [ y [ k ] ] = E [ y [ 0 ] ] E[y[k]] = E[y[0]] E[y[k]]=E[y[0]];
- R y y ( i , j ) = R y y [ j − i ] R_{yy}(i, j) = R_{yy}[j-i] Ryy(i,j)=Ryy[j−i].
维纳滤波可以应用于很多场景, 但是这里只讨论下面的去噪的情形:
x
[
n
]
=
y
[
n
]
+
v
[
n
]
,
x[n] = y[n] + v[n],
x[n]=y[n]+v[n],
且假设
y
,
v
y, v
y,v之间相互独立,
E
[
v
]
=
0
E[v]=0
E[v]=0.
我们的目标是找到一个滤波
h
h
h, 得到一个估计
y
^
[
n
]
=
h
⋆
x
[
n
]
,
\hat{y}[n] = h\star x [n],
y^[n]=h⋆x[n],
使得下式最小
E
[
e
2
[
n
]
]
.
E[e^2[n]].
E[e2[n]].
维纳滤波需要分情况讨论, 这里只关注离散的情形, 包括
- non-causal: y ^ [ n ] = ∑ k = − ∞ + ∞ h [ k ] y [ n − k ] \hat{y}[n] = \sum_{k=-\infty}^{+\infty} h[k]y[n-k] y^[n]=∑k=−∞+∞h[k]y[n−k];
- FIR: y ^ [ n ] = ∑ k = 0 N − 1 h [ k ] y [ n − k ] \hat{y}[n] = \sum_{k=0}^{N-1} h[k]y[n-k] y^[n]=∑k=0N−1h[k]y[n−k].
causal的情况这里就不写了.
滤波的推导
自然地, 寻找驻点:
∂
E
[
e
2
[
n
]
]
∂
h
[
m
]
=
2
E
[
e
[
n
]
∂
e
[
n
]
∂
h
[
m
]
]
=
2
E
[
e
[
n
]
x
[
n
−
m
]
]
=
2
(
h
⋆
R
x
x
[
m
]
−
R
y
x
[
m
]
)
=
0.
\begin{array}{ll} \frac{\partial E[e^2[n]]}{\partial h[m]} &=2E[e[n]\frac{\partial e[n]}{\partial h[m]}] \\ &=2E[e[n]x[n-m]] \\ &=2(h\star R_{xx}[m] -R_{yx}[m]) \\ &=0. \end{array}
∂h[m]∂E[e2[n]]=2E[e[n]∂h[m]∂e[n]]=2E[e[n]x[n−m]]=2(h⋆Rxx[m]−Ryx[m])=0.
于是必须满足
h
⋆
R
x
x
[
m
]
=
R
y
x
[
m
]
.
(1)
\tag{1} h \star R_{xx}[m] = R_{yx}[m].
h⋆Rxx[m]=Ryx[m].(1)
相应的在频率域内存在(假设DFT存在, 参考文献用的z变换, 这个不是特别了解):
H
[
u
]
S
x
x
[
u
]
=
S
y
x
[
u
]
.
(2)
\tag{2} H[u] S_{xx}[u] = S_{yx}[u].
H[u]Sxx[u]=Syx[u].(2)
在FIR情况下, 可以通过(1)推导出一个线性方程组从而求解, non-causal下可用(2)推导出结果.
注: H [ u ] H[u] H[u]用了[]是为了保持一致, 在non-causal情况下 H ( z ) H(z) H(z)可能更加妥当.
故最优解为:
H
=
S
y
x
/
S
x
x
.
(*)
\tag{*} H = S_{yx} / S_{xx}.
H=Syx/Sxx.(*)
特别的情况
进一步, 由于
S
x
x
=
F
{
R
x
x
}
=
F
{
E
[
(
y
[
n
]
+
v
[
n
]
)
(
y
[
n
−
m
]
+
v
[
n
−
m
]
)
]
}
=
S
y
y
+
S
v
v
,
S_{xx} = \mathcal{F}\{R_{xx}\} = \mathcal{F}\{E[(y[n]+v[n])(y[n-m]+v[n-m])]\} =S_{yy} + S_{vv},
Sxx=F{Rxx}=F{E[(y[n]+v[n])(y[n−m]+v[n−m])]}=Syy+Svv,
最后一步成立的原因是
E
v
=
0
E{v}=0
Ev=0, 且
v
,
y
v, y
v,y相互独立.
同时
S
y
y
=
S
y
x
.
S_{yy} = S_{yx}.
Syy=Syx.
所以:
H
=
S
y
y
S
y
y
+
S
v
v
.
(3)
\tag{3} H = \frac{S_{yy}}{S_{yy} + S_{vv}}.
H=Syy+SvvSyy.(3)
特别的例子
在图像数字处理中, 给出的这样的情形(FIR):
x
[
n
]
=
g
⋆
y
[
n
]
+
v
[
n
]
,
x[n] = g \star y [n] + v[n],
x[n]=g⋆y[n]+v[n],
记
r
[
n
]
=
g
⋆
y
[
n
]
r[n] = g \star y [n]
r[n]=g⋆y[n],
则
S
r
r
[
u
]
=
F
{
R
r
r
[
m
]
}
[
u
]
=
F
{
E
[
r
[
n
]
r
[
n
−
m
]
]
}
[
u
]
=
E
[
r
[
n
]
F
{
r
[
n
−
m
]
}
]
[
u
]
=
E
[
r
[
n
]
G
[
−
u
]
Y
[
−
u
]
e
−
j
2
π
n
u
/
N
]
=
G
[
−
u
]
E
[
r
[
n
]
Y
[
−
u
]
e
−
j
2
π
n
u
/
N
]
=
G
[
−
u
]
F
{
R
r
y
[
m
]
}
[
u
]
=
G
[
−
u
]
F
{
E
[
r
[
n
+
m
]
y
[
n
]
]
}
[
u
]
=
G
[
−
u
]
E
[
F
{
r
[
n
+
m
]
}
[
u
]
y
[
n
]
]
=
G
[
−
u
]
E
[
G
[
u
]
Y
[
u
]
e
−
j
2
π
n
u
/
N
y
[
n
]
]
=
G
[
−
u
]
G
[
u
]
F
{
R
y
y
[
m
]
}
[
u
]
=
G
[
−
u
]
G
[
u
]
S
y
y
[
u
]
\begin{aligned} S_{rr}[u] &= \mathcal{F}\{R_{rr}[m]\}[u] \\ &= \mathcal{F}\{E[r[n] r[n-m]]\}[u] \\ &= E[r[n] \mathcal{F}\{r[n-m]\}][u] \\ &= E[r[n] G[-u] Y[-u] e^{-j2\pi nu/N}] \\ &= G[-u] E[r[n]Y[-u] e^{-j2\pi nu/N}] \\ &= G[-u] \mathcal{F}\{R_{ry}[m]\}[u] \\ &= G[-u] \mathcal{F}\{E[r[n+m]y[n]]\}[u] \\ &= G[-u] E[\mathcal{F}\{r[n+m]\}[u]y[n]] \\ &= G[-u] E[G[u]Y[u]e^{-j2\pi nu/N}y[n]] \\ &= G[-u]G[u] \mathcal{F}\{R_{yy}[m]\}[u] \\ &= G[-u]G[u] S_{yy}[u] \end{aligned}
Srr[u]=F{Rrr[m]}[u]=F{E[r[n]r[n−m]]}[u]=E[r[n]F{r[n−m]}][u]=E[r[n]G[−u]Y[−u]e−j2πnu/N]=G[−u]E[r[n]Y[−u]e−j2πnu/N]=G[−u]F{Rry[m]}[u]=G[−u]F{E[r[n+m]y[n]]}[u]=G[−u]E[F{r[n+m]}[u]y[n]]=G[−u]E[G[u]Y[u]e−j2πnu/Ny[n]]=G[−u]G[u]F{Ryy[m]}[u]=G[−u]G[u]Syy[u]
由于
S
y
x
=
S
y
r
=
G
[
−
u
]
S
y
y
[
u
]
,
S_{yx} = S_{yr} = G[-u]S_{yy}[u],
Syx=Syr=G[−u]Syy[u],
证明和上面是类似的,
S
x
x
=
S
r
r
+
S
v
v
,
S_{xx} = S_{rr} + S_{vv},
Sxx=Srr+Svv,
于是
H
[
u
]
=
1
G
[
u
]
1
1
+
S
v
v
[
u
]
/
(
G
[
−
u
]
G
[
u
]
S
y
y
[
u
]
)
.
H[u] = \frac{1}{G[u]}\frac{1}{1 + S_{vv}[u] / (G[-u]G[u]S_{yy}[u])}.
H[u]=G[u]11+Svv[u]/(G[−u]G[u]Syy[u])1.
当
g
g
g为实的时候, 有
G
[
−
u
]
G
[
u
]
=
∣
G
[
u
]
∣
2
G[-u]G[u]=|G[u]|^2
G[−u]G[u]=∣G[u]∣2, 在数字图像处理书中, 给出的公式中:
S
y
y
[
u
]
=
∣
F
[
u
]
∣
2
,
S
v
v
[
u
]
=
∣
V
[
u
]
∣
2
,
S_{yy}[u] = |F[u]|^2, S_{vv}[u] = |V[u]|^2,
Syy[u]=∣F[u]∣2,Svv[u]=∣V[u]∣2,
个人觉得是这里的期望
E
E
E用
R
y
y
[
m
]
=
∑
n
=
0
N
−
1
y
[
n
]
y
[
n
−
m
]
,
R_{yy}[m] = \sum_{n=0}^{N-1} y[n]y[n-m],
Ryy[m]=n=0∑N−1y[n]y[n−m],
代替了,
所以
F
{
R
y
y
[
m
]
}
=
F
∗
(
u
)
F
(
u
)
=
∣
F
(
u
)
∣
2
,
\mathcal{F}\{R_{yy}[m]\} = F^*(u)F(u) = |F(u)|^2,
F{Ryy[m]}=F∗(u)F(u)=∣F(u)∣2,
这里假设
y
[
n
]
∈
R
y[n] \in \mathbb{R}
y[n]∈R.
S v v [ u ] S_{vv}[u] Svv[u]也是类似的.
5133
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
