保姆级教程：用Python从零复现Pan-Tompkins算法（含MIT-BIH数据库验证）
保姆级教程：用Python从零复现Pan-Tompkins算法（含MIT-BIH数据库验证）
在生物医学信号处理领域，心电信号（ECG）分析一直是研究热点。而QRS波群的准确检测，则是整个ECG分析流程中最关键的环节之一。今天，我们将一起动手实现经典的Pan-Tompkins算法，这个诞生于1985年却至今仍在使用的算法瑰宝。
1. 环境准备与数据获取
1.1 Python环境配置
首先确保你的Python环境已安装以下必要库：
pip install numpy matplotlib scipy wfdb pyqtgraph推荐使用Python 3.8+环境，某些库的最新版本可能存在兼容性问题。
1.2 MIT-BIH数据库下载
我们将使用业界标准的MIT-BIH心律失常数据库进行验证：
import wfdb # 下载100号记录（包含正常和异常心跳） record = wfdb.rdrecord('100', pb_dir='mitdb', sampto=3000) ecg_signal = record.p_signal[:,0] # 获取第一导联信号 fs = record.fs # 采样率（通常为360Hz）注意：首次运行会自动下载数据，约3.6MB。完整数据库包含48条记录，但我们的教程先使用单条记录演示。
2. 算法核心模块实现
2.1 带通滤波器设计
Pan-Tompkins算法使用5-15Hz的带通滤波器组合：
from scipy import signal def bandpass_filter(ecg, fs=360): # 低通滤波器（11Hz） b_low = [1, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 1] a_low = [1, -2, 1] # 高通滤波器（5Hz） b_high = [-1/32, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] a_high = [1, -1] # 级联滤波 lowpass = signal.lfilter(b_low, a_low, ecg) return signal.lfilter(b_high, a_high, lowpass)滤波器特性验证技巧：
- 用
signal.freqz()绘制频率响应曲线 - 测试正弦波输入观察衰减效果
- 注意处理边界效应（建议保留200ms过渡区）
2.2 微分与平方处理
五点微分器实现斜率增强：
def derivative(signal): diff = np.zeros_like(signal) for i in range(2, len(signal)-2): diff[i] = (-signal[i-2] - 2*signal[i-1] + 2*signal[i+1] + signal[i+2]) / 8 return diff def squaring(signal): return signal ** 2实际工程中会用卷积优化：
np.convolve(signal, [1, 2, 0, -2, -1]/8, 'same')
2.3 滑动窗口积分
关键参数是窗口宽度（通常150ms）：
def moving_window(signal, window_size=30): integrated = np.zeros_like(signal) for i in range(window_size, len(signal)): integrated[i] = np.sum(signal[i-window_size:i]) / window_size return integrated窗口大小对检测效果的影响：
| 窗口大小 | 优点 | 缺点 |
|---|---|---|
| 100ms | 保留QRS细节 | 可能产生多峰 |
| 150ms | 最佳平衡点 | 略模糊波形 |
| 200ms | 平滑效果好 | 可能合并QRS-T |
2.4 自适应阈值检测
这是算法的核心智能所在：
class AdaptiveThreshold: def __init__(self): self.SPKI = 0 self.NPKI = 0 self.THRESHOLD1 = 0 self.THRESHOLD2 = 0 def update(self, peak, is_qrs): if is_qrs: self.SPKI = 0.125 * peak + 0.875 * self.SPKI else: self.NPKI = 0.125 * peak + 0.875 * self.NPKI self.THRESHOLD1 = self.NPKI + 0.25 * (self.SPKI - self.NPKI) self.THRESHOLD2 = 0.5 * self.THRESHOLD13. 完整流程组装与优化
3.1 主处理流水线
def pan_tompkins(ecg, fs=360): # 1. 滤波 filtered = bandpass_filter(ecg, fs) # 2. 微分+平方 diff = derivative(filtered) squared = squaring(diff) # 3. 积分 integrated = moving_window(squared, window_size=int(0.15*fs)) # 4. 阈值检测 detector = AdaptiveThreshold() peaks = [] for i, value in enumerate(integrated): if value > detector.THRESHOLD1: peaks.append(i) detector.update(value, True) else: detector.update(value, False) return peaks3.2 实时处理技巧
对于连续ECG流处理，需要：
- 维护200ms的缓冲区
- 实现重叠窗口处理
- 阈值状态持久化
- 添加搜索回溯机制
class RealTimeProcessor: def __init__(self, fs=360): self.buffer = np.zeros(int(0.5*fs)) # 500ms缓冲区 self.threshold = AdaptiveThreshold() def process_chunk(self, chunk): # 实现类似pan_tompkins()但带状态保持 pass4. 结果验证与可视化
4.1 性能评估指标
使用MIT-BIH标注数据计算：
def evaluate(peaks, annotations, fs=360): TP = 0 # 正确检测 FP = 0 # 误检 FN = 0 # 漏检 # 允许±100ms的误差窗口 tolerance = int(0.1 * fs) for ann in annotations: found = any(abs(ann - p) < tolerance for p in peaks) if found: TP += 1 else: FN += 1 for p in peaks: if not any(abs(p - ann) < tolerance for ann in annotations): FP += 1 sensitivity = TP / (TP + FN) precision = TP / (TP + FP) return sensitivity, precision4.2 可视化诊断
使用pyqtgraph实现交互式查看：
import pyqtgraph as pg def plot_results(ecg, filtered, integrated, peaks): win = pg.GraphicsLayoutWidget() p1 = win.addPlot(title="原始ECG") p1.plot(ecg, pen='b') p2 = win.addPlot(title="处理后信号") p2.plot(filtered, pen='g') p2.plot(integrated, pen='r') for peak in peaks: p1.addItem(pg.InfiniteLine(pos=peak, angle=90, pen='r')) win.show()典型问题诊断表：
| 现象 | 可能原因 | 解决方案 |
|---|---|---|
| 漏检R波 | 阈值过高 | 降低THRESHOLD1系数 |
| T波误检 | 窗口过宽 | 减小积分窗口至100ms |
| 噪声敏感 | 滤波不足 | 检查滤波器频率响应 |
5. 高级优化方向
5.1 基于心率的动态调整
def dynamic_adjustment(rr_intervals): recent_rr = np.mean(rr_intervals[-8:]) if recent_rr < 0.6: # 心动过速 return 0.8 # 提高灵敏度 elif recent_rr > 1.2: # 心动过缓 return 1.2 # 提高特异性 else: return 1.05.2 多导联融合
def multi_lead_detection(leads): all_peaks = [] for lead in leads: peaks = pan_tompkins(lead) all_peaks.extend(peaks) # 聚类相近的峰值 return cluster_peaks(all_peaks)5.3 现代改进思路
- 机器学习增强：用SVM区分真假峰值
- 小波变换：取代传统滤波环节
- 深度学习：端到端QRS检测
- 硬件加速：用Numba优化计算
@njit def accelerated_filter(ecg): # 使用Numba加速的滤波器实现 pass6. 工程实践建议
实时性优化：
- 预分配数组内存
- 使用循环缓冲区
- 并行化独立计算单元
临床适配技巧：
- 针对运动伪影增加运动传感器融合
- 对不同年龄段设置基础阈值
- 添加起搏脉冲检测模块
调试检查清单：
- [ ] 验证每个处理阶段的信号形态
- [ ] 检查采样率一致性
- [ ] 确认峰值对齐正确
- [ ] 测试不同噪声条件下的鲁棒性
在真实项目中，我们还需要考虑：
- 边缘设备上的内存限制
- 电池供电时的计算功耗
- 不同厂商ECG设备的信号差异
- 临床环境中的电磁干扰特性