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), 即

  1. (E[y[k]] = E[y[0]]);
  2. (R_{yy}(i, j) = R_{yy}[j-i]).

维纳滤波可以应用于很多场景, 但是这里只讨论下面的去噪的情形:

[x[n] = y[n] + v[n], ]

且假设(y, v)之间相互独立, (E[v]=0).

我们的目标是找到一个滤波(h), 得到一个估计

[hat{y}[n] = hstar x [n], ]

使得下式最小

[E[e^2[n]]. ]

维纳滤波需要分情况讨论, 这里只关注离散的情形, 包括

  1. non-causal: (hat{y}[n] = sum_{k=-infty}^{+infty} h[k]y[n-k]);
  2. FIR: (hat{y}[n] = sum_{k=0}^{N-1} h[k]y[n-k]).

causal的情况这里就不写了.

滤波的推导

自然地, 寻找驻点:

[egin{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(hstar R_{xx}[m] -R_{yx}[m]) \ &=0. end{array} ]

于是必须满足

[ ag{1} h star R_{xx}[m] = R_{yx}[m]. ]

相应的在频率域内存在(假设DFT存在, 参考文献用的z变换, 这个不是特别了解):

[ ag{2} H[u] S_{xx}[u] = S_{yx}[u]. ]

在FIR情况下, 可以通过(1)推导出一个线性方程组从而求解, non-causal下可用(2)推导出结果.

注: (H[u])用了[]是为了保持一致, 在non-causal情况下(H(z))可能更加妥当.

故最优解为:

[ ag{*} 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}. ]

所以:

[ ag{3} H = frac{S_{yy}}{S_{yy} + S_{vv}}. ]

特别的例子

在图像数字处理中, 给出的这样的情形(FIR):

[x[n] = g star y [n] + v[n], ]

(r[n] = g star y [n]),

[egin{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^{-j2pi nu/N}] \ &= G[-u] E[r[n]Y[-u] e^{-j2pi 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^{-j2pi 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])也是类似的.

原文地址:https://www.cnblogs.com/MTandHJ/p/15039061.html