简介这是一套面向生物医学信号处理开发者与算法学习者的Python实现心电图R峰检测算法源码包重点针对MIT-BIH等心电数据库的心跳定位、HRV计算与心律失常指标评估适合有Python基础并希望深入ECG分析的读者。包内共31个文件约600KB包含10个Python脚本和7个C程序Python侧覆盖Pan-Tompkins、Jafarnejadsani等经典检测器及测试、绘图流程C程序负责滤波与心率统计另有CSV/TSV结果文件和不同采样率模板可直接复现基准测试。整个工程自带安装脚本、Makefile、示例数据和运行入口从环境配置到结果输出层次清晰能帮助读者快速跑通检测流程并对照源码理解阈值自适应、峰谷判定与噪声抑制等细节。已有2007人学习适合正在开展可穿戴监测、远程心电预警或生物信号课程设计的开发者作为参考。1. 拿到心电图数据之后第一件事为什么是找 R 波手里有一份心电图数据几万甚至几十万个采样点真正能用的信息其实只有一小段每一次心跳对应一个 QRS 波群其中最陡、最高、最醒目的那个尖峰叫 R 波。所谓 Python 心电图心跳检测算法核心就是把这个尖峰的位置找出来剩下的心率计算、心律不齐分析、疲劳监测全都建立在 R 波定位之上。很多人一开始就想着上深度学习但我个人建议先把经典算法跑通因为这类算法不依赖标注数据也能调试参数透明出了问题一眼能看出是滤波还是阈值的问题非常适合低成本落地和做技术预研。这篇文章要讲的是我从零搭起一组用 Python 实现的心跳检测流程先做带通滤波把噪声压下去再用斜率、幅度和滑动窗口积分把 R 波从波形里抠出来最后用公开数据库核对精度并整理出我在实际数据上踩过的几个坑。如果你刚装好 Python 和 NumPy 还没想好从哪个方向入手或者已经用阈值法做过检测但总觉得误检很多这篇应该能给你一条可以直接照着走下去的路。2. 预处理做对一半为什么先带通滤波再检测而不是直接找峰值2.1 心电信号里的噪声都藏在哪些频段直接对原始心电信号做峰值检测结果通常惨不忍睹。原因很简单ECG 是一种低频信号但采集过程中会叠加上各种噪声。基线漂移来自电极移动和呼吸频率在 0.5 Hz 以下肌电干扰来自肌肉紧张频率在 20 Hz 以上工频干扰在 50 或 60 Hz幅度可能比 QRS 波本身还大。而 QRS 波群的主要能量集中在 5 到 20 Hz 之间。这一条频带信息就是整个预处理设计的出发点。常见做法是先设计一个带通滤波器只保留 5 到 20 Hz 这个频段把低频漂移和高频肌电一起切掉。我用的是 Butterworth 滤波器因为它的通带内最平坦参数少且行为可预期。很多人第一次写滤波代码时直接调用 scipy.signal.butter 拿到系数就上 lfilter结果发现波形整体往后平移了R 波位置全部偏移几十毫秒。这是因为普通滤波器都有相位延迟而心跳检测对位置敏感所以我这里用 filtfilt 做零相位滤波前向和反向各过一次相位延迟相互抵消。2.2 用 scipy 搭一版零相位带通滤波下面这段代码是我常用的预处理函数输入原始信号和采样率输出滤波后的信号。import numpy as np from scipy.signal import butter, filtfilt def bandpass_ecg(signal, fs, lowcut5.0, highcut20.0, order2): 对心电信号做零相位带通滤波。 signal: 一维数组原始 ECG 波形 fs: 采样率单位 Hz lowcut/highcut: 通带上下限 order: 滤波器阶数2 阶通常够用 nyquist 0.5 * fs low lowcut / nyquist high highcut / nyquist # 注意 butter 返回的是传递函数的分子分母系数 b, a butter(order, [low, high], btypeband) # filtfilt 做零相位滤波R 波位置不会漂移 filtered filtfilt(b, a, signal) return filtered逻辑说明butter 函数把频率参数归一化到 Nyquist 频率之后返回 IIR 滤波器的分子分母系数。filtfilt 先正向滤波一次再把结果反转后反向滤波一次抵消相位偏移代价是计算量翻倍但对于离线处理的数据完全可接受。order 取 2 是因为更高阶的滤波器虽然过渡带更陡但数值稳定性和波形畸变风险也会增加对心跳检测来说 2 到 4 阶足够。参数说明lowcut 建议设在 3 到 5 Hz太高会把 QRS 波的低频成分切掉导致波形变尖太低则挡不住基线漂移。highcut 建议设在 15 到 25 HzT 波和 P 波频率较低设到 20 Hz 时它们会被明显抑制对 R 波影响不大。fs 必须跟采集设备一致200 Hz 和 1000 Hz 的归一化计算完全不同这一步错了后面全是错的。2.3 滤波前后到底该看什么怎么确认滤波参数是对的我的习惯是把原始信号和滤波后信号画在同一张图上对比横坐标按秒而不是按采样点。画图的时候注意一个坑采样点太多横坐标会挤成一团根本看不出节奏所以要自己换算成秒再画。判断标准有三个。第一滤波后的基线应该基本是一条平线没有大幅上下起伏第二QRS 波群仍然保留明显的尖峰形态没有被磨圆第三T 波幅度应该明显小于 R 波如果 T 波还跟 R 波一样高说明高通截止频率设得太低了。3. 核心检测逻辑从滤波信号里把 R 波的位置抠出来3.1 斜率、平方、滑动窗口积分把尖峰变成一块突出的区域滤波之后的信号仍然是一个连续波形直接在原始幅度上做阈值检测会遇到一个问题QRS 波的幅度在不同导联、不同人身上差异很大有的 R 波 0.5 mV有的高达 2 mV一个固定阈值没法通用。经典 Pan-Tompkins 算法给出了一套不依赖绝对幅度的处理流程核心分三步求导、平方、滑动窗口积分。求导是为了突出斜率R 波的上升沿是心电信号中变化最快的地方求导之后在对应位置会产生一个明显的正峰值而 T 波和 P 波因为变化缓慢导数很小。平方操作把所有负值变成正值同时进一步放大高频成分的差异。最后做滑动窗口积分把一段短时间内的能量累加起来QRS 波群所在的位置会形成一个凸起的平台这个平台远比原始尖峰容易检测。下面这段代码把三步连在一起实现注意移动窗口积分我用了累计和做差分避免反复求和导致计算量爆炸。def qrs_energy_transform(filtered, fs, win_sec0.15): 对滤波后的信号做差分、平方、滑动窗口积分。 返回变换后的能量曲线和窗口长度。 # 一阶差分近似求导注意长度减少 1 diff np.diff(filtered) # 平方突出高频成分转为正值 squared diff ** 2 # 窗口长度按采样率换算150ms 是常见值 win_len max(1, int(win_sec * fs)) # 用 cumsum 高效算滑动窗口和 cum np.cumsum(np.concatenate(([0.0], squared))) integrated np.zeros_like(squared) for i in range(len(squared) - win_len 1): integrated[i] cum[i win_len] - cum[i] return integrated, win_len逻辑说明np.diff 得到的导数序列长度比原信号少 1所以后面所有操作的索引都会整体偏移一点定位 R 波时要补偿回来。cumsum 的用法是关键它先算累计和然后窗口和用两次索引相减得到这里看起来是 for 循环但实际只访问 O(n) 次比每次 sum 一个窗口快非常多。参数说明win_sec 是滑动窗口的物理时长我这里取 150 ms。这个值的物理意义是 QRS 波群本身的持续时间约 80 到 120 ms窗口太小会把一个 QRS 波切成两段产生双峰窗口太大又会让两个相邻的心跳融合在一起150 ms 是经验上最稳的选择。如果采样率是 250 Hzwin_len 就是 38 个采样点如果是 1000 Hzwin_len 就是 150 个采样点务必按 fs 换算而不是用固定点数。3.2 动态阈值和不应期两个最容易翻车的设计点能量变换之后R 波位置会变成一处处明显的凸起但这还不够你仍然需要回答一个问题多高的凸起算一次心跳如果设固定阈值信号幅值漂移时就会失效。我常用自适应动态阈值思路是维护两个量信号峰值水平 signal_peak 和噪声水平 noise_peak。每次检测到新峰值就按比例更新这两个水平阈值取两者的加权值。此外还有一个必须处理的问题R 波本身在积分后有一定宽度如果不做抑制同一个心跳会被连报好几次。解决办法是不应期检测到一次 R 波后在其后 200 ms 内不再接受新的检测同时对这个时间段内的最大峰值做一次比较保留最高峰作为真实位置。这就是把精度从“次”提升到“位置”的关键一步。def detect_r_peaks(integrated, fs, win_len, threshold_ratio0.4, refactory_sec0.2): 在能量曲线上做动态阈值检测返回 R 波位置索引。 integrated: 上一步滑动窗口积分的输出 fs: 采样率 win_len: 滑动窗口长度用于补偿差分导致的位置偏移 peaks [] signal_peak 0.0 noise_peak 0.0 refractory_samples int(refactory_sec * fs) last_peak_idx -refractory_samples # 初始化为足够远 for i in range(len(integrated)): val integrated[i] # 更新噪声/信号峰值估计 if val signal_peak: signal_peak 0.8 * signal_peak 0.2 * val else: noise_peak 0.8 * noise_peak 0.2 * val threshold noise_peak threshold_ratio * (signal_peak - noise_peak) if val threshold and (i - last_peak_idx) refractory_samples: # 进入不应期在这个时间窗内找局部最大值 search_start i search_end min(len(integrated), i refractory_samples) local_max search_start int(np.argmax(integrated[search_start:search_end])) peaks.append(local_max) last_peak_idx local_max # 把搜索范围跳过去直接以局部峰为下一次起点 i local_max # 差分导致波形整体往前偏移补偿回来得到原始信号上的位置 peaks [p win_len // 2 1 for p in peaks] return peaks逻辑说明signal_peak 和 noise_peak 都用了指数滑动平均0.8 和 0.2 的权重决定了跟踪速度。权重越大跟踪越慢抗干扰能力越强但对幅值突变反应迟钝。阈值本身是噪声水平加上信号与噪声差值的 40%这个比例意味信号峰只要比噪声高出一定比例就会被接受不会因为整体幅度变大变小而失效。参数说明threshold_ratio 的常见取值在 0.3 到 0.5 之间取 0.4 是比较均衡的起点。refactory_sec 设 200 ms 对应正常人最大心率约 300 bpm超过 300 的检测结果应当视为噪声。如果你处理的是运动员静息数据200 ms 完全没有问题如果是运动负荷测试数据心率可能到 180 bpm对应 R-R 间期约 333 ms200 ms 依然安全。3.3 这一版算法的默认参数一览参数默认值调节方向说明lowcut5 Hz基线漂移严重时抬高到 8 Hz太高会削掉 QRS 低频成分highcut20 Hz肌电干扰严重时降到 15 Hz太低会让 T 波变大filter order2波形仍毛糙时加到 34 阶以上容易数值不稳win_sec150 ms双峰时加大融合时减小必须跟随采样率换算threshold_ratio0.4误检多时调高漏检多时调低这个值是全局灵敏度旋钮refactory_sec200 ms心率极高的数据调小到 120 ms不要小于 100 ms4. 检测完怎么算对用公开数据库核对检测精度4.1 有标注数据的评估MIT-BIH 数据库的读取思路算法写完不能只看两张图感觉差不多就行要用带标注的数据算准确率。最常见的公开数据集是 MIT-BIH Arrhythmia Database每个记录包含两路心电信号和一组由专家标注的心跳位置。标注信息按 beat 类型分类N 代表正常心跳V 代表室性早搏但心跳检测这一步不用区分类型只需要把标注里的绝对位置读出来即可。这个数据库的原始格式不是常见的那种 CSV而是 .dat、.hea、.atr 三个文件配合使用需要从 PhysioNet 提供的工具或现成解析库读取。不过你现在只需要明确一点不管格式怎么解析最终你手里应该有两组数字一组是算法报告的 R 波位置索引另一组是专家标注的 R 波位置索引接下来做的就是两组数字的匹配。4.2 写一个评估函数TP、FP、FN 怎么算匹配规则我通常用容差窗口算法检测位置与某个标注位置之间距离小于 150 ms就认为是一次真阳性 TP。一个标注位置只能匹配一次一个检测位置也只能匹配一次这样算出的指标才严格。def evaluate_detection(detections, annotations, fs, tolerance_sec0.15): 计算检测结果的 TP、FP、FN 和灵敏度、阳性预测值。 detections: 算法检测到的 R 波位置索引列表 annotations: 专家标注的 R 波位置索引列表 fs: 采样率 tolerance int(tolerance_sec * fs) detections sorted(detections) annotations sorted(annotations) tp 0 matched_det set() matched_ann set() for d_idx, d in enumerate(detections): best_ann None best_dist float(inf) for a_idx, a in enumerate(annotations): if a_idx in matched_ann: continue dist abs(d - a) if dist tolerance and dist best_dist: best_dist dist best_ann a_idx if best_ann is not None: tp 1 matched_det.add(d_idx) matched_ann.add(best_ann) fp len(detections) - tp fn len(annotations) - tp sensitivity tp / len(annotations) if annotations else 0.0 ppv tp / len(detections) if detections else 0.0 return {TP: tp, FP: fp, FN: fn, Sensitivity: sensitivity, PPV: ppv}逻辑说明这里做了双重去重一个标注只能被一个检测匹配一个检测也只能匹配一个标注。最朴素的双重循环在数据量小时没有问题一条记录通常只有几百次心跳所以不需要优化。容差 150 ms 是一个在临床上被普遍接受的误差范围它既没有严格到把一点点相位偏移都算错也没有宽松到把误差很大的检测也放过去。参数说明tolerance_sec 我习惯取 0.15但如果你做的是心率变异性分析对 R 波位置的精确度要求更高建议收紧到 0.05 到 0.08 秒。另外如果你用的滤波器阶数较高或波形形态特殊检测位置可能系统性地偏移几十毫秒这种情况先检查滤波相位而不是放大容差掩盖问题。4.3 没有标注的数据怎么验证现实项目里大部分数据没有专家标注这时候可以换一种思路不验证每个 R 波位置而是验证 R-R 间期序列的合理性。正常心跳的 R-R 间期应该在 600 到 1500 ms 之间且相邻间期变化不会超过 30% 左右。如果算法测出来的 R-R 间期出现大量小于 300 ms 的值说明检测到了 T 波如果间期频繁大于 2000 ms说明漏检了心跳。def check_rr_intervals(r_peaks, fs, min_rr0.3, max_rr2.0): 用 R-R 间期做合理性检查返回可疑片段。 r_peaks: R 波位置索引列表 fs: 采样率 rr np.diff(r_peaks) / fs # 单位换算为秒 problems [] for i, interval in enumerate(rr): if interval min_rr or interval max_rr: problems.append((i, interval)) # 统计异常占比 anomaly_ratio len(problems) / len(rr) if len(rr) else 0.0 return problems, anomaly_ratio逻辑说明这段代码把相邻 R 波位置的差值从采样点换算成秒然后检查是否落在合理范围内。异常占比如果超过 5%通常说明检测参数不适合这条数据而不仅仅是个别噪声。这个方法不能帮你确认检测到的位置百分百正确但能快速把明显不合格的结果筛选掉是拿到新数据后第一个应该跑的脚本。5. 用 Python 实现心跳检测的 5 个避坑记录现象、原因、解决5.1 基线漂移导致阈值整体失效现象滤波后的信号整体仍然像一条波浪线R 波在波峰处被检测到在波谷处被漏掉检测结果忽多忽少。原因高通截止频率设太低假设设成 0.5 Hz其实呼吸引起的基线漂移仍在通带内滤波没有把它们抑制干净。解决把 lowcut 从 0.5 Hz 提到 5 Hz。基线漂移频率集中在 0.5 Hz 以下5 Hz 的高通虽然会损失一部分 ST 段信息但对 R 波检测毫无影响。这个坑最容易在采样率 500 Hz 以上的数据里出现因为低频成分在长记录中特别明显。5.2 T 波被当成 R 波心率直接翻倍现象检测到的心率大约是真实值的两倍画出的检测点落在波形的小鼓包上。原因T 波在滤波后仍然有一定幅度如果信号能量曲线中 T 波的峰值超过了阈值且不在不应期内就会被当作一次心跳。T 波离 R 波通常有 300 ms 左右的间隔恰好落在 200 ms 不应期之外。解决优先检查 highcut如果 highcut 是 20 HzT 波幅度明显、误检频繁就把 highcut 降到 12 Hz如果还有问题再把 threshold_ratio 从 0.4 提高到 0.5。T 波比 R 波变化缓慢斜率特征上天然有区分度这也是为什么先做差分而不是直接用原始波形来找峰。5.3 滤波后 R 波位置整体偏移精度怎么都提不上去现象评估结果中 TP 很多但位置平均偏差超过 30 ms拿去算心率没问题算心率变异性就超差。原因用了 lfilter 而不是 filtfilt。lfilter 是因果滤波输出必然滞后于输入滞后量跟滤波器阶数相关阶数越高延迟越大。这个偏差在单个波形上肉眼看不出来但和标注一对比立刻暴露。解决换成 filtfilt 做零相位滤波代码只有一处变化但精度提升是质的。如果你用的是实时流式数据filtfilt 因为需要未来数据不可用就只能接受延迟并在时间戳上做补偿。5.4 同一个心跳被报了两次双峰问题现象算法检测出的 R 波之间有大量 80 到 100 ms 的极短间隔画出来发现同一个 QRS 波上升沿和下降沿各触发了一次。原因差分平方之后QRS 波群在上升沿和下降沿都产生了明显的能量峰如果滑动窗口时间小于 QRS 本身宽度这两个峰就不会融合成一个。解决把 win_sec 从 0.10 提到 0.15 秒。如果数据采样率特别高比如 1000 Hz还要确认 win_len 至少覆盖 100 ms 以上。另外不应期的存在能挡住第二峰但不应期会把漏在两个峰之间的真实心跳也挡掉所以根本解法还是窗口宽度正确。5.5 采样率假设写死换数据就翻车现象同一套代码在一条数据上跑得很好换一条数据后异常率飙升检查发现采样率从 360 Hz 换成了 250 Hz而滤波窗口和滑动积分窗口全部用固定点数。原因代码里写死了窗口长度没有除以实际的 fs。200 ms 不应期在 360 Hz 下是 72 个点在 250 Hz 下是 50 个点差 22 个点足以导致完全不同的检测行为。解决所有时间相关参数全部用 fs 换算不要在函数里写 magic number。这是我见过发生频率最高的问题尤其是从网上复制代码时别人的数据采样率跟你不一定相同。6. 进阶验证技巧用合成信号自检再批量扫参数最有效的自检工具是先构造一段标准的仿真心电信号比如用几组高斯函数模拟 P 波、QRS 波和 T 波再叠加噪声和基线漂移然后用你的检测算法去测。合成信号的好处是每个 R 波位置是已知的算法表现可以被精确量化这比在真实数据上反复调整阈值效率高得多。要注意合成信号只是验证流程正确性不能替代真实数据评估因为合成波形往往比真实数据干净。批量调参是另一个值得投入的小工具。把 threshold_ratio 从 0.2 到 0.6 按步长 0.05 扫描把 win_sec 从 0.08 到 0.20 扫描每次跑完整条记录记录灵敏度和阳性预测值最后选两个指标同时最高的组合。这个扫描在单条记录上花费的时间只有几秒到几十秒换来的却是参数选择的确定性。我自己做每个新项目时都会把这条数据跑一遍记录参数和结果而不是靠感觉和运气去调因为真实数据噪声形态千差万别你没法凭空猜出最合适的阈值。最后再补一句任何检测算法跑新数据集都要先画一次波形亲眼确认检测点落在心拍上这个过程不能省希望帮到你。本文还有配套的精品资源点击获取
