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