无需Scipy在Pan-Tompkins算法中实现ECG带通滤波器

ECG信号的带通滤波核心是高通去除基线漂移（截止0.5-1Hz）和低通滤除高频噪声（截止30-40Hz），下面直接给出纯Python实现方案，分一阶简易版和二阶精准版，全程不依赖Scipy：

一阶IIR滤波（快速实现，适配作业需求） #

一阶IIR是最易手动实现的滤波方式，计算量小，完全满足ECG预处理的基本要求。

1. 一阶高通滤波器（去除基线漂移） #

基于离散化RC高通电路的递归公式实现：

import numpy as np

def iir_highpass(x, cutoff_hz, sample_rate):
    alpha = np.exp(-2 * np.pi * cutoff_hz / sample_rate)
    filtered = np.zeros_like(x)
    # 初始条件：第一个点无前置数据，设为0（高通天然滤除直流分量）
    filtered[0] = 0.0
    for n in range(1, len(x)):
        filtered[n] = alpha * filtered[n-1] + alpha * (x[n] - x[n-1])
    return filtered

2. 一阶低通滤波器（滤除高频噪声） #

基于离散化RC低通电路的递归公式实现：

def iir_lowpass(x, cutoff_hz, sample_rate):
    alpha = np.exp(-2 * np.pi * cutoff_hz / sample_rate)
    filtered = np.zeros_like(x)
    # 初始条件：第一个点直接取原始信号值
    filtered[0] = x[0]
    for n in range(1, len(x)):
        filtered[n] = alpha * filtered[n-1] + (1 - alpha) * x[n]
    return filtered

3. 级联成带通滤波器 #

先通过高通去除基线漂移，再通过低通滤除高频噪声：

def ecg_bandpass(x, high_cutoff=40, low_cutoff=0.5, sample_rate=250):
    x_high = iir_highpass(x, low_cutoff, sample_rate)
    x_band = iir_lowpass(x_high, high_cutoff, sample_rate)
    return x_band

二阶巴特沃斯滤波（更陡峭的滤波特性） #

如果作业要求更精准的滤波效果，可通过双线性变换法手动实现二阶巴特沃斯滤波器，完全避开Scipy的butter等函数：

二阶低通滤波器 #

def butter2_lowpass(x, cutoff_hz, sample_rate):
    # 双线性变换预畸变，避免频率混叠
    fc_warp = 2 * sample_rate * np.tan(np.pi * cutoff_hz / sample_rate)
    # 计算巴特沃斯二阶系数
    a0 = fc_warp**2 + 2 * np.sqrt(2) * fc_warp * sample_rate + (2 * sample_rate)**2
    a1 = 2 * (fc_warp**2 - (2 * sample_rate)**2)
    a2 = fc_warp**2 - 2 * np.sqrt(2) * fc_warp * sample_rate + (2 * sample_rate)**2
    b0 = fc_warp**2
    b1 = 2 * fc_warp**2
    b2 = fc_warp**2

    filtered = np.zeros_like(x)
    # 初始化前两个点的输出
    filtered[0] = (b0 * x[0]) / a0
    filtered[1] = (b0 * x[1] + b1 * x[0] - a1 * filtered[0]) / a0
    # 迭代计算后续点
    for n in range(2, len(x)):
        filtered[n] = (b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*filtered[n-1] - a2*filtered[n-2]) / a0
    return filtered

二阶高通滤波器 #

同理推导二阶巴特沃斯高通的系数，实现代码：

def butter2_highpass(x, cutoff_hz, sample_rate):
    fc_warp = 2 * sample_rate * np.tan(np.pi * cutoff_hz / sample_rate)
    a0 = fc_warp**2 + 2 * np.sqrt(2) * fc_warp * sample_rate + (2 * sample_rate)**2
    a1 = 2 * (fc_warp**2 - (2 * sample_rate)**2)
    a2 = fc_warp**2 - 2 * np.sqrt(2) * fc_warp * sample_rate + (2 * sample_rate)**2
    b0 = (2 * sample_rate)**2
    b1 = -2 * (2 * sample_rate)**2
    b2 = (2 * sample_rate)**2

    filtered = np.zeros_like(x)
    filtered[0] = (b0 * x[0]) / a0
    filtered[1] = (b0 * x[1] + b1 * x[0] - a1 * filtered[0]) / a0
    for n in range(2, len(x)):
        filtered[n] = (b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*filtered[n-1] - a2*filtered[n-2]) / a0
    return filtered

二阶带通滤波器 #

级联二阶高通和低通即可：

def ecg_bandpass_butter2(x, high_cutoff=40, low_cutoff=0.5, sample_rate=250):
    x_high = butter2_highpass(x, low_cutoff, sample_rate)
    x_band = butter2_lowpass(x_high, high_cutoff, sample_rate)
    return x_band

注意事项 #

• 采样率sample_rate需与你的ECG信号匹配（常用250Hz或360Hz）
• 一阶滤波足以应付作业中的QRS检测需求，二阶适合对滤波效果要求更高的场景
• 所有代码仅依赖numpy（作业中通常允许使用numpy处理数组），完全未调用Scipy的滤波相关函数

内容的提问来源于stack exchange，提问作者Armand Surbakti