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

Question

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

上述数据的一个问题是心电图的基线上下跳动很大。如果您不熟悉心电图或心跳分析，它们应该是平坦的，带有一些“QRS 复合体”（又称实际心跳）尖峰，如下所示：

制造心电图监视器的医疗设备公司使用一些平坦化和过滤功能来使心电图更平滑且更易于读取，因为自然的人体运动是不可避免的，并且会导致心电图信号跳跃，如上图所示。但我不知道他们使用哪些过滤器/功能。

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

我尝试过的

我尝试阅读SciPy 信号处理文档，但尚未找到任何可以“展平”数据或使其达到基线的函数。我是数字信号处理方面的新手，想知道是否有更好或官方的方法来做到这一点。

我尝试过将值向上或向下移动一些移动平均值，但 Y 值的简单加法或减法不可能是正确的方法。这太hacky了。

感谢您的帮助。

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

评论中的另一位用户（Christoph Rackwitz）建议使用高通滤波器。在谷歌搜索高通滤波器时，我发现了一篇关于心电图的研究论文中的图像：

他们说第一张图像有“基线漂移”，这就是我真正想用这个问题回答的问题。如何解决基线漂移问题？

两篇有用的研究论文

考虑到我在上一节中提到的研究论文，以及我刚刚发现的另一篇关于修复心电图基线漂移的研究论文，我将阅读这些论文并看看我发现了什么。

Answer 1

这是一些低通滤波器的展示。

您需要了解“因果”滤波器与非因果滤波器，以及 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()