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

Question

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

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

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

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

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

在那次“开始”（开始时的海拔）之后，我的过滤信号似乎很好，但我不知道为什么会开始。谁能告诉我巴特沃斯过滤器的“开始”是否正常？如果是，有什么方法可以解决吗？

感谢您的帮助。

我的代码：

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，你会得到一个更三角的图形，例如：