信号数据有噪声？试试这个Python滤波器

几乎所有真实世界的数据都是有噪声的，数学上可以用下面的简单方程来表示，观测量 = 真实数据 + 噪声 。

所有现实世界的传感器都会受到一定程度的噪声的影响。因此需要使用滤波器从被测量的、有噪声的数据中恢复原始信号。

滤波器的作用是去除信号中不想要的部分，或者提取信号的有用部分。滤波器可以根据信号的频率特性，选择性地通过或者衰减某些频率成分，从而实现对信号的改善或者分析。

滤波器的应用有很多，例如：

• 在通信系统中，滤波器可以抑制谐波、干扰、噪声等，提高信号的质量和可靠性。

• 在音频处理中，滤波器可以调节声音的音色、音量、混响等，增强音乐的效果和表现力。

• 在图像处理中，滤波器可以去除图像的噪点、模糊、锯齿等，提高图像的清晰度和美观度。

• 在数据分析中，滤波器可以消除数据的异常值、趋势、周期等，提取数据的特征和规律。

理论基础

假设我有一组数据集，其中包含连续采样的时间序列数据（如下图所示），每个数据点之间的采样时间（dt）被认为是恒定的。

Savgol滤波器将接收测量数据（以有N个数据点的数组形式），然后输出一个新的数组（也是长度为N），其中有原始数据的过滤/平滑版本。请注意，下面的星号表示输出是对基础信号的估计。

该滤波器通过循环遍历测量数据数组中的每个单独数据点，并生成相应的已过滤输出数据点来工作。下面的图表阐明了这一过程。对于任何给定的数据点（例如，数据点k），在该数据点周围选择一个窗口/子集（用红色括号表示）的数据。然后，在此窗口内拟合多项式，并通过在时间tk（输入数据点k的时间）处评估最佳拟合多项式来计算相应的已过滤输出数据点。然后，将此过程重复所有其他数据点上。

这样可以通过反复曲线拟合和评估多项式来平滑数据。

代码实现

这里假设数据为带有偏移量的阻尼正弦波，该函数可以表示很多物理问题，例如弹簧-阻尼系统、RLC电路等。

数据生成：

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter

np.random.seed(1)

plt.style.use('default')
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 12

fs = 200 
t = np.arange(0, 5, 1/fs) 

#定义信号参数
A = 20 
w = 2*np.pi*1
tau = 2
B = -3 

#原始数据及其一阶导数
y_signal = A*np.exp(-t/tau)*np.cos(w*t) + B
y_deriv = -A*np.exp(-t/tau)*((1/tau)*np.cos(w*t) + w*np.sin(w*t))

#在信号中添加噪声
y_noise = np.random.normal(0,1.5,size=t.shape)
y_measured = y_signal + y_noise

绘制原始数据图像：

ylabel = 'y'
figsize=(8, 4)
fig, ax = plt.subplots(figsize=figsize)
ax.plot(t,y_signal ,color='C2')
ax.set_xlabel('Time [s]')
ax.set_ylabel(ylabel)
fig.tight_layout()

绘制带有噪声的数据图像：

fig, ax = plt.subplots(figsize=figsize)
ax.plot(t,y_measured,ls='', marker='x', label='Noisy Data',color='C2')
ax.set_xlabel('Time [s]')
ax.set_ylabel(ylabel)
fig.tight_layout()

调用savgol_filter滤波器并输入测量的数据集、窗口长度和适合的多项式的阶数。窗口长度是在每个步骤中用来拟合多项式的数据点的数量。

注意，window_length必须是奇数，因为Savgol过滤器通过将窗口围绕特定数据点居中来工作，即给定数据点的左右必须有相同数量的数据点。

运行下述代码对数据进行滤波，并将其与原始信号和噪声测量结果进行对比

#进行滤波
y_filtered = savgol_filter(y_measured, window_length=51, polyorder=2)
#绘图
fig, ax = plt.subplots(figsize=figsize)
ax.plot(t,y_measured,ls='', marker='x', label='Noisy Data',color='C2')
ax.plot(t,y_signal, label='Original Signal',lw=6, ls='--')
ax.plot(t,y_filtered, lw=3,label='Filtered Signal',color='C1')
ax.set_xlabel('Time [s]')
ax.set_ylabel(ylabel)
ax.legend(frameon=True, facecolor='white', fancybox=True)
fig.tight_layout()

最后，我们可以计算噪声数据的导数。这需要对savgol_filter函数进行第二次调用，但是这一次我们还必须传入额外的参数sderiv=1，该参数指定我们要查找的内容和delta=1/fs，它指定采样之间的时间间隔。

y_filtered_deriv = savgol_filter(y_measured, window_length=51,polyorder=2, deriv=1, delta=1/fs)
#绘图
fig, ax = plt.subplots(figsize=figsize)
ax.plot(t,y_deriv, label='Actual',lw=3, ls='--')
ax.plot(t,y_filtered_deriv,ls='-', marker='', label='Estimate',color='C2')
ax.set_xlabel('Time [s]')
ax.set_ylabel(ylabel)
ax.legend(frameon=True, facecolor='white', fancybox=True)
fig.tight_layout()

当对数据进行数值微分时，噪声往往会被放大，通过下图可以看到savgol滤波器面对微分后的数据依旧十分准确。

—— end ——