Wiener Filtering

文章目录

基本
滤波的推导
- 特别的情况
特别的例子

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]也是类似的.