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

Question

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

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

我还尝试使用带有 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了 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你会得到一个更三角形的图形，例如：