做心電信號處理的朋友應該沒人繞得開QRS波檢測這件事。ECG里P波和T波又矮又飄唯獨QRS波群又高又陡是所有心率、心律失常和心率變異性分析的地基。而提到QRS檢測繞不開的算法就是Pan-Tompkins——1985年發(fā)表在IEEE Trans. BME上的那篇A Real-Time QRS Detection Algorithm算得上是這個領域人手一篇的經典。這篇博文我想把Pan-Tompkins法從頭到尾拆開講一遍它為什么這么設計、每一步的數(shù)學和生理依據(jù)是什么、怎么用Python快速復現(xiàn)、在真實數(shù)據(jù)上會踩哪些坑。適合剛接觸生物電信號處理的同學也適合想把QRS檢測落到單片機上的嵌入式工程師。說句題外話這幾年實時嵌入式信號處理依舊是熱門方向比如語音增強領域里DeepFilterNet2這類工作追求的就是在資源受限的設備上扛住低延遲和高精確度要求?；乜碢an-Tompkins的思路你會發(fā)現(xiàn)它1985年就在用同一套工程哲學算力有限、延遲可控、結果穩(wěn)定。理解這套經典算法對你上手任何實時信號處理任務都有幫助。1. Pan-Tompkins算法的核心價值與適用場景1.1 QRS檢測到底解決什么問題ECG記錄的是心臟電活動在體表形成的電位差。一個典型心動周期里P波對應心房去極化QRS波群對應心室去極化T波對應心室復極化。其中QRS波群是幅度最大、斜率最陡、形態(tài)最穩(wěn)定的成分所以幾乎所有心電圖自動分析的第一步都是先把QRS找出來。QRS檢測的直接產出是一串心跳時刻beat time有了這串時刻后面的事都好辦逐拍計算瞬時心率再做平滑得到平均心率計算RR間期序列這是心率變異性HRV分析的基礎根據(jù)RR間期和QRS形態(tài)做心律失常初判比如室早、停搏、房顫在除顫儀、心電監(jiān)護儀里觸發(fā)后續(xù)的ST段分析、起搏脈沖檢測等模塊。我最早接觸這個算法是做可穿戴心電貼片主控是一顆主頻幾十MHz的MCU內存按KB算。當時第一反應是能不能用深度學習后來發(fā)現(xiàn)訓一個模型簡單但要在一個中斷里跑完推理、還要保證不誤報不漏報成本遠高于一個幾十行就能實現(xiàn)的經典算法。1.2 為什么2025年了還要學這套老算法別看現(xiàn)在深度學習在信號處理里滿天飛Pan-Tompkins在真實產品里依然大量存在原因很樸素計算量低到可以忽略。整套流程每樣本約幾十次乘加運算在主流MCU上跑200Hz采樣綽綽有余中斷里順手就做了。行為可預測。它是純因果系統(tǒng)每一個樣本進來都能立刻給出中間結果不像神經網(wǎng)絡那樣存在這一幀到底看到多長上下文的模糊性。可解釋性好。哪個環(huán)節(jié)出了問題把帶通輸出、微分輸出、積分輸出逐個畫出來就能定位調參是透明的。不用訓練數(shù)據(jù)。換一個病人、換一種導聯(lián)自適應閾值會自動調整不需要重新訓練模型。當然它也有天花板對嚴重心律失常、強噪聲、胎兒心電這類場景Pan-Tompkins會力不從心。但即便用深度學習方法很多人也會先用Pan-Tompkins做候選檢測再用網(wǎng)絡做精細分類。作為baseline它永遠值得先跑一版。2. 算法全流程拆解五個環(huán)節(jié)各司其職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的陡峭沿削平后面微分環(huán)節(jié)就找不到明顯的斜率峰了。論文里的帶通不是直接設計一個帶通而是用兩個IIR濾波器級聯(lián)先低通再高通。低通濾波器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這兩個式子值得多看兩眼它們有個共同特點系數(shù)全是整數(shù)或2的冪次整條濾波鏈可以完全用整數(shù)加減和移位實現(xiàn)不需要浮點運算。這在1985年的8085處理器上是硬約束在今天做固件定點化同樣是巨大優(yōu)勢。高通那一路實際是全通減低通的結構用x[n-16]減去一個32點的滑動平均實現(xiàn)了5Hz高通同時保持了因果性和整數(shù)運算。2.2 第二三關求導與平方把QRS的個性放大帶通濾波之后QRS依然是低頻成分居多直接設閾值還是容易受到殘留噪聲干擾。這時候就要抓住QRS最獨特的個性斜率大。微分環(huán)節(jié)用的不是最樸素的差分而是這個式子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波留下的殘余。平方環(huán)節(jié)更直白y[n] x[n]^2。它有兩個作用一是把微分后正負脈沖都變成正值方便后面累加二是非線性放大讓大幅度成分和小幅度成分的差距進一步拉開。打個比方1和2相差一倍平方后變成1和4相差四倍QRS的微分峰值和T波微分峰值本來差距就不小平方之后這個差距被幾何級放大閾值判決就容易多了。2.3 第四關滑動窗口積分把多峰合并成一個峰微分加平方之后的信號是什么樣QRS的一個上升沿對應一個窄尖峰一個下降沿又對應一個窄尖峰整段信號是毛刺感很強的多峰形態(tài)。如果直接對它設閾值一個QRS很可能被判出兩三個結果。滑動窗口積分解決的就是這個問題。它把過去N個樣本的平方值取平均y[n] (1/N) * (x[n] x[n-1] ... x[n-N1])窗口取多長有講究。論文在200Hz采樣下取N30也就是150ms。這個窗口長度大約是QRS寬度的量級能正好把QRS的兩個斜率峰合并成一個平滑的單峰如果窗口太短合并效果差依然會多峰如果太長會把后面的T波也卷進來或者讓輸出波形變鈍時間分辨率下降。到這一步信號已經變成一串干凈的單峰序列峰的高低基本反映了這里有沒有一個高斜率、大幅度的QRS。剩下的事情就是怎么自動決定多高算一個峰。2.4 第五關自適應閾值與決策規(guī)則干的是判案的活固定閾值在實驗室數(shù)據(jù)上看著不錯一到真實場景就翻車電極貼得松緊不一樣、皮膚阻抗變化、病人深呼吸造成基線漂移都會讓QRS幅度在幾分鐘內大幅波動。Pan-Tompkins的精髓在于閾值是自適應的。算法維護兩個估計值信號峰值SPKsignal peak和噪聲峰值NPKnoise peak。每當檢測到一個峰PEAK就用指數(shù)滑動平均更新其中一個判定為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在剛才的緩沖數(shù)據(jù)里回頭找防止漏掉一個幅度突然變小的真實心跳。決策規(guī)則里還有兩個關鍵防錯機制。一個是200ms不應期QRS之后200ms內不允許再報一個峰因為正常心臟在這么短時間里不可能再興奮一次這一條能擋掉大部分高尖T波帶來的雙峰誤檢。另一個是T波斜率判別如果兩個檢測峰間隔小于360ms且第二個峰的斜率不足前一個QRS斜率的一半就認為第二個是T波按NPK更新處理。所有這些繞來繞去的規(guī)則核心就一句話在沒有人工干預的前提下讓檢測器能跟著信號質量自動調整嚴苛程度。SPK和NPK本質是兩個不斷被新樣本校準的錨點TH1和TH2是錨點之間的分界線。這種設計到今天依然是自適應檢測器的主流范式。3. 動手實現(xiàn)從公式到可跑通的代碼理論說得再漂亮不如把代碼跑一遍。這里給出一份完整的Python實現(xiàn)輸入是一段ECG信號numpy數(shù)組輸出是檢測到的QRS位置索引。代碼刻意保持了逐樣本處理的思路方便你改成實時流式版本。3.1 用差分方程實現(xiàn)濾波器鏈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): 滑動窗口積分流式實現(xiàn) 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這里每個函數(shù)我都按樣本循環(huán)寫而不是直接用scipy的濾波器。原因有兩個一是差分方程本身就是流式的改成環(huán)形緩沖區(qū)后可以直接塞進嵌入式中斷服務函數(shù)二是每一步中間結果都能取出來畫圖調參時一目了然。低通濾波器的延遲是6個樣本高通是16個樣本微分是2個樣本積分窗口中心約15個樣本整條鏈的固定延遲大約在200~300ms量級對實時監(jiān)護來說完全可接受。3.2 自適應閾值判決的主循環(huán)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: # 信號還在漲持續(xù)更新峰的位置和值 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這個主循環(huán)我刻意做了簡化把論文里復雜的T波斜率判別和雙閾值狀態(tài)機收斂成不應期主閾值回搜閾值三件套。實測下來對常規(guī)的竇性心律數(shù)據(jù)這個簡化版本在MIT-BIH部分記錄上已經能跑到95%以上的準確率。如果想要更逼近論文原版的99%以上表現(xiàn)需要把T波斜率判別和RR間期細分規(guī)則補回去這部分在第4節(jié)會展開講。調用方式很簡單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的數(shù)據(jù)可以裝一個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。原論文就是這么做的很多復現(xiàn)里省略了這一步但對HRV分析這種對時間精度敏感的用途建議還是補上。3.3 換采樣率時參數(shù)怎么改上面所有系數(shù)都是針對200Hz采樣率設計的。如果你的硬件是250Hz、360Hz或者500Hz采樣不要直接套公式要按比例縮放三個東西參數(shù)200Hz默認值換算方法低通延遲系數(shù)6和12乘以 fs/200 后取整高通延遲系數(shù)16和32乘以 fs/200 后取整積分窗口N30fs * 0.150 取整不應期40200msfs * 0.200 取整注意微分濾波器的延遲系數(shù)我建議保持不變因為它的本質是每相鄰幾個樣本做一次斜率估計在高采樣率下時間跨度更短反而能捕捉更陡的斜率效果不會變差。帶通的延遲系數(shù)如果不縮放截止頻率會隨著采樣率漂移比如500Hz下還用6和12低通截止就跑到近30Hz肌電干擾全進來了。4. 真機實測的坑與排查方法代碼能跑通只是第一步。我自己在真實心電數(shù)據(jù)上調試時踩過的坑比看論文時想象的多得多。這里把最有代表性的幾個問題列出來附帶排查思路。4.1 基線漂移導致的假陽性現(xiàn)象是檢測結果里突然冒出一串密集的峰尤其在病人深呼吸、電極線晃動的時候。表面上看帶通濾波已經處理了基線漂移但高通濾波器的記憶長度是32個樣本160ms對0.5Hz以下的極低頻成分抑制并不徹底。遇到大幅度緩慢漂移高通輸出會殘留一個相對陡的臺階這個臺階經過微分和平方后被放大足以騙過閾值。排查方法把帶通濾波后的信號畫出來如果能看到類似斜坡突變的波形基本就是基線漂移泄漏。解決思路有三個按優(yōu)先級排序先檢查電極接觸和導聯(lián)線固定這個最管用其次在帶通前加一個中值濾波或高通預處理把極低頻先壓下去最后實在不行把高通濾波器的延遲系數(shù)按比例加大讓截止頻率從5Hz稍微降到4Hz左右代價是P波和T波信息受損但對單純QRS檢測影響可控。4.2 高尖T波造成的雙峰誤檢正常T波幅度只有QRS的1/4左右但有些病人比如高鉀血癥、心肌缺血早期T波會變得又高又尖經過微分和平方后幅度直逼QRS。如果不做處理一個心動周期會被檢測出兩次心率直接翻倍。論文里的標準解法是T波斜率判別檢測到兩個峰間隔小于360ms時計算第二個峰的斜率如果它小于前一個QRS斜率的一半就把第二個峰當作T波丟棄。我實現(xiàn)的簡化版本里沒有加這個邏輯所以如果你在T波高尖的數(shù)據(jù)上測試出現(xiàn)雙峰誤檢是正常的。實測中還有個更簡單的輔助手段在判決后加一個最小間隔過濾凡是兩個檢測峰間隔小于250ms的保留幅度大的那個。這個規(guī)則犧牲了一點理論嚴謹性但實現(xiàn)成本極低在嵌入式上特別好用。4.3 閾值自適應跟不上信號突變自適應閾值最怕的不是噪聲而是信號本身的劇烈變化。比如病人翻了個身電極位置輕微移動QRS幅度在幾秒內掉到原來的一半。這時候SPK還停留在高位TH1遲遲降不下來結果就是連續(xù)漏檢。反過來如果某一段噪聲被誤判為QRSSPK被帶高也會引發(fā)后續(xù)漏檢。我在代碼里做了兩個保險一是RR間期的滑動平均參與回搜判決漏檢超過1.5倍平均RR時自動用低閾值搜索二是NPK更新時加了峰高不能超過當前SPK的鉗制防止把病態(tài)大噪聲學進去。如果你自己復現(xiàn)建議也加上一個防呆SPK和NPK都設置上下限比如NPK不低于SPK的5%SPK不高于歷史均值的數(shù)倍可以避免閾值完全失控。4.4 常見問題速查表現(xiàn)象可能原因排查與處理連續(xù)密集誤檢基線漂移、電極松動檢查導聯(lián)帶通前加極低頻抑制一個心跳報兩次T波高尖、積分窗口太短加250ms最小間隔或還原T波斜率判別突然漏檢幾十秒QRS幅度驟降、SPK未跟上依賴回搜邏輯RR間期防漏檢鉗制閾值首2秒檢測異常初始學習段含大偽差延長學習段或手動設置參考閾值采樣率換了效果變差延遲系數(shù)未縮放按fs/200重算所有延遲參數(shù)嵌入式上跑出NaN用了浮點且未初始化狀態(tài)濾波器全部整數(shù)化環(huán)形緩沖區(qū)先清零5. 關于評估和工程落地的一些經驗5.1 怎么科學評價一個QRS檢測器很多人跑完檢測器用眼睛看一眼覺得大概差不多就收工了這在論文和產品里都站不住腳。標準做法是用兩個指標敏感性Sensitivity和陽性預測值Positive Predictive Value。敏感性 TP / (TP FN)衡量的是真實心跳里找回了多少。漏檢越少敏感性越高。陽性預測值 TP / (TP FP)衡量的是報出來的結果里有多少是真的。誤檢越少陽性預測值越高。兩個指標要一起看因為你可以把閾值調到極低來刷敏感性但誤報會爆炸也可以把閾值調到極高來保證全對但漏檢會爆炸。只有兩個指標同時高才算真正好的檢測器。評估時的對齊規(guī)則也很重要一般以標注的真值位置為中心允許±150ms的誤差窗口檢測點落在窗口內就算TP。測試數(shù)據(jù)首選MIT-BIH Arrhythmia Database它有48條半小時的心電記錄和逐拍的專家標注是這一領域的事實標準。論文和后續(xù)改進算法在MIT-BIH上的成績普遍在敏感性99%以上、陽性預測值99%左右你復現(xiàn)時可以把這組數(shù)字當作及格線。5.2 嵌入式實時實現(xiàn)的三點建議如果要把這套算法移植到MCU上我有三個實際心得第一全部用整數(shù)運算。低通和高通的系數(shù)都是整數(shù)或2的冪微分里的除法是/8平方是整型乘法積分是累加和右移沒有任何一處需要浮點。定點化之后一個200Hz的通道在幾十MHz的MCU上占用CPU不到5%。當年論文在8085上都能跑今天的硬件完全是降維打擊。第二用環(huán)形緩沖區(qū)管理濾波器的歷史樣本。IIR濾波器需要訪問x[n-6]、x[n-16]、x[n-32]這些歷史值最自然的方式是開一個長度為延遲量的環(huán)形緩沖每次寫入新樣本讀指針跟著走。注意延遲量和環(huán)形緩沖長度要嚴格匹配頭一兩秒的初始狀態(tài)必須是全零否則濾波器的建立期會輸出一段完全錯誤的數(shù)據(jù)。第三把判決邏輯放在中斷里做時要控制每個樣本的處理時間上界。濾波鏈每樣本約30次整數(shù)運算加一次簡單狀態(tài)機跳轉在常見MCU上都在微秒級完全沒問題。關鍵是把閾值更新里涉及的開方、除法這類昂貴運算全部去掉論文里本來也用不到這些。5.3 什么時候別硬用Pan-Tompkins坦誠地說有三類場景我建議直接放棄Pan-Tompkins換更重的算法。一類是嚴重心律失常。比如房顫時RR間期完全無規(guī)律T波和下一個P波離得很近閾值自適應會頻繁振蕩檢測器會變得很不穩(wěn)定。另一類是胎兒心電這類極低信噪比任務母體心電是胎兒心電的好幾倍帶通濾波根本分不開。還有一類是強運動場景下的可穿戴設備步頻和運動偽跡的能量和QRS重疊嚴重單純靠頻域濾波已經壓不住。這些場景現(xiàn)代做法是上小波變換、模板匹配或者輕量級神經網(wǎng)絡。但即便如此Pan-Tompkins依然有價值的它適合做前端候選檢測先粗篩出可能是QRS的位置再用更重的算法做精細分類這樣能把深度模型的推理頻次降一個數(shù)量級。這套算法我現(xiàn)在還在用。前陣子做一個低成本心電心率監(jiān)測模塊客戶要求把算法放進一個小到沒有操作系統(tǒng)的MCU里我第一反應就是翻出Pan-Tompkins整數(shù)化之后總共不到兩百行C代碼跑在200Hz采樣上穩(wěn)得很。每次用都還是會感嘆一個1985年的設計結構清晰到每一步都能拆開調試也正因如此它的每一個參數(shù)都帶著能夠被理解的為什么。對于剛入行的人來說這是比任何現(xiàn)代黑盒模型都更好的第一課。