Wiener Filtering
2024-09-06 20:41:21
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}
\]
\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].
\]
h \star R_{xx}[m] = R_{yx}[m].
\]
相应的在频率域内存在(假设DFT存在, 参考文献用的z变换, 这个不是特别了解):
\[\tag{2}
H[u] S_{xx}[u] = S_{yx}[u].
\]
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}.
\]
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},
\]
= \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}}.
\]
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_{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]\)也是类似的.
最新文章
- mongoose数据库连接和操作
- Android—实现自定义相机倒计时拍照
- sys.sysprocesses视图的使用小结
- C语言实现单链表-03版
- 【原】macbook不睡眠的排查与解决
- IE6 BUG 汇总
- BZOJ 1071组队
- Linux 修改hostname 文件
- SPOJ CNTPRIME 13015 Counting Primes (水题,区间更新,求区间的素数个数)
- Error 1937.An error occurred during the installation of assembly...
- WildFly 9.0.2+mod_cluster-1.3.1 集群配置
- BZOJ 1668: [Usaco2006 Oct]Cow Pie Treasures 馅饼里的财富
- 基于FPGA的cordic算法的verilog初步实现
- Rails 执行 rails server 报错 Could not find a JavaScript runtime
- (办公)json报错的解决问题的小经验.
- pycharm安装配置
- Confluence 6 配置手动备份
- maven项目update报错
- Iowait的成因、对系统影响及对策--systemtap
- 转:如何捕获winform程序全局异常?
热门文章
- A Child's History of England.39
- day14搭建博客系统项目
- java9 模块化 jigsaw
- Spark(八)【利用广播小表实现join避免Shuffle】
- [php安全]原生类的利用
- 【编程思想】【设计模式】【创建模式creational】Borg/Monostate
- 【Java多线程】Java 原子操作类API(以AtomicInteger为例)
- 【Java 调优】Java性能优化
- 30个类手写Spring核心原理之MVC映射功能(4)
- 4、Linux下安装tomcat