维纳滤波器python 维纳滤波 python

文章目录

• 原理
• scipy调用维纳滤波

原理

从信号与系统的角度出发，有一部分噪声是系统的固有噪声，另一部分，则是对信号的某种响应，换言之，这部分噪声可以理解为一个噪声系统。所以滤除后者，可以理解为去除噪声系统的影响，换言之，就是针对噪声系统做反卷积。

如果噪声系统可以测量，那么反卷积自然可以顺利执行，否则那就要对这个噪声系统进行估计，维纳滤波履行的就是这个思路。

现有一组观测量$\{x_i\}$，由信号和噪声部分组成，即$x_i=y_i+v_i$，$v_i$就是噪声分量。所谓去噪过程，就是找出一组$y_i$的逼近值$\hat y_i$。

所以接下来的问题是，这组量应该怎么找？考虑到求比证难，不妨把这个问题再简化一下，如果找到了一组值，那么如何确定这组值是对的？

正常思路一定是，定义一个误差量$e_i=\hat y_i-y_i$，然后找到这个误差量的某个统计值的极值，考虑到正负号的问题，这个统计量有可能是方差。

而学过信号与系统，就一定会定义一个$h_i$，作为滤波器的冲击响应，则$\hat y_i$可以表示为$h_i$对信号的卷积

$\hat y_i=\sum_k h_kx_{i-k}$

为了便于阅读，可以将上面的式子写为矩阵形式，就是$\hat y=hx$，$e=\hat y-y$。

接下来就可以用$h$来表示误差量$e$的方差

$\begin{aligned} \varepsilon&=E(e^2)=E[(\hat y-y)^2]\\ &=E[(hx-y)^2]\\ &=E[(y^2+hxx^Th^T-2yhx)^2] \end{aligned}$

取极值，那就少不了求导，这个$h$就是要求导的量

$\frac{\partial\varepsilon}{\partial h} =E[2xx^Th^T-2yx]$

极值点处，偏导数为0，从而得到

$h=E[(xx^T)^{-1}]E(xy)$

有了$h$，相当于有了噪声系统的作用形式，也就可以基于此进行反卷积了。

scipy调用维纳滤波

scipy中封装了维纳滤波，其形式为

scipy.signal.wiener(im, mysize=None, noise=None)

其中，im为待滤波数据，mysize为滤波模板的尺寸，noise为系统的噪声，如果为None，那么将自动估计系统噪声。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import wiener
x = np.arange(0,5,0.01)
y = np.sin(x)
yhat = y + np.random.rand(len(yhat))*0.1
yFilt = wiener(yhat, 5) # 维纳滤波
plt.scatter(x, yhat, marker='.')
plt.plot(x, yFilt, c='r')
plt.show()

效果如下