做心电ECG信号处理的人几乎都绕不开QRS波峰检测这一步。心率计算、心率变异性分析、早搏识别、房颤筛查任何一种下游分析都要先知道每一次心跳发生的确切时刻。很多人一开始会觉得这事不简单吗在信号里找一个最高峰就行但真实心电一进来基线漂移、肌电噪声、形态各异的P波和T波会立刻把你这个想法打碎。Pan-Tompkins算法算得上是这个问题最经典也最实用的答案1985年发表在IEEE Transactions on Biomedical Engineering上用一串简单的级联滤波和一个自适应阈值决策在MIT-BIH数据库上把QRS检测做到了99%以上的敏感度。这篇文章我把它的每个模块拆开讲清楚给出200Hz采样率下的Python实现再聊几个只有真正动手跑过才会遇到的坑。适合刚接触心电信号处理的同学也适合想把现有检测方案替换成这个经典baseline的工程师。1. 为什么QRS检测是个难题噪声和波形都在和你作对1.1 心电波形里有哪些成分常规心电图里一次心跳主要由P波、QRS波群、T波三段组成。P波是心房除极产生的幅度很低通常在0.1到0.25mV之间上升沿平缓频谱能量主要分布在低频段。QRS波群是心室除极产生的幅度高一般能达到0.5到1.5mV时间宽度约60到110ms能量集中区大致在5到15Hz这一段是心电信号里最“陡”的部分。T波是心室复极产生的幅度中等频率比QRS低得多形态宽而缓。听起来很清楚但实际数据远比教科书复杂。不同导联、不同病人的QRS形态差异很大有的R波特别高有的R波很小S波很深还有双相、切迹、束支阻滞造成的宽大畸形QRS。如果只靠“找最高点”或者“找斜率最大的点”任何单一规则都会在形态变化面前失效。这也是为什么QRS检测不能简单套一个固定模板而是要先做信号处理把“QRS能量”和“其他东西”在特征空间里拉开距离。1.2 三类噪声怎么影响检测除了形态变化真实心电信号还叠加了大量干扰。最常见的三类是基线漂移、肌电干扰和工频干扰。干扰类型主要来源频谱范围对检测的影响基线漂移呼吸、电极松动1Hz幅度可达数mV直接把整个波形抬高或压低肌电干扰肌肉收缩、身体活动20Hz以上高频毛刺容易被当成假的波峰工频干扰电源环境50Hz/60Hz及谐波波形上叠一层细密锯齿基线漂移最隐蔽因为它的幅度可能比QRS还大。如果不去掉它阈值检测会随着基线忽高忽低而失灵。肌电干扰和运动伪差在可穿戴设备里尤其常见表现为高频的尖峰脉冲经过滤波后仍然可能形成杂散峰。工频干扰相对好处理Pan-Tompkins原算法干脆把低通截止频率设在11Hz附近直接把它压掉根本不需要额外陷波器。1.3 真正的敌人是P波和T波三噪声再麻烦本质上是“外来干扰”可以通过滤波器压制。真正麻烦的是心电信号里那些和QRS长得像的“内部成员”——尤其是T波。T波的幅度在正常情况下比R波低但在某些病理或生理状态下T波可以异常高尖甚至接近或超过R波比如高钾血症、早复极综合征等情况。P波虽然幅度低但如果滤波参数没调好或者信号质量差它也可能在积分波形上形成一个足以触发阈值的峰。更麻烦的是P波和T波在形态上有一定“欺骗性”单靠幅度、斜率都很难和低幅QRS区分。所以QRS检测从本质上说是在“低信噪比形态多变生理干扰”的三重压力下做事件检测。Pan-Tompkins的解决思路不是用一个更聪明的判断器而是先用级联滤波把QRS的能量集中成一个干净的单峰再用带生理约束的自适应阈值做决策。这个思路到今天看依然非常扎实。2. 五级处理链拆解每一级在压制什么、放大什么2.1 带通滤波用两个简单递归滤掉基线和肌电Pan-Tompkins算法的第一步是构造一个约5到15Hz的带通滤波器。这个频段是QRS能量最集中的区域同时又能抑制低频基线漂移、高频肌电和工频干扰。实现上作者没有用通用的FIR或IIR滤波器工具而是设计了一对极其简单的递归差分方程。低通部分的差分方程是y[n] 2*y[n-1] - y[n-2] x[n] - 2*x[n-6] x[n-12]这个滤波器在200Hz采样率下的截止频率约11Hz直流增益36。它的z域传递函数本质上是H(z) ((1 - z^-6) / (1 - z^-1))^2也就是一个6点滑动平均连续做了两次。滑动平均本身就是最简单的低通两次级联后带外衰减更好。最关键的是这个滤波器实现起来只需要加法和移位连乘法都不需要非常适合当时的实时硬件今天的MCU同样喜欢这种结构。高通部分的差分方程是y[n] y[n-1] - x[n]/32 x[n-16] - x[n-17] x[n-32]/32截止频率约5Hz作用是彻底干掉直流分量和0.5Hz附近的呼吸基线漂移。它的实现思路也很有意思构造了一个截止频率约5Hz的低通再用原始信号减去这个低通得到高通输出。1/32这个系数同样可以用移位完成。两个滤波器级联后输出信号里剩下的是以5到15Hz为主的QRS成分P波和T波的大部分能量在这个阶段已经被明显压制。注意这些系数都是针对200Hz采样率设计的如果采样率变了系数里的延迟项对应的时间就变了必须重新整定不能直接套。2.2 微分把斜率变成音量滤波之后信号平滑了但QRS和P/T波的区分度还不够。下一步是微分目的是提取波形的斜率信息。QRS的上升沿和下降沿非常陡微分后会形成一个大幅度的正脉冲和负脉冲P波和T波宽而缓微分后幅度明显小。论文里用的五点差分公式是y[n] (1/8) * (2*x[n] x[n-1] - x[n-3] - 2*x[n-4])这个式子相当于在4个采样间隔上做了近似中心差分在200Hz采样率下能较好地覆盖20到30Hz范围的斜率成分而且系数都是2的幂除以8用移位就能完成。很多开源实现里直接用最简单的中心差分(x[n1] - x[n-1]) / 2替代实际差别不大。我自己对比过两种方式R峰定位结果差异在毫秒级对于绝大多数应用可以忽略但要注意微分环节本身会引入几个采样点的相位延迟。2.3 平方与积分让每个QRS在波形上“变成一个峰”微分后的信号有正有负而且一个QRS会产生两个方向的尖峰。接下来做逐点平方这一步有两个作用一是把所有值变成非负方便后续的峰检测二是非线性放大差异。如果QRS微分后的幅度是P波微分幅度的3到4倍平方之后能量差距就变成9到16倍QRS的优势被显著拉开。平方后的信号仍然很碎一个QRS可能产生多个相邻的尖峰直接检测容易出现重复计数。所以最后一步是移动窗口积分y[n] (1/N) * sum(x[n-i] for i in range(N))N在200Hz采样率下取30对应150ms窗口。这个窗口的物理含义很关键它要把一个QRS对应的多个子峰合并成一个宽的平滑峰同时尽量不把T波包进来。窗口太小比如20ms合并效果差还是会出现多个峰窗口太大比如250ms会把T波甚至下一个P波卷进同一个峰里导致两个心拍被合并成一次。经过这五级处理后原始ECG被转换成一个近似单峰的“事件指示信号”后续的检测就完全在这条平滑波形上进行了。这个逐步抽象的思路值得记住不要在原始信号上直接设阈值先用信号处理把问题简化。3. 阈值与决策为什么偏偏是两套阈值加一个不应期3.1 单阈值在真实信号上必然失效滤波链只能让QRS能量更集中不能保证每个QRS幅度都一样。不同病人、不同导联、不同电极接触状态QRS幅度可能相差几倍。同一个病人也会因为呼吸、体位变动QRS幅度发生周期性变化。如果用一个固定阈值阈值设高了低幅QRS漏检设低了噪声峰和T波峰全部触发。所以Pan-Tompkins的阈值必须是自适应的而且要区分“信号峰”和“噪声峰”两套统计量分别跟踪、分别更新。这是整个算法里工程味道最浓的部分也最容易被初学者忽略。3.2 SPK/NPK动态更新规则算法维护两个关键状态量SPKsignal peak信号峰估计和NPKnoise peak噪声峰估计。它们的更新公式是SPK 0.125 * PEAK 0.875 * SPK # 当前峰判定为QRS时 NPK 0.125 * PEAK 0.875 * NPK # 当前峰判定为噪声时PEAK是积分信号上检测到的局部峰幅度。0.125和0.875这一组系数本质上是低通滤波等价于只给当前观测值12.5%的权重历史估计值87.5%的权重。这样设计的好处是阈值的更新速度跟得上QRS幅度的渐变但又不会被某一次异常的大噪声峰瞬间带偏。每一轮峰判定结束后用新的SPK和NPK计算两个阈值TH1 NPK 0.25 * (SPK - NPK) TH2 0.5 * TH1TH1是主检测阈值插在噪声水平和信号水平之间25%的位置偏向“宁可多检也不漏检”。TH2是TH1的一半作为低灵敏度档位的备用阈值后面会用到。这个双阈值结构是整个决策逻辑的基础。3.3 不应期和360ms低阈值回看第一个关键生理约束是不应期。一旦检测到一个QRS接下来200ms内不再接受任何新的峰。原因是人心脏的生理上限约300bpm对应RR间期200ms如果两个峰相隔不到200ms物理上不可能是两次独立心跳只能是同一个心拍留下的过冲或T波残留。这个约束直接消灭了滤波器无法完全消除的T波误检和重复计数。第二个关键逻辑是低阈值回看。检测过程中有些峰落在TH2和TH1之间单看不满足主阈值条件但可能是真实的QRS。怎么区分看时间。如果这个峰距离上一个QRS已经超过360ms那它大概率的的确确是一个心跳因为正常的T波不会离上一个QRS这么远。此时用TH2这个低阈值再判一次接受这个峰为QRS。这个机制解决的是高T波场景下的漏检问题。T波如果很高它在积分波形上也是一个明显的峰但因为它紧跟QRS会被不应期挡住。如果不应期结束后没有马上检测到下一个QRS而T波峰又刚好落在TH2~TH1之间那么360ms回看就能把可能是真实心跳的事件捞回来。决策层这几条规则必须配合使用缺一个都可能在某些场景下出问题。3.4 RR间隔异常时的search-back还有一类更隐蔽的漏检某段信号噪声很大或者幅度突变把SPK抬得过高导致后续真正QRS幅度低于TH1连续好几个心拍都没检测到。这时候需要第三层保护——RR间隔回溯。算法维护一个平均RR间期估计更新方式和SPK一样采用一阶低通RR_avg 0.125 * (RR_new - RR_avg) RR_avg RR_low 0.92 * RR_avg RR_high 1.66 * RR_avg如果距离上一次QRS的时间已经超过RR_high说明大概率漏检了。此时回到最近的信号缓冲区在里面找积分波形的最大峰。如果这个峰超过TH2就把它认定为漏掉的QRS并据此更新阈值。这个机制能在心率突然变化或者信号暂时劣化时自动恢复避免一整段信号全部丢失。这里有个容易误解的地方RR_high的1.66倍和前面说的360ms回看不是一回事。360ms回看是处理“单峰落在低阈值区间”的日常情况search-back是处理“长时间没有新峰”的系统性漏检。前者是细调后者是救场两者相互补充。4. 200Hz采样率下的Python实现从滤波链到R峰定位4.1 滤波器系数与初始化下面这套代码可以直接作为基础版本使用。用scipy.signal.lfilter实现差分方程逻辑清晰也方便后续改成流式处理。import numpy as np from scipy.signal import lfilter fs 200 # 低通截止约11Hzfs200 b_lp np.zeros(13) b_lp[[0, 6, 12]] [1, -2, 1] a_lp [1, -2, 1] # 高通截止约5Hzfs200 b_hp np.zeros(33) b_hp[[0, 16, 17, 32]] [-1/32, 1, -1, 1/32] a_hp [1, -1] def bandpass(x): y lfilter(b_lp, a_lp, x) y lfilter(b_hp, a_hp, y) return y # 中心差分 def derivative(x): y np.zeros_like(x) y[1:-1] (x[2:] - x[:-2]) / 2 return y def classify_signal(x): x_lp bandpass(x) x_diff derivative(x_lp) x_sq x_diff ** 2 window np.ones(30) / 30 mwi np.convolve(x_sq, window, modesame) return mwi初始化阶段我习惯用处理后信号前2秒的最大值作为SPK初始值NPK初始为0然后让算法自己快速收敛。这跟论文原文的初始值取法略有差异不同复现版本也不完全一致但实测差别不大前提是前2秒信号质量要可靠。4.2 峰检测主循环怎么写主循环是决策层的核心。下面是一个便于阅读的示意版本省略了search-back的完整实现保留了不应期、双阈值判断和阈值更新def detect_qrs(mwi, x_raw, fs200): blank int(0.200 * fs) # 不应期 rr_back int(0.360 * fs) # 低阈值回看时间 spk np.max(mwi[:fs * 2]) npk 0.0 th1 npk 0.25 * (spk - npk) th2 0.5 * th1 last_qrs -blank peaks [] peak_val 0.0 peak_idx -1 for n in range(len(mwi)): if mwi[n] peak_val: peak_val mwi[n] peak_idx n elif mwi[n] peak_val and n - peak_idx 2: # 已经跨过局部峰开始判定 if peak_idx - last_qrs blank: if peak_val th1: accept True elif peak_val th2 and peak_idx - last_qrs rr_back: accept True else: accept False if accept: peaks.append(peak_idx) spk 0.125 * peak_val 0.875 * spk last_qrs peak_idx else: npk 0.125 * peak_val 0.875 * npk th1 npk 0.25 * (spk - npk) th2 0.5 * th1 peak_val 0.0 peak_idx n return np.array(peaks)代码里有一个细节要注意局部峰判定必须等到波形开始下降才确认不能一出现新高就清零状态否则会漏掉大峰前面的小峰。这个“滞后确认”的方法虽然简单但在流式处理里很关键。4.3 把积分峰位置换算成真实R峰千万不要把积分波形的峰位置直接当作R峰时间。移动窗口积分本身会带来约150ms/2的滞后再加上带通滤波器的群延迟直接用的结果比真实R峰位置偏后几十毫秒对心率计算影响不大但对HRV分析、QT间期测量是致命的。工程上最常用的做法是检测到积分峰后回到原始信号或带通滤波后的信号在积分峰位置附近一个短窗口内找幅值最大的极值点作为真正的R峰位置。窗口长度一般取60到100ms。论文里的做法更精细一些会先找斜率最大点再回溯到尖峰但对多数场景最大幅值搜索已经足够。def refine_r_position(mwi_peak, x_raw, fs): half_win int(0.05 * fs) lo max(0, mwi_peak - half_win) hi min(len(x_raw), mwi_peak half_win) return lo int(np.argmax(np.abs(x_raw[lo:hi])))注意滤波后信号和原始信号在R峰位置附近可能极性不同有的导联R波是负向的所以要取绝对值最大。这也是“找极大值”时要格外小心的地方。4.4 流式处理与离线批处理的区别Pan-Tompkins算法从设计之初就是流式处理的样本逐个进入滤波链不断更新状态变量阈值和不应期始终在线计算。真正的嵌入式实现里低通需要保存最近12个输入输出样本高通需要保存最近32个样本积分窗需要保存最近30个平方值。这些状态量加起来不到一百个变量内存开销极小。很多Python版本的实现是批处理模式也就是先把整段信号滤波完再做峰检测好处是调参方便坏处是会让初学者忽略状态保存和因果性要求。如果你做离线分析用lfilter和convolve完全没问题但如果想验证算法在实时场景下的行为建议把代码改成逐样本处理保留滤波状态别用filtfilt这种零相位滤波——它引入了未来信息不是实时的语义确定的R峰时间会和在线系统不一致。5. 在MIT-BIH数据上的实测记录哪些参数最容易让人翻车5.1 评估指标怎么算Se与P验证Pan-Tompkins实现最常用的公开数据是MIT-BIH Arrhythmia DatabasePhysioNet上可以获取包含48条约30分钟的双通道心电图记录有两位专家独立标注的R峰位置。评估时把检测到的R峰和标注峰做匹配落在容忍窗口内的才算命中窗口常用±100ms或±150ms。两个核心指标Se TP / (TP FN) # 敏感度漏检越少越高 P TP / (TP FP) # 阳性预测值误检越少越高原作者论文在MIT-BIH上报告的结果在99.3%敏感度、99.7%左右阳性预测值。如果你的复现跑下来正常质量记录也能在99%附近说明实现基本正确如果明显偏低问题大概率出在阈值更新或不应期逻辑而不是滤波链。5.2 采样率一变整个滤波器系数都要重新设计这是我见过最多人踩的坑。有人拿到360Hz采样率的数据直接把200Hz的滤波器系数带进去跑结果一团糟。原因很简单低通差分方程里的6和12是采样点数量不是物理时间。200Hz下它们对应30ms和60ms到了360Hz下就变成16.7ms和33.3ms滤波器实际截止频率完全变了QRS的主要能量可能直接被滤掉。正确做法有两个一是用scipy.signal.resample或resample_poly把数据重采样到200Hz再跑原版系数二是用scipy.signal.butter重新设计5到15Hz的带通滤波器替代原来的整数系数滤波链。前者保住了Pan-Tompkins的参数体系后者更灵活但失去了整数运算的便利。在工程上我一般优先重采样除非存储和算力实在不允许。5.3 初始化前2秒的质量决定了整段检测的上限阈值初始化依赖前2秒的处理信号。如果这2秒里包含大量噪声、导联脱落段、或者T波异常高初始SPK和NPK就会偏离正常范围后面很长一段时间阈值都难以回到正轨。我排查问题时发现很多“算法偶尔整段失效”的场景根因都是初始化段数据质量不好。一个实用的检查方法把积分信号画出来确认前2秒至少有一个清晰可见的QRS峰。如果初始化段本身就有问题直接把这2秒丢弃用后面信号重新初始化效果会好很多。产品化时更稳妥的做法是配合信号质量指数SQI判断初始化期间信号质量太差就一直等待不启动检测。5.4 高T波和低幅QRS会把阈值带偏我在一条记录里遇到过T波和R波几乎一样高的病例。简化版算法把T波连续误判成QRS结果SPK被T波峰撑高真正的QRS反而低于TH1又造成整段漏检。这种正反馈式的错误非常危险一旦SPK被污染QRS就会被持续漏掉。解决办法就是严格执行三层决策不应期挡住紧跟QRS的T波TH2回看把可能被遗漏的QRS捞回来search-back处理整段漏检。另外PEAK的取值源也值得注意。用积分波形的峰值比用原始信号峰值稳定得多因为积分已经平滑掉了高频毛刺。我第一个版本用原始信号峰值更新阈值遇到一段噪声稍大的数据就震荡换回积分波形峰值后稳定很多。5.5 R峰定位的时间偏移被很多人忽略了滤波链里每个环节都有延迟低通约6个采样点、高通约16个、微分数点、积分窗15点。累计下来积分波形峰相对于真实R峰可能偏后几十毫秒。如果直接拿积分峰当R峰心率计算不受影响但HRV分析里的RR间期序列会出现系统性偏置QT间期估计就更不可靠。我自己的习惯是每次检测到积分峰后必须回到原始信号做局部精细搜索把修正后的R峰位置存下来然后用MIT-BIH标注抽查几条记录确认位置误差在±10ms以内。这一步虽然简单却是从“能检测”到“能被专业分析使用”的分水岭。6. 落到硬件和真实产品前还要想清楚的几件事6.1 采样率适配重采样到200Hz还是重新设计滤波器如果设备ADC采样率不是200Hz而且无法修改面临一个选择重采样还是重设计滤波链。我的建议是优先重采样到200Hz。原因很简单Pan-Tompkins算法里除滤波系数外积分窗口、不应期、360ms回看时间、RR阈值更新等一系列逻辑都以200Hz为隐含前提。全部统一到200Hz之后调参只需考虑物理时间不容易出错。只有两种情况我会选择重新设计滤波器一是目标平台没有重采样条件二是在MCU上希望减少一次重采样带来的额外算力和内存。这时候可以用scipy.signal.butter(2, [5, 15], btypebandpass, fsfs)设计新的带通但要把系数转成适合定点的格式再进硬件。6.2 嵌入式定点化这个算法为什么适合裸机Pan-Tompkins能在四十多年后仍然是嵌入式心电检测的常青树很大程度上因为它的计算复杂度极低。整套滤波链的系数不是整数就是2的负幂次微分除8积分除30可以近似成除32再用小修正全部可以用整数移位实现不需要浮点单元。在Cortex-M3这个级别的MCU上逐样本处理200Hz数据的开销几乎可以忽略。只需要为低通保存12个历史样本、高通保存32个历史样本、积分窗保存30个历史平方值再加上SPK、NPK、RR相关状态内存占用不超过几十字节。相比跑一个神经网络模型这个资源开销完全是另一个量级。6.3 什么情况下该换深度学习方案Pan-Tompkins再经典也有边界。当信号来源特别复杂比如强运动场景下的可穿戴心电、胎儿心电、多病种畸形QRS或者需要同步精确定位P波起点、QRS起点、T波终点这些精细边界时经典算法的框架很难满足。工程上比较务实的做法是组合使用用Pan-Tompkins做第一道候选生成保证延迟低、覆盖率高再用一个小型神经网络或模板匹配做第二道分类过滤掉T波误检和噪声峰。这样既保留了经典算法的实时性和可解释性又能借用深度模型在复杂形态上的泛化能力。但无论怎样先用经典算法拿到baseline永远是最值的投入它能让你清楚知道问题到底出在检测层还是出在信号质量层这是直接上深度学习时很难获得的信息。整套算法跑通之后我最深的感受是Pan-Tompkins看起来只是几条差分方程和几个阈值公式但它背后把“信号特征生理约束工程可实现性”三个层面都考虑全了。1985年没有GPU没有Python生态作者靠的是对信号本质的理解。后来我自己做心电监护模块第一版检测逻辑都是先复现这个算法再在它的框架上加SQI判断和专科特征。如果你也是第一次接触QRS波峰检测我建议先别急着换方向把这个算法在标准数据库上完整跑一遍把阈值更新的过程画出来看你会从中得到比跑通一个深度学习模型更多的信号处理直觉。