簡介一套用Python實現的心電算法工程面向生物醫(yī)學工程學習者、算法入門者與醫(yī)療數據分析人員解決心電信號去噪、R波定位、心率計算及心律失常識別等核心問題。代碼涵蓋巴特沃茲與卡爾曼濾波器、小波變換R波檢測、R-R間期心率計算并集成SVM、隨機森林和神經網絡分類器可完成房顫、室顫室速等病理狀態(tài)識別與偽差干擾研究。包內共109個文件以py源碼、dat/hea/atr心電數據文件為主要組成輔以pdf說明文檔、可視化圖表和測試用例壓縮包大小34.06MB。已有96人學習下載適合結合MIT-BIH數據庫驗證完整算法流程也可作為課程項目或畢業(yè)設計的參考實現。1. 這份 Python 心電資源包里裝了什么從濾波到房顫識別的完整鏈路看到資源列表里那串200.atr、203.atr、105.atr老讀者應該直接反應過來了——這是 MIT-BIH 心律失常數據庫的注釋文件。我拆過不少 Python 心電算法開源工程這套資源的價值不在于某個單點算法多炫而在于它把濾波、R 波檢測、心率計算、特征提取、心律失常分類到房顫/室顫識別這條鏈路完整打通了。適合三類人剛接觸生物醫(yī)學信號處理的學生、要做心電預處理和搏動檢測的算法工程師、以及想拿公開數據驗證自己分類模型的從業(yè)者。它能讓你少走至少兩周彎路因為每一環(huán)都有可運行代碼和對應的驗證文件。2. 濾波不是玄學給 ECG 信號做零相位帶通與去基線漂移2.1 為什么不用滑動平均而是巴特沃斯 filtfilt心電信號的主要噪聲來源就三類基線漂移頻率通常低于 0.5Hz、工頻干擾50Hz國內電網、肌電偽差高頻能量集中在 100Hz 以上。很多人一上來就套滑動窗口濾波那東西做平滑可以但會直接把 QRS 波群邊緣抹掉后面 R 波檢測的定位精度會肉眼可見地變差。心電處理里最穩(wěn)妥的方案是巴特沃斯帶通濾波器再用filtfilt做零相位濾波。filtfilt和普通lfilter的關鍵差別在于lfilter會引入與濾波器階數相關的群延遲相當于把整段信號在時間軸上平移了一段而filtfilt先正向過一遍再反向過一遍相位畸變互相抵消R 波的位置不會因為濾波而產生系統(tǒng)性偏移。在心率變異性分析這種對時間精度敏感的環(huán)節(jié)這是硬性要求。import numpy as np from scipy.signal import butter, filtfilt, iirnotch def bandpass_ecg(ecg, fs360.0, lowcut0.5, highcut50.0, order4): 零相位帶通濾波器默認按 MIT-BIH 的 360Hz 采樣率設置 nyq fs / 2.0 b, a butter(order, [lowcut / nyq, highcut / nyq], btypeband) return filtfilt(b, a, ecg) def notch_50hz(ecg, fs360.0, f050.0, q30): 50Hz 工頻陷波q 值越大陷波帶寬越窄 b, a iirnotch(f0, q, fs) return filtfilt(b, a, ecg)帶通范圍選 0.5–50Hz 的依據0.5Hz 以上的范圍能保留大部分 ST 段和 T 波形態(tài)同時壓掉呼吸引起的基線漂移50Hz 以上的信號里 QRS 高頻分量已經很少切掉不影響 R 波定位反而能削弱肌電噪聲。如果你的數據是 250Hz 采樣把highcut改到 40Hz 更穩(wěn)妥因為 50Hz 已經在奈奎斯特頻率附近濾波器邊緣會變得很陡容易出現振鈴。2.2 去除基線漂移的另一種思路中值濾波與移動平均巴特沃斯高通只能做全局去漂移如果碰到電極接觸不良導致的大幅低頻擺動帶通濾波之后可能依然殘留一段弧線。常見的補充做法是用中值濾波估計基線再從原信號里減掉。中值濾波的優(yōu)勢是它不假設基線符合某個固定頻率范圍對突變型漂移更魯棒。from scipy.ndimage import median_filter def remove_baseline(ecg, fs360.0, window_ms200): 用中值濾波估計基線并扣除 window_ms 控制窗口大小200ms 對 QRS 太窄不會把波形吃掉 win int(fs * window_ms / 1000) // 2 * 2 1 # 強制奇數窗 baseline median_filter(ecg, sizewin, modenearest) return ecg - baselinewindow_ms這個參數我一般取 150–250ms。窗口小于 QRS 寬度時中值濾波會把 QRS 本身也當成基線的一部分濾完波形就平了窗口太大又跟不上快速漂移。modenearest是為了避免邊緣效應采樣點前段和中段的濾波效果差異明顯這類細節(jié)就是調參中的“血淚經驗”。2.3 濾波結果怎么驗證頻響曲線與波形對照濾波不是跑完就完事的一定要做兩步驗證。第一步是畫頻響曲線確認通帶確實平坦、阻帶有足夠衰減import matplotlib.pyplot as plt from scipy.signal import sosfreqz b, a butter(4, [0.5 / 180.0, 50 / 180.0], btypeband) w, h sosfreqz([b], [a], worN2048, fs360) plt.semilogx(w, 20 * np.log10(abs(h)))第二步是直接把濾波前后的信號疊在一起畫。重點看兩件事R 波峰值是否被削平、濾波后的信號在起始段有沒有明顯上沖或下沖。filtfilt的邊緣效應雖然比lfilter小但當濾波器階數超過 5 時信號的頭尾幾十個采樣點仍然可能出現瞬時波動。我習慣先把信號前后各延拓 2 秒再濾波濾完切掉徹底規(guī)避邊緣問題。3. 抓住 R 波Pan-Tompkins 實現、閾值自適應與心率計算3.1 Pan-Tompkins 的四個階段與參數拆解MIT-BIH 注釋文件里的 R 波位置是人工逐搏校正過的算法上最貼近它精度的經典方法還是 Pan-Tompkins。它分四步帶通濾波5–15Hz、差分、平方、滑動窗口積分。帶通聚焦 QRS 的主要能量頻帶讓 P 波和 T 波大幅衰減差分突出 R 波的斜率特征平方讓所有值轉正并放大高頻分量滑動窗口積分把單個尖峰變成一個光滑的駝峰方便用閾值判斷。def pan_tompkins(ecg, fs360.0, refractory0.2, win_ms150, thresh_factor0.35): 完整 Pan-Tompkins R 波檢測 返回 R 峰位置采樣點下標和對應的濾波信號 b, a butter(4, [5 / (fs / 2), 15 / (fs / 2)], btypeband) filtered filtfilt(b, a, ecg) diff np.diff(filtered, prependfiltered[0]) squared diff ** 2 win int(fs * win_ms / 1000) window np.ones(win) / win integrated np.convolve(squared, window, modesame) # 主角自適應閾值 peak_thresh 0.3 * np.max(integrated[:int(fs * 2)]) r_peaks [] ref_period int(refractory * fs) i 0 while i len(integrated): if integrated[i] peak_thresh: # 以積分信號局部極大值為準回找原始濾波信號上的最大點 j min(i ref_period, len(integrated) - 1) seg integrated[i:j] loc i np.argmax(seg) r_loc i np.argmax(filtered[i:j]) r_peaks.append(r_loc) # 用最新峰值動態(tài)更新閾值自適應核心 peak_thresh 0.5 * peak_thresh 0.5 * integrated[loc] i j else: i 1 return np.array(r_peaks), filteredrefractory不應期取 0.2 秒——正常心率下兩次 R 波間隔不會小于 0.2 秒這段時間內即使有高幅噪聲也不會被當作新搏動。win_ms取 150ms 是經驗值它是 QRS 典型寬度 80–120ms 的 1.2–1.5 倍能讓積分窗口剛好覆蓋一個完整 QRS 駝峰。thresh_factor最初用前 2 秒信號最大值乘 0.3 作為初始閾值這段信號一般包含 2–4 個搏動能適應前 2 秒的幅值水平。這段代碼我做了個簡化處理初始閾值用固定比例等檢測到 5 個搏動后再進入“0.5 舊閾值 0.5 新峰值”的更新模式。這個策略在實際跑 MIT-BIH 的200.atr時表現不錯但跑203.atr這類有大量室早的記錄時單閾值跟不上形態(tài)突變第二年就有論文提出雙閾值方案。后面避坑章節(jié)會展開。3.2 初始閾值失敗怎么辦雙閾值與回掃機制Pan-Tompkins 最經典的翻車場景是前 2 秒信號里有個巨大偽差峰值閾值被抬得太高導致后續(xù)真正的 R 波全部漏檢或者前 2 秒全是低幅信號閾值太低后面高頻噪聲全被當成 R 波。處理辦法是雙閾值回掃高閾值檢測出保真的 R 波集合對高閾值漏檢的片段用低閾值通常取高閾值的 0.5 倍重新掃一遍看能否找到時間位置合理的疑似峰。def detect_with_dual_threshold(integrated, filtered, fs, thresh_high0.3, thresh_low0.15, refractory0.2, methodpeak): 高閾值初篩 低閾值回掃緩解 R 波幅值突變導致的漏檢 thresh_high/thresh_low 可以傳絕對值也可以傳比例因子 rr_min int(refractory * fs) # 第一階段高閾值 high_peaks find_peaks(integrated, heightthresh_high, distancerr_min) # 第二階段每個相鄰高閾值峰之間的間隙用低閾值再掃 final_peaks list(high_peaks[0]) for k in range(len(final_peaks) - 1): gap_start final_peaks[k] rr_min gap_end final_peaks[k 1] - rr_min if gap_end - gap_start 0: continue seg integrated[gap_start:gap_end] local_max np.max(seg) if local_max thresh_low: loc_rel np.argmax(seg) final_peaks.append(gap_start loc_rel) final_peaks.sort() return final_peaks回掃不是無腦加峰兩個硬約束必須同時滿足候選峰與前后已確認峰的距離必須大于不應期回掃峰對應的原始濾波信號幅值至少要達到前后 R 波平均幅值的三分之一。203.atr里室早和正常搏動的形態(tài)差異極大高閾值抓到的是正常 QRS室早的積分幅值可能只有正常的一半這時候低閾值回掃就是救場的那個角色。3.3 心率計算的正確姿勢RR 間期異常值過濾R 波位置拿到之后心率計算本身不復雜但很容易算錯。直接對所有 RR 間期取平均是一種常見錯誤——因為漏檢和誤檢會制造異常 RR 間期比如一次漏檢會把兩個搏動合并成一個 2 秒的間隔一次誤檢會制造一個 0.1 秒的尖峰。正確的流程是先算全部 RR 間期把小于 0.4 秒或大于 1.5 秒的間期過濾掉再用剩余間期計算平均心率。import pandas as pd def compute_hr(r_peaks, fs360.0): 基于 R 峰位置計算平均心率和 RR 間期序列 rr np.diff(r_peaks) / fs # 單位秒 rr rr[(rr 0.4) (rr 1.5)] # 濾掉異常間隔 hr 60.0 / np.mean(rr) rr_df pd.DataFrame({rr_time: np.cumsum(rr), rr_value: rr}) return hr, rr_df濾波范圍 0.4–1.5 秒對應心率 40–150 次/分鐘。健康成年人安靜狀態(tài)下心率基本在這個范圍里如果實測心率超過上限先把 R 波檢測的漏檢率檢查一遍——多數“爆表”的假心率都是漏檢導致的。這里pandas的使用不只是錦上添花后續(xù)做心率變異性HRV分析時RR 間期序列的時間對齊、缺失值填充、滾動統(tǒng)計都需要 DataFrame 的接口。4. 特征提取與心律失常分類從 QRS 形態(tài)到隨機森林4.1 提取哪些特征QRS 幅度、寬度、QT 間期與 RR 間期搏動級分類和心拍級分類用的特征集不同。心拍級分類區(qū)分正常、室早、房早關心的是單個 QRS 的形態(tài)特征心律失常類型判斷房顫 vs 竇性關心的是 RR 間期的動態(tài)特征。這套資源的核心是前者但它的特征提取函數是通用的完全可以直接復用。常見做法是以每個 R 峰為中心往前取 100ms 估 Q 波起點往后取 300ms 估 T 波終點然后計算下面這組特征特征名計算方式臨床意義PeakAmpR 峰幅值 - Q 波谷幅值QRS 整體振幅QRSDurationS 波終點 - Q 波起點正常 80–120ms增寬提示束支阻滯QTIntervalT 波終點 - Q 波起點QT 延長與惡性心律失常相關PRIntervalR 峰前 120ms 處到 Q 波起點房室傳導時間RRPrev前一個 RR 間期長度識別代償間歇RRAfter后一個 RR 間期長度識別早搏后停頓T_AmpT 波峰值T 波倒置/高尖提示缺血def extract_beat_features(ecg, r_peaks, fs360.0): 以 R 峰為錨點在固定窗內提取 QRS 形態(tài)特征 features [] q_win int(0.08 * fs) # Q 波搜索窗R 前 80ms s_win int(0.12 * fs) # S 波搜索窗R 后 120ms for r in r_peaks: if r - q_win 0 or r s_win len(ecg): continue q_min np.min(ecg[r - q_win:r]) s_min np.min(ecg[r:r s_win]) r_val ecg[r] peak_amp r_val - q_min qrs_dur (np.argmin(ecg[r:r s_win]) np.argmin(ecg[r - q_win:r])) / fs rr_prev (r_peaks[np.where(r_peaks r)[0][-1]] if any(r_peaks r) else r) - r if any(r_peaks r) else 0 features.append([peak_amp, qrs_dur, rr_prev, r_val, q_min, s_min]) return np.array(features)這些特征直接喂給分類器不夠還要做標準化。QRS 幅值跨度可以從 0.5mV 到 3mV而 PR 間期以毫秒為單位數值量級差幾百倍。sklearn的StandardScaler會在內部把各特征拉到同一量綱這一步千萬別省否則 SVM 的核函數計算會被大數值特征主導小特征等于白提。4.2 SVM 與隨機森林的完整訓練流程搏動分類最穩(wěn)的組合是上述 6–8 維特征 隨機森林。隨機森林對特征量綱不敏感、對缺失值有容忍度而且能輸出特征重要性排序方便你反查哪些特征對分類貢獻最大。SVM 在小樣本高維場景下表現出色但需要調核函數和 C 參數如果只是為了快速拿到一個可用的分類器隨機森林是第一選擇。from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler from sklearn.metrics import classification_report # X: (n_samples, n_features), y: 搏動標簽, 0正常, 1室早, 2房早 X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.3, stratifyy, random_state42 ) scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train) X_test_scaled scaler.transform(X_test) clf RandomForestClassifier( n_estimators300, max_depth10, class_weightbalanced, # 處理類別不均衡 random_state42 ) clf.fit(X_train_scaled, y_train) y_pred clf.predict(X_test_scaled) print(classification_report(y_test, y_pred))class_weightbalanced是必選項。MIT-BIH 的118系列記錄里正常搏動可能占 90% 以上如果不加權分類器只要把所有樣本都預測為正常類就能拿到 90% 準確率但室早一個都發(fā)現不了。stratifyy保證訓練集和測試集里各類別比例與原數據一致避免隨機劃分時測試集恰好沒有某類樣本。4.3 房顫與室顫識別從 RR 間期變異到特征統(tǒng)計房顫的核心特征是 RR 間期高度不規(guī)律——絕對不齊。這個特性決定了它不需要太復雜的模型用 RR 間期序列的統(tǒng)計量就能做到相當靠譜的初篩。我習慣提取三個指標RR 間期標準差SDNN、連續(xù) RR 差值的均方根RMSSD、以及 Poincaré 散點圖的寬度——正常竇性心律的 Poincaré 圖呈彗星狀房顫則呈寬散形。def af_detection_features(rr_array): 房顫的 RR 間期特征SDNN、RMSSD、pNN50 rr np.asarray(rr_array, dtypefloat) sdnn np.std(rr) diff_rr np.abs(np.diff(rr)) rmssd np.sqrt(np.mean(diff_rr ** 2)) pnn50 np.mean(diff_rr 0.05) * 100 # 相鄰 RR 差超過 50ms 的占比 return np.array([sdnn, rmssd, pnn50])這三個特征喂給邏輯回歸或隨機森林對 MIT-BIH 里的房顫記錄比如105.atr附近的多條數據通常能拿到 0.85 以上的 AUC。室顫則是另一個模式信號退化成完全不規(guī)則的細碎波形R 波檢測基本失效這時候看的不再是 RR 間期而是信號的幅度概率密度和頻率分布。真遇到這個狀態(tài)算法層面已經沒必要繼續(xù)走搏動分類流程當務之急是觸發(fā)惡性心律失常報警——這是個完全不同的檢測分支資源包里把它單列出來是有道理的。5. 避坑錦囊六個會翻車的 ECG 細節(jié)與排查手段5.1 R 波檢測結果里出現大量誤檢現象檢測出的 R 峰數量明顯多于實際搏動數量且位置集中在 T 波或噪聲段。原因T 波在某些導聯(lián)上幅值不低經過 5–15Hz 帶通后仍有殘留更常見的是滑動窗口積分閾值設置過低T 波被當成候選峰。還有一種隱蔽情況是refractory寫成了秒數忘記換算成采樣點數。解決先手動畫圖把 R 峰標注疊加在濾波后的信號上看誤檢峰與正常 R 波之間的間隔。如果誤檢峰與前面真峰的間隔小于 0.25 秒說明refractory失效或沒有生效如果是規(guī)律性出現在每個 T 波位置把帶通低端從 5Hz 提到 8Hz衰減低斜率和 T 波的低頻分量。做完這兩步誤檢通常能降一個數量級。5.2 濾波后 QRS 形態(tài)整體變形現象濾波前后的 QRS 波形對不上R 波寬度變大或出現雙峰。原因帶通濾波器的order設置過高。階數超過 6 時濾波器在截止頻率附近會有嚴重的相位畸變雖然filtfilt補償了相位延遲但幅頻響應的過沖會導致 QRS 邊緣出現振鈴偽跡。解決把order降回 4優(yōu)先保證波形形態(tài)完整性。如果對高頻噪聲的抑制不夠用級聯(lián)的方式先 4 階低通 50Hz再 4 階高通 0.5Hz兩級串聯(lián)對信號的相位影響比單級 8 階小得多。5.3 心率計算結果忽高忽低平均心率卻正?，F象瞬時心率曲線劇烈跳變但平均心率在合理范圍。原因漏檢和誤檢同時存在。漏檢把兩個 RR 間期合并成一個長間隔拉低瞬時心率誤檢又在里面插了一個短間隔把瞬時心率抬高。平均后正負抵消掩蓋了問題。解決瞬時心率曲線必須經過中值濾波再展示。用長度為 5 的中值窗口能干凈利落地去掉單點異常同時回到 R 波檢測環(huán)節(jié)把漏檢的片段找出來單獨調閾值。我在實際項目里發(fā)現很多所謂的“心率算法不穩(wěn)定”問題根源不在心率計算而在 R 波檢測。5.4 ATR 文件讀進來后與自己的檢測結果對不上現象用 WFDB 工具讀200.atr得到的 R 波位置和自己檢測出的位置總是差 5–10 個采樣點。原因MIT-BIH 的注釋位置標注的是 QRS 波群的某個參考點通常是最大斜率點或峰值附近不同版本的標注工具參考點定義不完全一致另一個常見原因是自己的信號預處理過程中用了np.diff沒有注意到diff會使信號長度減 1后續(xù)所有索引整體偏移。解決先確認采樣率是否一致——MIT-BIH 是 360Hz如果你用自己的采樣率需要resample_poly做重采樣。然后允許檢測結果與注釋之間有一個 ±15ms約 5 個采樣點的誤差窗口在這個窗口內都記為正確檢出。臨床上評價 QRS 檢測器的標準是靈敏度和陽性預測值本身就容忍 10ms 級的偏差。5.5 特征提取時切片越界程序崩潰現象提取特征時IndexError集中在信號尾部。原因最后一個 R 峰距離信號末尾不足一個特征窗口長度向后取 S 波或 T 波窗口時溢出信號開頭的 R 峰則可能因為向前取 Q 波窗口導致負索引。解決特征提取函數里必須加邊界保護最粗暴但有效的方式是直接丟棄前后不足一個窗口長度的 R 峰。欠采樣時丟失幾個搏動完全不影響分類器性能但越界崩潰會直接中斷整個流程。類似q_min np.min(ecg[r - q_win:r])這行代碼看起來簡單不加if r - q_win 0就是定時炸彈。5.6 訓練集分類準確率 95%測試集只有 60%現象模型在訓練集上表現優(yōu)異獨立測試集上性能明顯下降。原因特征提取時使用了整個數據集的信息。最常見的泄漏是特征標準化用了全數據的均值和方差或者打亂訓練集和測試集時沒有按記錄分組——同一條記錄里相鄰搏動的特征高度相關被分到訓練集和測試集兩側時模型相當于“見過”了同類樣本。解決按記錄分組劃分數據同一條記錄的所有搏動只能出現在訓練集或測試集中不能兩邊都出現。標注數據時也要注意是否存在時間上的重疊窗口。這條規(guī)則適用于心電、腦電、肌電等所有時間序列相關的分類任務。6. 用 ATR 文件做驗證把注釋解析成標簽并量化檢測精度這套資源里的.atr文件是 MIT-BIH 的心律注釋文件它記錄了每個搏動的類型代碼N 表示正常、V 表示室早、A 表示房早等和精確時間位置。對我們的價值是可以用它作為金標準量化自己 R 波檢測器的靈敏度和陽性預測值。import wfdb def load_atr_labels(record_name): 讀取 MIT-BIH 注釋文件返回搏動位置與類型標簽 ann wfdb.rdann(record_name, atr) beat_idx [i for i, s in enumerate(ann.symbol) if s in NLVARJES] beat_pos np.array(ann.sample)[beat_idx] beat_type np.array(ann.symbol)[beat_idx] return beat_pos, beat_type拿到金標準后匹配邏輯是這樣的對自己的每個檢測點在 ±15ms 容差范圍內尋找最近的金標準 R 波位置。找到則記為真陽性TP找不到則記為假陽性FP反過來金標準位置附近沒有自己的檢測點記為假陰性FN。靈敏度 TP/(TPFN) 和陽性預測值 TP/(TPFP) 都超過 0.95這算法才算合格。我的最后一步永遠是可視化驗證把 R 波峰和注釋位置同時畫在信號上用不同顏色區(qū)分真陽性、假陽性、漏檢。一套完整的心電分析流程從濾波到房顫識別中間任何一環(huán)出問題都會在圖上暴露出來。從那以后我每次跑新數據集都會強制走一遍這個流程——先把濾波前后波形疊畫檢查再跑 R 波檢測最后對標簽算指標缺一步都不安心。希望幫到你。本文還有配套的精品資源點擊獲取