使用 Python 对原始信号应用合适的巴特沃斯滤波器

Question

我从我的 TI AFE4490 获得了 10 秒的原始 PPG（光体积描记图）信号。我的硬件已校准，我每秒使用 250 个样本来记录这些信号。我最后获得了2500点

我使用了 lowcut=0.5 、 highcut=15 和 order=2 的 Butterworth 带通滤波器。您可以在下面看到我的原始信号和过滤后的信号：

我还尝试使用具有 lowcut=15 和 order=2 的 Butterworth 低通滤波器对此进行过滤。如您所见，我的原始信号和过滤后的信号如下：

我在一些文章中读到 0.5Hz 和 15Hz 是此类信号的良好低切和高切频率。

在应用过滤器之前，我使用了 Scipy Butterworth（来自 scipy 文档）算法来显示过滤器响应，效果很好。

之后我的过滤信号似乎很好"start"（开始时的高度），但我不知道为什么那个开始。谁能告诉我 "start" 在 Butterworth 过滤器中是否正常？如果是，有什么方法可以解决吗？

感谢您的帮助。

我的代码：

RED, IR, nSamples, sRate = getAFESignal()

period = 1/sRate # Signal period.
# Desired cutoff frequency (in Hz) and filter order.
lowcut = 0.5
highcut = 15
orders = 2

plt.figure(1)

x = np.linspace(0, nSamples*period, nSamples, endpoint=True)

plt.subplot(2,1,1)
y = IR
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.plot(x,y, label='Noisy signal')


plt.subplot(2,1,2)
yf = butter_bandpass_filter(IR, lowcut, highcut, nSamples, order=orders)
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.plot(x, yf, label='Filtered signal')

plt.grid()
plt.show()

函数getAFEsignal()只是一个读取.txt文件并将所有内容放入两个numpy数组的函数。

Answer 1

您在图表中看到的初始瞬态是滤波器的阶跃响应，因为突然输入应用于处于静止状态的滤波器。如果您刚刚连接了一个包含此类带通滤波器的物理仪器，则该仪器的传感器可能会拾取从 0（当探头断开连接时）跳到第一个样本值 ~0.126V 的输入数据样本。仪器滤波器的响应会显示出类似的瞬变。

然而，您可能更感兴趣的是仪器在不再受到这些外部因素（例如连接的探头）的干扰后，有时间适应以下特性后的稳态响应感兴趣的信号。

实现此目的的一种方法是使用足够长的数据样本并丢弃初始瞬态。另一种方法是强制滤波器的初始内部状态接近于在您的第一个输入样本之前应用类似幅度的信号一段时间时可能预期的状态。例如，这可以通过使用 scipy.signal.lfilter_zi.

设置初始条件来完成

现在，我假设您已经使用了 butter_bandpass_filter from SciPy Cookbook，它不考虑过滤器的初始条件。幸运的是，可以很容易地为此修改它：

from scipy.signal import butter, lfilter, lfilter_zi

def butter_bandpass_filter_zi(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    zi = lfilter_zi(b, a)
    y,zo = lfilter(b, a, data, zi=zi*data[0])
    return y

此时要注意的另一件事是，您将 butter_bandpass_filter 调用为：

yf = butter_bandpass_filter(IR, lowcut, highcut, nSamples, order=orders)

传递 nSamples（样本总数，在您的情况下为 2500）作为第四个参数，而函数期望采样率（在您的情况下为 250）。两个量之间的因数 10 的效果相当于将滤波范围从 [0.5,15] Hz 减小到 [0.05,1.5] Hz。要获得预期的带通频率范围，您应该将 sRate 作为第四个参数传递：

yf = butter_bandpass_filter_zi(IR, lowcut, highcut, sRate, order=orders)

最后，您可能会注意到最后一个输出比输入的三角形少一点。这是由于滤掉了 0.5Hz 附近的一些低频内容造成的。如果那是您所期待的，那就太好了。否则，您仍然可以尝试使用截止频率来获得您认为会产生最佳结果的任何频率。例如（我并不是说这是一个更好的频率范围），如果你设置 lowcut=0.25 你会得到一个更三角形的图形，例如：

使用 Python 对原始信号应用合适的巴特沃斯滤波器

Applying an suitable butterworth filter on raw signal using Python

python

scipy

signal-processing

digital-filter

butterworth