Python心电图心跳检测：Pan-Tompkins算法实现与避坑指南

简介这是一套面向生物医学信号处理初学者与算法开发者的Python心电图R峰检测源码合集聚焦心率变异性、心律失常等健康指标计算中的心跳定位问题。资源包共31个文件、约600KB以10个py脚本为核心辅以7个c语言实现、6个csv与1个tsv测试数据以及Makefile、install、license等构建与说明文件覆盖Pan-Tompkins、匹配滤波等多种检测思路并附带MIT-BIH等数据库的基准测试脚本与结果文件。已有2007人学习下载。读者可从中获得完整的算法实现细节、可直接运行的示例数据与评估流程理解预处理、特征提取到R峰判定的全链路并借助HRV时域分析、绘图与统计脚本快速验证效果适合作为可穿戴健康监测、远程监护等场景的算法参考与二次开发起点。1. 从一段 30 秒心电信号说起Python 心跳检测到底在解决什么问题拿到一段 30 秒的单导联心电信号采样率 250 Hz一共 7500 个点。你把它画出来能看到一串规律起伏的波形每个大尖峰就是一次心跳。问题在于人眼能数程序怎么数更麻烦的是真实信号里混着基线漂移、工频干扰、肌电噪声还有幅度忽高忽低的 R 波直接找最大值会漏检也会误检。这就是 Python 心电图心跳检测算法要干的事把原始采样点变成一串准确的心搏位置再算出心率、RR 间期这些指标。它适合做可穿戴设备固件验证、健康类 App 后端、科研信号处理也适合刚学 Python 想找个真实项目练手的人。整套流程不依赖深度学习用 NumPy 和 SciPy 就能跑通核心是滤波、差分、平方、积分、阈值判定这几步。下面按我实际做过的顺序把每一步的参数和坑讲清楚。2. 先把信号洗干净滤波与预处理的参数怎么定2.1 为什么不能直接对原始信号找峰值原始心电信号里R 波幅度大概 1 mV但基线漂移能到几 mV频率在 0.5 Hz 以下工频干扰在 50 Hz 或 60 Hz肌电噪声分布在 20 到 100 Hz。如果直接对原始信号做峰值检测漂移会让阈值失效噪声会制造假峰。所以第一步必须带通滤波把有用频段留下来。心电信号的主要能量集中在 5 到 15 HzQRS 波群的频率成分大约在 10 到 25 Hz。常见做法是设计一个 5 到 15 Hz 的带通滤波器或者更宽一点 0.5 到 40 Hz 保留波形形态。我一般用 5 到 20 Hz因为后面要做差分和平方窄带能让 R 波更突出。import numpy as np from scipy.signal import butter, filtfilt def bandpass_filter(signal, fs, lowcut5.0, highcut20.0, order3): 对心电信号做带通滤波 signal: 一维数组原始采样点 fs: 采样率单位 Hz lowcut/highcut: 通带上下限 order: 滤波器阶数3 阶足够太高会振铃 nyq 0.5 * fs # 奈奎斯特频率 low lowcut / nyq high highcut / nyq b, a butter(order, [low, high], btypeband) # filtfilt 零相位滤波避免波形时移 return filtfilt(b, a, signal)这段代码的关键在filtfilt而不是lfilter。lfilter会引入相位延迟R 波位置会偏移后面算 RR 间期就不准。filtfilt正向反向各滤一次相位抵消代价是计算量翻倍但 30 秒数据无所谓。阶数选 3 是因为阶数越高过渡带越陡但振铃越明显R 波前后会出现虚假波动反而干扰检测。参数上lowcut不要低于 0.5否则漂移没滤干净highcut不要高于 40否则肌电噪声进来。如果采样率只有 100 Hzhighcut要压到 40 以下因为奈奎斯特频率只有 50 Hz。2.2 归一化与去均值让阈值有可比性滤波之后信号幅度还受电极接触、个体差异影响。不同人、不同导联R 波幅度可能差十倍。如果阈值写死一个绝对值换一段数据就翻车。所以要做归一化把信号缩放到统一量纲。def normalize(signal): 去均值后按最大绝对值缩放 signal signal - np.mean(signal) max_abs np.max(np.abs(signal)) if max_abs 0: return signal return signal / max_abs去均值消除直流偏置除以最大绝对值把幅度压到 [-1, 1]。这样后面阈值可以设成相对值比如 0.3、0.5换数据不用改。注意如果信号里有特别大的伪迹最大绝对值会被拉高归一化后 R 波变小。稳妥做法是用 99 百分位代替最大值或者先做一次粗筛。提示归一化前先检查有没有 NaN 或 Inf真实设备采集的数据偶尔会出现直接算会污染整段结果。3. 从波形到心搏Pan-Tompkins 思路的 Python 实现3.1 差分、平方、积分三步为什么有效Pan-Tompkins 是心电 QRS 检测里最经典的算法1985 年提出到现在嵌入式设备还在用。它的逻辑很直白R 波是陡峭的尖峰差分后斜率最大平方后全部变正同时放大高斜率部分再滑动积分把 QRS 波群能量聚成一个包络。这样原本尖锐的 R 波变成一个宽而平滑的峰用阈值就好判断了。def derivative(signal): 五点差分比两点差分更抗噪 return np.ediff1d(signal, to_begin0) def square(signal): return signal ** 2 def moving_window_integration(signal, fs, window_ms150): 滑动窗口积分 window_ms: 窗口长度一般取 QRS 宽度的 1.5 倍左右 window_size int(fs * window_ms / 1000) kernel np.ones(window_size) / window_size return np.convolve(signal, kernel, modesame)差分用np.ediff1d是最简单的但实际我更推荐五点差分公式(-x[n-2] - 2x[n-1] 2x[n1] x[n2]) / 8它对高频噪声抑制更好。平方就是逐点平方没有参数。积分窗口是这里最关键的参数窗口太短包络不平滑一个 QRS 可能出多个峰窗口太长相邻心搏的包络会连在一起快心率时漏检。150 ms 是常用值对应心率 60 到 100 比较稳。如果做运动心率心率可能到 180RR 间期只有 333 ms窗口要缩到 100 ms 左右。3.2 自适应阈值与峰值判定积分后的包络要用阈值找峰。固定阈值不行因为信号幅度会变。自适应阈值根据历史峰值动态调整常见公式是threshold 0.25 * recent_peak 0.75 * running_threshold或者用信号和噪声两个估计值。def detect_peaks(envelope, fs, min_rr_ms200): 在积分包络上检测峰值 min_rr_ms: 最小 RR 间期防止把 T 波当成 R 波 min_distance int(fs * min_rr_ms / 1000) peaks [] threshold 0.3 * np.max(envelope) # 初始阈值 last_peak -min_distance for i in range(1, len(envelope) - 1): if envelope[i] threshold and envelope[i] envelope[i-1] and envelope[i] envelope[i1]: if i - last_peak min_distance: peaks.append(i) last_peak i # 动态更新阈值 threshold 0.25 * envelope[i] 0.75 * threshold return np.array(peaks)min_rr_ms设 200 对应最大心率 300这是生理极限实际用 250 到 300 更稳。阈值更新系数 0.25 和 0.75 是经验值新峰值权重低一点避免一个异常大峰把阈值拉太高导致后续漏检。如果信号质量差可以加一个回溯机制检测到连续漏检时把阈值降 20% 重新扫一遍。3.3 把包络峰映射回原始 R 波位置积分包络的峰和原始 R 波位置有偏移因为滤波、差分、积分都引入了延迟。需要在包络峰附近的一个窗口内回原始信号找绝对值最大的点那才是真正的 R 波。def refine_peaks(raw_signal, envelope_peaks, fs, search_ms100): 在包络峰附近搜索原始信号的最大绝对值点 search_size int(fs * search_ms / 1000) refined [] for p in envelope_peaks: start max(0, p - search_size) end min(len(raw_signal), p search_size) segment raw_signal[start:end] # 找绝对值最大点R 波可能向上也可能向下 local_idx np.argmax(np.abs(segment)) refined.append(start local_idx) return np.array(refined)搜索窗口 100 ms 足够覆盖滤波和积分带来的延迟。注意有些导联 R 波是倒置的所以用绝对值而不是最大值。这一步做完得到的就是最终的心搏位置序列。4. 避坑与排查五个真实踩过的坑4.1 现象心率算出来是实际的两倍原因T 波被当成 R 波检测了。T 波幅度有时能达到 R 波的 30% 到 50%如果阈值偏低或者积分窗口太长把 T 波能量也聚起来就会误检。解决把min_rr_ms从 200 提到 300强制两次检测间隔至少 300 ms。同时检查阈值更新逻辑初始阈值不要低于最大包络的 0.3 倍。如果还不行在检测后加一步规则如果相邻两个峰的间期小于 400 ms 且后一个峰幅度低于前一个的 60%删掉后一个。4.2 现象信号开头一段总是漏检原因自适应阈值初始值依赖np.max(envelope)如果开头恰好没有心搏或者幅度低阈值会偏高。另外filtfilt在信号两端有边缘效应前几十个点不可靠。解决丢弃前 1 秒的数据或者用前 2 秒的包络均值初始化阈值。更稳的做法是先用固定阈值粗检一遍用检出的峰值中位数初始化自适应阈值再正式检测。4.3 现象换了个人数据检测全乱原因归一化用了全局最大值如果新数据里有一个巨大伪迹归一化后 R 波被压得很小阈值相对就高了。解决归一化改用 99 百分位或者分段归一化。也可以在归一化前先做一次简单的伪迹剔除计算信号标准差把超过 5 倍标准差的点截断。4.4 现象程序跑得特别慢原因用了 Python 原生循环逐点处理7500 个点还好如果做 24 小时动态心电几百万个点循环就扛不住了。解决差分、平方、积分全部用 NumPy 向量化。峰值检测的循环可以用scipy.signal.find_peaks替代它底层是 C 实现。如果还要更快用 Numba 的jit装饰器或者把核心逻辑用 Cython 编译。4.5 现象RR 间期算出来有负值或异常大值原因峰值序列没有排序或者检测到了重复峰。另外如果信号中间有段丢失两个峰间隔会异常大。解决检测完先np.sort去重相邻峰间隔小于min_rr_ms的只保留幅度大的那个。RR 间期算完后做中位数滤波把超过中位数 3 倍的间期标记为异常不参与心率计算。5. 验证与进阶用 MIT-BIH 数据检验再谈实时化算法写完不能只看一段数据。公开的 MIT-BIH 心律失常数据库是心电检测的标准测试集里面有专家标注的 R 波位置。你可以下载记录 100、101 这些常见片段用你的算法跑一遍和标注对比。验证指标用敏感度和阳性预测值。敏感度 正确检测数 / 实际心搏数阳性预测值 正确检测数 / 检测总数。好的算法在 MIT-BIH 上敏感度能到 99% 以上。如果低于 95%回去检查滤波频段和阈值逻辑。def evaluate(detected, annotated, tolerance150): 对比检测结果和标注 tolerance: 允许的毫秒误差一般 150 ms fs 360 # MIT-BIH 采样率 tol_samples int(fs * tolerance / 1000) tp 0 matched set() for d in detected: for i, a in enumerate(annotated): if i not in matched and abs(d - a) tol_samples: tp 1 matched.add(i) break fp len(detected) - tp fn len(annotated) - tp sensitivity tp / (tp fn) if (tp fn) 0 else 0 ppv tp / (tp fp) if (tp fp) 0 else 0 return sensitivity, ppv实时化的话把整段处理改成滑动窗口。维护一个 2 秒的缓冲区每次新数据进来对缓冲区做滤波和检测只输出缓冲区中间部分的峰避免边缘效应。积分窗口和阈值更新都要改成流式。这一步做完算法就能跑在树莓派或者手机上。我自己做这个方向最大的教训是不要一上来就调阈值。先把滤波后的波形画出来看确认 R 波清晰、漂移被压住再去调检测参数。很多时候检测不准问题出在滤波频段选错了而不是阈值。另外任何参数都不要写死在代码里做成配置项换数据时先跑一遍可视化心里有数再批量处理。希望帮到你。本文还有配套的精品资源点击获取