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