代码之家 › 专栏 › 技术社区 › aero-man

如何用Python压平基线上下跳跃的数字信号?

ekg signal-processing scipy python

aero-man · 技术社区 · 3 年前

我正在分析SciPy的心电图(EKG)内置数据集,其中一部分如下所示:

上面数据的一个问题是心电图的基线上下跳动了很多。如果你不熟悉心电图或心跳分析,它们应该是平坦的,有一些“QRS波群”(也称为实际心跳)的尖峰,如下所示:

生产心电图监护仪的医疗设备公司使用一些平坦和过滤功能,使心电图更平滑、更容易读取,因为人体的自然运动是不可避免的,会导致心电图信号如上图所示四处跳动。但我不知道他们使用了哪些过滤器/函数。

我如何使用SciPy或编写自定义函数来压平上面的SciPy EKG数据集?

我尝试过的

我试着阅读 SciPy documentation for signal processing ,但尚未找到任何可以“压平”数据或将其带到基线的函数。我是数字信号处理的新手,想知道是否有更好或官方的方法来做到这一点。

我试过通过移动平均值来上下移动这些值,但不可能简单地对Y值进行加法或减法是正确的方法。这太离谱了。

谢谢你的帮助。

问题澄清:修复“基线漂移”

评论中的另一位用户(Christoph Rackwitz)建议使用高通滤波器。在谷歌上搜索高通滤波器时,我发现了一张来自 research paper on ECGs :

他们说第一张照片有“基线漂移”,这就是我真正想回答的问题。如何修复基线漂移?

两篇有用的研究论文

考虑到我在上一节中提到的研究论文,以及 another one I just found about fixing baseline wandering in ECGs ,我会读这些,看看我能找到什么。

0 回复 | 直到 3 年前

Christoph Rackwitz 3 年前

下面展示一些低通滤波器。

你会想了解“因果”滤波器与非因果滤波器,以及FIR和IIR之间的区别。

加载数据:

signal = scipy.datasets.electrocardiogram()
fs = 360 # say the docs
time = np.arange(signal.size) / fs # for plotting only

探索信号:

fig, axs = plt.subplots(3, 1, figsize=(15, 15))
axs[0].plot(time[30000:31000], signal[30000:31000])
axs[1].plot(time[30000:40000], signal[30000:40000])
axs[2].plot(time, signal)
axs[0].set_xlabel('Time (s)')
axs[1].set_xlabel('Time (s)')
axs[2].set_xlabel('Time (s)')
plt.show()

正在尝试一些筛选器:

# Butterworth, first order, 0.5 Hz cutoff
lowpass = scipy.signal.butter(1, 0.5, btype='lowpass', fs=fs, output='sos')
lowpassed = scipy.signal.sosfilt(lowpass, signal)
highpassed = signal - lowpassed

fig, axs = plt.subplots(2, 1, figsize=(15, 10))
axs[0].plot(time[30000:32000], signal[30000:32000])
axs[0].plot(time[30000:32000], lowpassed[30000:32000])
axs[1].plot(time[30000:32000], highpassed[30000:32000])
axs[0].set_xlabel('Time (s)')
axs[1].set_xlabel('Time (s)')
axs[0].set_ylim([-3, +3])
axs[1].set_ylim([-3, +3])
plt.show()

# take note of these coefficients:
>>> scipy.signal.butter(1, 0.5, btype='lowpass', fs=fs, output='ba')
(array([0.00434, 0.00434]), array([ 1.     , -0.99131]))
# and compare to the following...

# Almost the same thing, different formulation: "exponential average"
# y += (x - y) * alpha  # applied to each value X of the signal to produce a new Y

alpha = 1/100 # depends on fs and desired cutoff frequency
lowpassed = scipy.signal.lfilter([alpha], [1, -(1-alpha)], signal)
highpassed = signal - lowpassed

fig, axs = plt.subplots(2, 1, figsize=(15, 10))
axs[0].plot(time[30000:32000], signal[30000:32000])
axs[0].plot(time[30000:32000], lowpassed[30000:32000])
axs[1].plot(time[30000:32000], highpassed[30000:32000])
axs[0].set_xlabel('Time (s)')
axs[1].set_xlabel('Time (s)')
axs[0].set_ylim([-3, +3])
axs[1].set_ylim([-3, +3])
plt.show()

# the first two filters were "causal", i.e. only using past samples.
# downside: lag, phase shift, i.e. the lowpass doesn't quite match/track the signal.
# "non-causal" filters can use future samples.
# this allows to remove the phase shift but the processing introduces a delay instead.
# this delay is irrelevant for offline processing or if it's considered "small enough".
# the following are non-causal.

# median filter. interpret peaks as outliers, so this reveals whatever can be considered "baseline".
# can be causal if the kernel window only covers the past but that introduces lag (noticeable when the signal drifts actively).
# might need another pass of smoothing, on the median filter, before subtracting.
# median filtering CAN be cheap, if using the right data structure. scipy implementation seems less smart, takes noticeable time.

lowpassed = scipy.signal.medfilt(signal, kernel_size=fs+1)
highpassed = signal - lowpassed

fig, axs = plt.subplots(2, 1, figsize=(15, 10))
axs[0].plot(time[30000:32000], signal[30000:32000])
axs[0].plot(time[30000:32000], lowpassed[30000:32000])
axs[1].plot(time[30000:32000], highpassed[30000:32000])
axs[0].set_xlabel('Time (s)')
axs[1].set_xlabel('Time (s)')
axs[0].set_ylim([-3, +3])
axs[1].set_ylim([-3, +3])
plt.show()

lowpassed = scipy.ndimage.gaussian_filter1d(signal, sigma=0.2 * fs, order=0)
highpassed = signal - lowpassed

fig, axs = plt.subplots(2, 1, figsize=(15, 10))
axs[0].plot(time[30000:32000], signal[30000:32000])
axs[0].plot(time[30000:32000], lowpassed[30000:32000])
axs[1].plot(time[30000:32000], highpassed[30000:32000])
axs[0].set_xlabel('Time (s)')
axs[1].set_xlabel('Time (s)')
axs[0].set_ylim([-3, +3])
axs[1].set_ylim([-3, +3])
plt.show()

dmedine 3 年前

为了补充Christoph Rackwitz相当广泛的答案(我从中摘录了绘图代码),您还可以使用非因果(相位保持)IIR高通滤波器来消除基线漂移。这实际上与使用低通来计算移动基线,然后从原始信号中减去它相同:

from scipy import signal, datasets
import matplotlib.pyplot as plt
import numpy as np

ecg = datasets.electrocardiogram()
fs = 360 
time = np.arange(ecg.size) / fs 
[b,a] = signal.butter(4, 0.5, btype='highpass', fs=fs)
filtered = signal.filtfilt(b, a, ecg)

fig, axs = plt.subplots(2, 1, figsize=(15, 10))
axs[0].plot(time[30000:32000], ecg[30000:32000])
axs[1].plot(time[30000:32000], filtered[30000:32000])
axs[0].set_xlabel('Time (s)')
axs[1].set_xlabel('Time (s)')
axs[0].set_ylim([-3, +3])
axs[1].set_ylim([-3, +3])
plt.show()

一定要使用 filtfilt 而不是 filt (或类似 sosfilt )或者IIR滤波器(例如巴特沃斯滤波器)的非线性相位响应将改变信号的形状,并可能破坏进一步的QRS分析。FIR滤波器不会做到这一点,但它也会将原始信号时移与滤波器系数一样多的样本——这将是许多(多达400+),以获得良好的紧密频率响应。