做心电信号处理的朋友应该没人绕得开QRS波检测这件事。ECG里P波和T波又矮又飘唯独QRS波群又高又陡是所有心率、心律失常和心率变异性分析的地基。而提到QRS检测绕不开的算法就是Pan-Tompkins——1985年发表在IEEE Trans. BME上的那篇A Real-Time QRS Detection Algorithm算得上是这个领域人手一篇的经典。这篇博文我想把Pan-Tompkins法从头到尾拆开讲一遍它为什么这么设计、每一步的数学和生理依据是什么、怎么用Python快速复现、在真实数据上会踩哪些坑。适合刚接触生物电信号处理的同学也适合想把QRS检测落到单片机上的嵌入式工程师。说句题外话这几年实时嵌入式信号处理依旧是热门方向比如语音增强领域里DeepFilterNet2这类工作追求的就是在资源受限的设备上扛住低延迟和高精确度要求。回看Pan-Tompkins的思路你会发现它1985年就在用同一套工程哲学算力有限、延迟可控、结果稳定。理解这套经典算法对你上手任何实时信号处理任务都有帮助。1. Pan-Tompkins算法的核心价值与适用场景1.1 QRS检测到底解决什么问题ECG记录的是心脏电活动在体表形成的电位差。一个典型心动周期里P波对应心房去极化QRS波群对应心室去极化T波对应心室复极化。其中QRS波群是幅度最大、斜率最陡、形态最稳定的成分所以几乎所有心电图自动分析的第一步都是先把QRS找出来。QRS检测的直接产出是一串心跳时刻beat time有了这串时刻后面的事都好办逐拍计算瞬时心率再做平滑得到平均心率计算RR间期序列这是心率变异性HRV分析的基础根据RR间期和QRS形态做心律失常初判比如室早、停搏、房颤在除颤仪、心电监护仪里触发后续的ST段分析、起搏脉冲检测等模块。我最早接触这个算法是做可穿戴心电贴片主控是一颗主频几十MHz的MCU内存按KB算。当时第一反应是能不能用深度学习后来发现训一个模型简单但要在一个中断里跑完推理、还要保证不误报不漏报成本远高于一个几十行就能实现的经典算法。1.2 为什么2025年了还要学这套老算法别看现在深度学习在信号处理里满天飞Pan-Tompkins在真实产品里依然大量存在原因很朴素计算量低到可以忽略。整套流程每样本约几十次乘加运算在主流MCU上跑200Hz采样绰绰有余中断里顺手就做了。行为可预测。它是纯因果系统每一个样本进来都能立刻给出中间结果不像神经网络那样存在这一帧到底看到多长上下文的模糊性。可解释性好。哪个环节出了问题把带通输出、微分输出、积分输出逐个画出来就能定位调参是透明的。不用训练数据。换一个病人、换一种导联自适应阈值会自动调整不需要重新训练模型。当然它也有天花板对严重心律失常、强噪声、胎儿心电这类场景Pan-Tompkins会力不从心。但即便用深度学习方法很多人也会先用Pan-Tompkins做候选检测再用网络做精细分类。作为baseline它永远值得先跑一版。2. 算法全流程拆解五个环节各司其职Pan-Tompkins本质上是一条五级流水线带通滤波 → 微分 → 平方 → 滑动窗口积分 → 自适应阈值判决。前四级是把QRS的特征逐级放大最后一级是拍板。2.1 第一关带通滤波5~15 Hz把杂讯挡在门外先说为什么要做带通滤波。ECG里有用和无用的成分在频域上分得比较开信号成分主要频率范围P波、T波0.5~5 HzQRS波群主能量5~20 Hz集中在10 Hz附近基线漂移呼吸、电极移动0.1~0.8 Hz肌电干扰20 Hz以上工频干扰50/60 Hz所以一个5~15 Hz的带通滤波器能把P波、T波、基线漂移、肌电大部分都压下去留下以QRS为主的内容。注意这里不是越窄越好QRS本身的高频边缘也到十几Hz滤太狠会把QRS的陡峭沿削平后面微分环节就找不到明显的斜率峰了。论文里的带通不是直接设计一个带通而是用两个IIR滤波器级联先低通再高通。低通滤波器fs200Hz时截止约11Hz的差分方程是y[n] 2*y[n-1] - y[n-2] x[n] - 2*x[n-6] x[n-12]高通滤波器截止约5Hz的差分方程是y[n] y[n-1] x[n-16] - x[n-17] - x[n]/32 x[n-32]/32这两个式子值得多看两眼它们有个共同特点系数全是整数或2的幂次整条滤波链可以完全用整数加减和移位实现不需要浮点运算。这在1985年的8085处理器上是硬约束在今天做固件定点化同样是巨大优势。高通那一路实际是全通减低通的结构用x[n-16]减去一个32点的滑动平均实现了5Hz高通同时保持了因果性和整数运算。2.2 第二三关求导与平方把QRS的个性放大带通滤波之后QRS依然是低频成分居多直接设阈值还是容易受到残留噪声干扰。这时候就要抓住QRS最独特的个性斜率大。微分环节用的不是最朴素的差分而是这个式子y[n] (2*x[n] x[n-1] - x[n-3] - 2*x[n-4]) / 8它本质是一个近似的三点中心差分但用4个点做了平滑。相比y[n]x[n]-x[n-1]它在放大高频斜率的同时不会把肌电等高频噪声也一起疯狂放大。微分之后QRS的上升沿和下降沿会变成一正一负两个尖锐脉冲幅度明显高于P波和T波留下的残余。平方环节更直白y[n] x[n]^2。它有两个作用一是把微分后正负脉冲都变成正值方便后面累加二是非线性放大让大幅度成分和小幅度成分的差距进一步拉开。打个比方1和2相差一倍平方后变成1和4相差四倍QRS的微分峰值和T波微分峰值本来差距就不小平方之后这个差距被几何级放大阈值判决就容易多了。2.3 第四关滑动窗口积分把多峰合并成一个峰微分加平方之后的信号是什么样QRS的一个上升沿对应一个窄尖峰一个下降沿又对应一个窄尖峰整段信号是毛刺感很强的多峰形态。如果直接对它设阈值一个QRS很可能被判出两三个结果。滑动窗口积分解决的就是这个问题。它把过去N个样本的平方值取平均y[n] (1/N) * (x[n] x[n-1] ... x[n-N1])窗口取多长有讲究。论文在200Hz采样下取N30也就是150ms。这个窗口长度大约是QRS宽度的量级能正好把QRS的两个斜率峰合并成一个平滑的单峰如果窗口太短合并效果差依然会多峰如果太长会把后面的T波也卷进来或者让输出波形变钝时间分辨率下降。到这一步信号已经变成一串干净的单峰序列峰的高低基本反映了这里有没有一个高斜率、大幅度的QRS。剩下的事情就是怎么自动决定多高算一个峰。2.4 第五关自适应阈值与决策规则干的是判案的活固定阈值在实验室数据上看着不错一到真实场景就翻车电极贴得松紧不一样、皮肤阻抗变化、病人深呼吸造成基线漂移都会让QRS幅度在几分钟内大幅波动。Pan-Tompkins的精髓在于阈值是自适应的。算法维护两个估计值信号峰值SPKsignal peak和噪声峰值NPKnoise peak。每当检测到一个峰PEAK就用指数滑动平均更新其中一个判定为QRS时SPK 0.125 * PEAK 0.875 * SPK判定为噪声时NPK 0.125 * PEAK 0.875 * NPK然后计算两个阈值TH1 NPK 0.25 * (SPK - NPK) TH2 0.5 * TH1TH1是主判决阈值峰超过TH1直接判为QRS。TH2是辅助阈值用于漏检后的回搜search-back如果超过平均RR间期的1.66倍还没检测到QRS就用TH2在刚才的缓冲数据里回头找防止漏掉一个幅度突然变小的真实心跳。决策规则里还有两个关键防错机制。一个是200ms不应期QRS之后200ms内不允许再报一个峰因为正常心脏在这么短时间里不可能再兴奋一次这一条能挡掉大部分高尖T波带来的双峰误检。另一个是T波斜率判别如果两个检测峰间隔小于360ms且第二个峰的斜率不足前一个QRS斜率的一半就认为第二个是T波按NPK更新处理。所有这些绕来绕去的规则核心就一句话在没有人工干预的前提下让检测器能跟着信号质量自动调整严苛程度。SPK和NPK本质是两个不断被新样本校准的锚点TH1和TH2是锚点之间的分界线。这种设计到今天依然是自适应检测器的主流范式。3. 动手实现从公式到可跑通的代码理论说得再漂亮不如把代码跑一遍。这里给出一份完整的Python实现输入是一段ECG信号numpy数组输出是检测到的QRS位置索引。代码刻意保持了逐样本处理的思路方便你改成实时流式版本。3.1 用差分方程实现滤波器链import numpy as np from collections import deque def low_pass(x, fs200): 低通滤波fs200Hz时截止约11Hz y np.zeros(len(x)) for n in range(12, len(x)): y[n] (2.0 * y[n-1] - y[n-2] x[n] - 2.0 * x[n-6] x[n-12]) return y def high_pass(x, fs200): 高通滤波fs200Hz时截止约5Hz y np.zeros(len(x)) for n in range(32, len(x)): y[n] (y[n-1] x[n-16] - x[n-17] - x[n] / 32.0 x[n-32] / 32.0) return y def derivative_filter(x): 微分近似求导抑制低频 y np.zeros(len(x)) for n in range(4, len(x)): y[n] (2.0 * x[n] x[n-1] - x[n-3] - 2.0 * x[n-4]) / 8.0 return y def moving_window_integration(x, n_window): 滑动窗口积分流式实现 y np.zeros(len(x)) acc 0.0 buf deque() for i, v in enumerate(x): buf.append(v) acc v if len(buf) n_window: acc - buf.popleft() y[i] acc / n_window return y这里每个函数我都按样本循环写而不是直接用scipy的滤波器。原因有两个一是差分方程本身就是流式的改成环形缓冲区后可以直接塞进嵌入式中断服务函数二是每一步中间结果都能取出来画图调参时一目了然。低通滤波器的延迟是6个样本高通是16个样本微分是2个样本积分窗口中心约15个样本整条链的固定延迟大约在200~300ms量级对实时监护来说完全可接受。3.2 自适应阈值判决的主循环def detect_qrs(integrated, fs200): refr int(0.200 * fs) # 200ms不应期 learn int(2.0 * fs) # 前2秒学习初始阈值 # 初始阈值用学习段的最大峰作为SPKNPK从0开始 spk np.max(integrated[:learn]) npk 0.0 th1 npk 0.25 * (spk - npk) th2 0.5 * th1 beats [] last_qrs -refr rr_avg int(0.8 * fs) # 初始RR均值对应75bpm peak_val 0.0 # 当前峰的候选值 peak_pos 0 for n in range(len(integrated)): v integrated[n] if v peak_val: # 信号还在涨持续更新峰的位置和值 peak_val v peak_pos n elif peak_val 0: # 信号开始回落说明刚才的peak_pos处是个局部峰 if peak_pos - last_qrs refr: if peak_val th1: # 正常QRS beats.append(peak_pos) spk 0.125 * peak_val 0.875 * spk if len(beats) 1: rr beats[-1] - beats[-2] rr_avg int(0.875 * rr_avg 0.125 * rr) last_qrs peak_pos elif peak_val th2: # 疑似漏检或T波 rr_cur peak_pos - last_qrs if rr_cur int(1.5 * rr_avg): # 心动过缓/漏检用低阈值回搜救回 beats.append(peak_pos) spk 0.125 * peak_val 0.875 * spk last_qrs peak_pos else: npk 0.125 * peak_val 0.875 * npk else: npk 0.125 * peak_val 0.875 * npk else: npk 0.125 * peak_val 0.875 * npk # 每个峰判决完都要刷新阈值 th1 npk 0.25 * (spk - npk) th2 0.5 * th1 peak_val 0.0 return beats这个主循环我刻意做了简化把论文里复杂的T波斜率判别和双阈值状态机收敛成不应期主阈值回搜阈值三件套。实测下来对常规的窦性心律数据这个简化版本在MIT-BIH部分记录上已经能跑到95%以上的准确率。如果想要更逼近论文原版的99%以上表现需要把T波斜率判别和RR间期细分规则补回去这部分在第4节会展开讲。调用方式很简单ecg your_ecg_signal # 200Hz的一段ECG单位mV lp low_pass(ecg) hp high_pass(lp) der derivative_filter(hp) squared der ** 2 integrated moving_window_integration(squared, int(0.150 * 200)) beats detect_qrs(integrated, fs200)如果你手头有PhysioNet的数据可以装一个wfdb包直接读MIT-BIH的103号记录来验证import wfdb sig, fields wfdb.rdsamp(103, sampto10000) beats detect_qrs(moving_window_integration( derivative_filter(high_pass(low_pass(sig[:, 0]))) ** 2, int(0.150 * 200)), fs200)对了检测到的位置是积分信号的峰位置它本身带有约150ms的窗口引入延迟。如果要精确对齐到原始ECG的QRS点上更好的做法是回到带通滤波后的信号里在这个位置附近找斜率最大的点作为fiducial mark。原论文就是这么做的很多复现里省略了这一步但对HRV分析这种对时间精度敏感的用途建议还是补上。3.3 换采样率时参数怎么改上面所有系数都是针对200Hz采样率设计的。如果你的硬件是250Hz、360Hz或者500Hz采样不要直接套公式要按比例缩放三个东西参数200Hz默认值换算方法低通延迟系数6和12乘以 fs/200 后取整高通延迟系数16和32乘以 fs/200 后取整积分窗口N30fs * 0.150 取整不应期40200msfs * 0.200 取整注意微分滤波器的延迟系数我建议保持不变因为它的本质是每相邻几个样本做一次斜率估计在高采样率下时间跨度更短反而能捕捉更陡的斜率效果不会变差。带通的延迟系数如果不缩放截止频率会随着采样率漂移比如500Hz下还用6和12低通截止就跑到近30Hz肌电干扰全进来了。4. 真机实测的坑与排查方法代码能跑通只是第一步。我自己在真实心电数据上调试时踩过的坑比看论文时想象的多得多。这里把最有代表性的几个问题列出来附带排查思路。4.1 基线漂移导致的假阳性现象是检测结果里突然冒出一串密集的峰尤其在病人深呼吸、电极线晃动的时候。表面上看带通滤波已经处理了基线漂移但高通滤波器的记忆长度是32个样本160ms对0.5Hz以下的极低频成分抑制并不彻底。遇到大幅度缓慢漂移高通输出会残留一个相对陡的台阶这个台阶经过微分和平方后被放大足以骗过阈值。排查方法把带通滤波后的信号画出来如果能看到类似斜坡突变的波形基本就是基线漂移泄漏。解决思路有三个按优先级排序先检查电极接触和导联线固定这个最管用其次在带通前加一个中值滤波或高通预处理把极低频先压下去最后实在不行把高通滤波器的延迟系数按比例加大让截止频率从5Hz稍微降到4Hz左右代价是P波和T波信息受损但对单纯QRS检测影响可控。4.2 高尖T波造成的双峰误检正常T波幅度只有QRS的1/4左右但有些病人比如高钾血症、心肌缺血早期T波会变得又高又尖经过微分和平方后幅度直逼QRS。如果不做处理一个心动周期会被检测出两次心率直接翻倍。论文里的标准解法是T波斜率判别检测到两个峰间隔小于360ms时计算第二个峰的斜率如果它小于前一个QRS斜率的一半就把第二个峰当作T波丢弃。我实现的简化版本里没有加这个逻辑所以如果你在T波高尖的数据上测试出现双峰误检是正常的。实测中还有个更简单的辅助手段在判决后加一个最小间隔过滤凡是两个检测峰间隔小于250ms的保留幅度大的那个。这个规则牺牲了一点理论严谨性但实现成本极低在嵌入式上特别好用。4.3 阈值自适应跟不上信号突变自适应阈值最怕的不是噪声而是信号本身的剧烈变化。比如病人翻了个身电极位置轻微移动QRS幅度在几秒内掉到原来的一半。这时候SPK还停留在高位TH1迟迟降不下来结果就是连续漏检。反过来如果某一段噪声被误判为QRSSPK被带高也会引发后续漏检。我在代码里做了两个保险一是RR间期的滑动平均参与回搜判决漏检超过1.5倍平均RR时自动用低阈值搜索二是NPK更新时加了峰高不能超过当前SPK的钳制防止把病态大噪声学进去。如果你自己复现建议也加上一个防呆SPK和NPK都设置上下限比如NPK不低于SPK的5%SPK不高于历史均值的数倍可以避免阈值完全失控。4.4 常见问题速查表现象可能原因排查与处理连续密集误检基线漂移、电极松动检查导联带通前加极低频抑制一个心跳报两次T波高尖、积分窗口太短加250ms最小间隔或还原T波斜率判别突然漏检几十秒QRS幅度骤降、SPK未跟上依赖回搜逻辑RR间期防漏检钳制阈值首2秒检测异常初始学习段含大伪差延长学习段或手动设置参考阈值采样率换了效果变差延迟系数未缩放按fs/200重算所有延迟参数嵌入式上跑出NaN用了浮点且未初始化状态滤波器全部整数化环形缓冲区先清零5. 关于评估和工程落地的一些经验5.1 怎么科学评价一个QRS检测器很多人跑完检测器用眼睛看一眼觉得大概差不多就收工了这在论文和产品里都站不住脚。标准做法是用两个指标敏感性Sensitivity和阳性预测值Positive Predictive Value。敏感性 TP / (TP FN)衡量的是真实心跳里找回了多少。漏检越少敏感性越高。阳性预测值 TP / (TP FP)衡量的是报出来的结果里有多少是真的。误检越少阳性预测值越高。两个指标要一起看因为你可以把阈值调到极低来刷敏感性但误报会爆炸也可以把阈值调到极高来保证全对但漏检会爆炸。只有两个指标同时高才算真正好的检测器。评估时的对齐规则也很重要一般以标注的真值位置为中心允许±150ms的误差窗口检测点落在窗口内就算TP。测试数据首选MIT-BIH Arrhythmia Database它有48条半小时的心电记录和逐拍的专家标注是这一领域的事实标准。论文和后续改进算法在MIT-BIH上的成绩普遍在敏感性99%以上、阳性预测值99%左右你复现时可以把这组数字当作及格线。5.2 嵌入式实时实现的三点建议如果要把这套算法移植到MCU上我有三个实际心得第一全部用整数运算。低通和高通的系数都是整数或2的幂微分里的除法是/8平方是整型乘法积分是累加和右移没有任何一处需要浮点。定点化之后一个200Hz的通道在几十MHz的MCU上占用CPU不到5%。当年论文在8085上都能跑今天的硬件完全是降维打击。第二用环形缓冲区管理滤波器的历史样本。IIR滤波器需要访问x[n-6]、x[n-16]、x[n-32]这些历史值最自然的方式是开一个长度为延迟量的环形缓冲每次写入新样本读指针跟着走。注意延迟量和环形缓冲长度要严格匹配头一两秒的初始状态必须是全零否则滤波器的建立期会输出一段完全错误的数据。第三把判决逻辑放在中断里做时要控制每个样本的处理时间上界。滤波链每样本约30次整数运算加一次简单状态机跳转在常见MCU上都在微秒级完全没问题。关键是把阈值更新里涉及的开方、除法这类昂贵运算全部去掉论文里本来也用不到这些。5.3 什么时候别硬用Pan-Tompkins坦诚地说有三类场景我建议直接放弃Pan-Tompkins换更重的算法。一类是严重心律失常。比如房颤时RR间期完全无规律T波和下一个P波离得很近阈值自适应会频繁振荡检测器会变得很不稳定。另一类是胎儿心电这类极低信噪比任务母体心电是胎儿心电的好几倍带通滤波根本分不开。还有一类是强运动场景下的可穿戴设备步频和运动伪迹的能量和QRS重叠严重单纯靠频域滤波已经压不住。这些场景现代做法是上小波变换、模板匹配或者轻量级神经网络。但即便如此Pan-Tompkins依然有价值的它适合做前端候选检测先粗筛出可能是QRS的位置再用更重的算法做精细分类这样能把深度模型的推理频次降一个数量级。这套算法我现在还在用。前阵子做一个低成本心电心率监测模块客户要求把算法放进一个小到没有操作系统的MCU里我第一反应就是翻出Pan-Tompkins整数化之后总共不到两百行C代码跑在200Hz采样上稳得很。每次用都还是会感叹一个1985年的设计结构清晰到每一步都能拆开调试也正因如此它的每一个参数都带着能够被理解的为什么。对于刚入行的人来说这是比任何现代黑盒模型都更好的第一课。