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

Stack Overflow用户

提问于 2018-08-27 22:16:22

回答 1查看 3.6K关注 0票数 1

我从我的TI AFE4490上获得了10秒的原始PPG (图像体积图)信号。我的硬件是经过校准的，我使用每秒250个样本来记录这些信号。我最终获得了2500分。

我使用巴特沃斯带通滤波器与lowcut=0.5，highcut=15和order=2。你可以看到我的原始信号和过滤后的信号如下：

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

我在一些文章中读到，对于这种信号，0.5 to和15 to是很好的低频和高频。

在我应用滤波器之前，我使用了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数组中的函数。

signal-processing

digital-filter

butterworth

python

scipy

回答 1

Stack Overflow用户

回答已采纳

发布于 2018-08-28 03:33:40

您在图中看到的初始瞬变是滤波器的阶跃响应，因为突然输入应用于处于静止状态的滤波器。如果您刚刚连接了一台包含此类带通滤波器的物理仪器，仪器的传感器可能会采集到从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.5 by附近的一些低频内容被过滤掉了。如果这就是你所期望的，那就太好了。否则，您仍然可以尝试使用截止频率，以获得您认为产生最佳结果的任何结果。例如(我并不是说这是一个更好的频率范围)，如果您设置lowcut=0.25，您将得到一个更具三角形的图，例如：