Python科学计算：数组💯数据生成💯数据交互💯微积分💯插值💯拟合💯FFT💯卷积

维纳滤波

信号经过系统之后，相当于进行了卷积操作，若想让其复原，只需再用系统进行反卷积即可。如果没有信号，系统却有了响应，那么这种噪声可以理解为系统的噪声。如果系统的数学形式是已知的，这种噪声就很容易滤掉，如果未知，那就需要进行估计，这就是维纳做的工作。

一个有限脉冲响应(finite impulse response, FIR)，其离散形式可通过卷积表示为

$$y_n=\sum _{k=0}^M{b}_k{x}_{n-k}$$

其中，$b_k$的值是未知的，需要通过观测量 { x i } \{x_i\} {xi​}来估计，而观测量包括信号和噪声部分，即$x_i={\hat{x}}_i+v_i$，$v_i$就是噪声分量。相应地，$y_n$也可以分解为信号响应${\hat{y}}_n$和噪声响应，所谓去噪过程，就是找出最佳${\hat{y}}_n$。

有了${\hat{y}}_n$的概念，就可以定义一个误差量$e_i={\hat{y}}_i-y_i$，为了便于计算，将其写为矩阵形式，即$\hat{y}=bx$，$e=\hat{y}-y$，当$e$的方差取极小值时，$\hat{y}$即可达到最佳。这就是维纳滤波的思路。

下面构造一个带有噪声的函数，通过维纳滤波进行处理

$$x=\sin (1.5\pi t(1-t)+2.1)+0.1\sin (2.5\pi t+1)+0.18\cos (7.6\pi t)$$

其滤波效果如下

实现代码如下

import numpy as np
import scipy.signal as ss

t = np.linspace(-1, 1, 201)
PI = 2*np.pi
x = (np.sin(PI*0.75*t*(1-t) + 2.1) +
     0.1*np.sin(PI*1.25*t + 1) +
     0.18*np.cos(PI*3.85*t))

# 原始数据添加噪声
np.random.seed(42)
xn = x + np.random.rand(len(t))

w = ss.wiener(xn, 9) # 维纳滤波

plt.scatter(t, xn, marker='.', label="original")
plt.plot(t, w, c = 'r', label="wiener")
plt.legend()
plt.show()

【wiener】是signal模块中的滤波函数，其输入参数分别是待滤波数据和滤波模板，此外还有一个noise，表示系统噪声，默认为None，表示自行估计噪声。

巴特沃斯滤波器

FIR的特点是无反馈，$y_n$完全由$x_n$决定，如果响应受到反馈的影响，便是无限脉冲响应(infinite impulse response, IIR)，其离散形式变为

$$a_0{y}_n=\sum _{k=0}^M{b}_k{x}_{n-k}-\sum _{k=1}^N{a}_k{y}_{n-k}$$

滤波器设计，就是对$a_k,b_k$具体形式的求解，signal模块中提供了一些函数，对这两种信号进行滤波。仍以函数$x$为例，在添加噪声之后，进行滤波，对于不同的滤波函数，其效果如下

代码为

import scipy.signal as ss
import matplotlib.pyplot as plt

b, a = ss.butter(3, 0.05)
z = ss.lfilter(b, a, xn)
z2 = ss.lfilter(b, a, z)
z3 = ss.filtfilt(b, a, xn)

# 下面为绘图代码
plt.plot(t, z, 'r--', label="lfilter, once")
plt.plot(t, z2, 'g--', label="lfilter, twice")
plt.plot(t, z3, 'b', label="filtfilt")
plt.scatter(t, xn, marker='.', alpha=0.75)

plt.grid()
plt.legend()
plt.show()

其中，

• 【butter】函数生成3阶巴特沃斯滤波器对应的 $a$和$b$值
• 【lfilter】是最基础的脉冲响应滤波器，从左侧开始进行滤波，故而会产生相位差
• 【filtfilt】从正反两个方向滤波，可消除了lfilter产生的相位差

中值滤波

中值滤波，就是挑选出将个滤波模板范围内数据的中位数，例如$[1,3,2,4]$这个数组，给定一个长度为$3$的滤波窗口，那么元素$3$所在位置的滤波范围就是$1,3,2$，其中位数是$2$，所以要把$3$更改为$2$。

import numpy as np
import scipy.signal as ss
x = [1,3,2,4]
ss.medfilt(x,3) # [1, 2, 3, 2]

二维的中值滤波在图像处理中非常常见，对椒盐噪声有着非常霸道的滤除效果。所谓椒盐噪声，如下方左图所示，就是图像中随机产生的黑色和白色的斑点。在使用二维的中值滤波之后，整张图片都变得清澈了。

绘图代码如下。

from scipy.misc import ascent
import matplotlib.pyplot as plt

img = ascent()
img = img[:256, :256]
r = np.random.rand(*img.shape)
img[r>0.96] = 255
img[r<0.04] = 0

plt.subplot(121)
plt.imshow(img, cmap='gray')
plt.axis('off')

plt.subplot(122)
imFilt = ss.medfilt2d(img, [3,3])
plt.imshow(imFilt, cmap='gray')
plt.axis('off')

plt.show()

排序滤波

排序滤波是中值滤波概念的扩充，和中值滤波的区别是，在对滤波窗口中的数据进行排序之后，可以指定用以替代当前数据的数值序号。下面四个矩阵，展示了以$3\times 3$单位矩阵为滤波模板，排序滤波在不同排序参数下的结果。

此滤波过程在scipy中的实现方式如下。

x = np.arange(25).reshape(5, 5).astype(float)
I = np.identity(3)

mats = {"original":x}
for i in range(3):
    mats[f"order_filter:{i}"] = ss.order_filter(x, I, i)

【order_filter】即为signal模块提供的排序滤波函数，以输入参数(x, I, i)为例，表示从矩阵$x$中选出单位阵$I$所覆盖区域中第$i$小的元素。$I$是一个单位阵，就实际情况来看，其覆盖的第一个子阵中，以0为中心，则只能覆盖到2x2的范围，对角元素$0,6$，最小值是0，最大值是6。如以6为中心，则可以完全覆盖3x3的内容，最小值为0，最大值为12。

下面是绘图代码。

def drawMat(x, ax=None):
    M, N = x.shape
    if not ax:
        ax = plt.subplot()
    arrM, arrN = np.arange(M), np.arange(N)
    plt.yticks(arrM+0.5, arrM)
    plt.xticks(arrN+0.5, arrN)
    ax.pcolormesh(x)
    ax.invert_yaxis()
    for i,j in product(range(M),range(N)):
        ax.text(j+0.2, i+0.6, f"{x[i,j]}")

for i,key in enumerate(mats,1):
    ax = plt.subplot(2,2,i)
    drawMat(mats[key], ax)
    plt.title(key)

plt.show()