使用Python对原始信号应用合适的butterworth滤波器

2024-07-05 14:54:37 发布

您现在位置:Python中文网/ 问答频道 /正文

我从我的TI AFE4490获得了一个10秒的原始PPG(光容积图)信号。我的硬件经过校准,我用每秒250个样本来记录这些信号。最后我得了2500分。在

我使用巴特沃斯带通滤波器,低切=0.5,高截=15,阶数=2。你可以看到我的原始和过滤信号如下:

Top = Raw Photoplethysmogram signal, Bottom = Filtered Photoplethysmogram signal

我也试着用巴特沃斯低通滤波器,低切=15,阶数=2。如您所见,我的原始和过滤信号如下:

Top = Raw Photoplethysmogram signal, Bottom = Filtered Photoplethysmogram signal

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

在我应用过滤器之前,我使用了一个scipybartworth(来自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数组的函数。在


Tags: 过滤器ir信号orderpltfilterperiodorders
1条回答
网友
1楼 · 发布于 2024-07-05 14:54:37

您在图中看到的初始瞬态是滤波器的阶跃响应,当一个突然的输入应用于处于静止状态的滤波器时。如果你刚刚连接了一个包括这样一个带通滤波器的物理仪器,仪器的传感器可能已经从0(探头断开时)跳到第一个采样值~0.126V的输入数据样本。然后仪器滤波器的响应将显示类似的瞬态。在

然而,在仪器不再受到这些外部因素(例如连接的探头)的干扰,并且有时间适应感兴趣信号的特性之后,您可能对仪器的稳态响应更感兴趣。在

实现这一点的一种方法是使用足够长的数据样本并丢弃初始瞬态。另一种方法是强制滤波器的初始内部状态接近于在第一个输入样本之前应用了一段时间的类似振幅的信号。这可以通过使用^{}设置初始条件来完成。在

现在,我假设您使用了^{} 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

{a3}

此时需要注意的另一件事是,您将butter_bandpass_filter调用为:

^{pr2}$

传递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)

enter image description here

最后,您可能会注意到最后的输出比输入的三角形要小一些。这是由于一些接近0.5Hz的低频内容被过滤掉。如果这是你所期待的,那就太好了。否则,你仍然可以玩的截止频率,以获得任何你觉得最好的结果。例如(我并不是说这是一个更好的频率范围),如果你设置lowcut=0.25,你会得到一个更三角形的图,比如:

enter image description here

相关问题 更多 >