简介一套用Python实现的心电算法工程面向生物医学工程学习者、算法入门者与医疗数据分析人员解决心电信号去噪、R波定位、心率计算及心律失常识别等核心问题。代码涵盖巴特沃兹与卡尔曼滤波器、小波变换R波检测、R-R间期心率计算并集成SVM、随机森林和神经网络分类器可完成房颤、室颤室速等病理状态识别与伪差干扰研究。包内共109个文件以py源码、dat/hea/atr心电数据文件为主要组成辅以pdf说明文档、可视化图表和测试用例压缩包大小34.06MB。已有96人学习下载适合结合MIT-BIH数据库验证完整算法流程也可作为课程项目或毕业设计的参考实现。1. 这份 Python 心电资源包里装了什么从滤波到房颤识别的完整链路看到资源列表里那串200.atr、203.atr、105.atr老读者应该直接反应过来了——这是 MIT-BIH 心律失常数据库的注释文件。我拆过不少 Python 心电算法开源工程这套资源的价值不在于某个单点算法多炫而在于它把滤波、R 波检测、心率计算、特征提取、心律失常分类到房颤/室颤识别这条链路完整打通了。适合三类人刚接触生物医学信号处理的学生、要做心电预处理和搏动检测的算法工程师、以及想拿公开数据验证自己分类模型的从业者。它能让你少走至少两周弯路因为每一环都有可运行代码和对应的验证文件。2. 滤波不是玄学给 ECG 信号做零相位带通与去基线漂移2.1 为什么不用滑动平均而是巴特沃斯 filtfilt心电信号的主要噪声来源就三类基线漂移频率通常低于 0.5Hz、工频干扰50Hz国内电网、肌电伪差高频能量集中在 100Hz 以上。很多人一上来就套滑动窗口滤波那东西做平滑可以但会直接把 QRS 波群边缘抹掉后面 R 波检测的定位精度会肉眼可见地变差。心电处理里最稳妥的方案是巴特沃斯带通滤波器再用filtfilt做零相位滤波。filtfilt和普通lfilter的关键差别在于lfilter会引入与滤波器阶数相关的群延迟相当于把整段信号在时间轴上平移了一段而filtfilt先正向过一遍再反向过一遍相位畸变互相抵消R 波的位置不会因为滤波而产生系统性偏移。在心率变异性分析这种对时间精度敏感的环节这是硬性要求。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 波形态同时压掉呼吸引起的基线漂移50Hz 以上的信号里 QRS 高频分量已经很少切掉不影响 R 波定位反而能削弱肌电噪声。如果你的数据是 250Hz 采样把highcut改到 40Hz 更稳妥因为 50Hz 已经在奈奎斯特频率附近滤波器边缘会变得很陡容易出现振铃。2.2 去除基线漂移的另一种思路中值滤波与移动平均巴特沃斯高通只能做全局去漂移如果碰到电极接触不良导致的大幅低频摆动带通滤波之后可能依然残留一段弧线。常见的补充做法是用中值滤波估计基线再从原信号里减掉。中值滤波的优势是它不假设基线符合某个固定频率范围对突变型漂移更鲁棒。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是为了避免边缘效应采样点前段和中段的滤波效果差异明显这类细节就是调参中的“血泪经验”。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 秒再滤波滤完切掉彻底规避边缘问题。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) # 用最新峰值动态更新阈值自适应核心 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这类有大量室早的记录时单阈值跟不上形态突变第二年就有论文提出双阈值方案。后面避坑章节会展开。3.2 初始阈值失败怎么办双阈值与回扫机制Pan-Tompkins 最经典的翻车场景是前 2 秒信号里有个巨大伪差峰值阈值被抬得太高导致后续真正的 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里室早和正常搏动的形态差异极大高阈值抓到的是正常 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 次/分钟。健康成年人安静状态下心率基本在这个范围里如果实测心率超过上限先把 R 波检测的漏检率检查一遍——多数“爆表”的假心率都是漏检导致的。这里pandas的使用不只是锦上添花后续做心率变异性HRV分析时RR 间期序列的时间对齐、缺失值填充、滚动统计都需要 DataFrame 的接口。4. 特征提取与心律失常分类从 QRS 形态到随机森林4.1 提取哪些特征QRS 幅度、宽度、QT 间期与 RR 间期搏动级分类和心拍级分类用的特征集不同。心拍级分类区分正常、室早、房早关心的是单个 QRS 的形态特征心律失常类型判断房颤 vs 窦性关心的是 RR 间期的动态特征。这套资源的核心是前者但它的特征提取函数是通用的完全可以直接复用。常见做法是以每个 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 形态特征 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 与随机森林的完整训练流程搏动分类最稳的组合是上述 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% 准确率但室早一个都发现不了。stratifyy保证训练集和测试集里各类别比例与原数据一致避免随机划分时测试集恰好没有某类样本。4.3 房颤与室颤识别从 RR 间期变异到特征统计房颤的核心特征是 RR 间期高度不规律——绝对不齐。这个特性决定了它不需要太复杂的模型用 RR 间期序列的统计量就能做到相当靠谱的初筛。我习惯提取三个指标RR 间期标准差SDNN、连续 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。室颤则是另一个模式信号退化成完全不规则的细碎波形R 波检测基本失效这时候看的不再是 RR 间期而是信号的幅度概率密度和频率分布。真遇到这个状态算法层面已经没必要继续走搏动分类流程当务之急是触发恶性心律失常报警——这是个完全不同的检测分支资源包里把它单列出来是有道理的。5. 避坑锦囊六个会翻车的 ECG 细节与排查手段5.1 R 波检测结果里出现大量误检现象检测出的 R 峰数量明显多于实际搏动数量且位置集中在 T 波或噪声段。原因T 波在某些导联上幅值不低经过 5–15Hz 带通后仍有残留更常见的是滑动窗口积分阈值设置过低T 波被当成候选峰。还有一种隐蔽情况是refractory写成了秒数忘记换算成采样点数。解决先手动画图把 R 峰标注叠加在滤波后的信号上看误检峰与正常 R 波之间的间隔。如果误检峰与前面真峰的间隔小于 0.25 秒说明refractory失效或没有生效如果是规律性出现在每个 T 波位置把带通低端从 5Hz 提到 8Hz衰减低斜率和 T 波的低频分量。做完这两步误检通常能降一个数量级。5.2 滤波后 QRS 形态整体变形现象滤波前后的 QRS 波形对不上R 波宽度变大或出现双峰。原因带通滤波器的order设置过高。阶数超过 6 时滤波器在截止频率附近会有严重的相位畸变虽然filtfilt补偿了相位延迟但幅频响应的过冲会导致 QRS 边缘出现振铃伪迹。解决把order降回 4优先保证波形形态完整性。如果对高频噪声的抑制不够用级联的方式先 4 阶低通 50Hz再 4 阶高通 0.5Hz两级串联对信号的相位影响比单级 8 阶小得多。5.3 心率计算结果忽高忽低平均心率却正常现象瞬时心率曲线剧烈跳变但平均心率在合理范围。原因漏检和误检同时存在。漏检把两个 RR 间期合并成一个长间隔拉低瞬时心率误检又在里面插了一个短间隔把瞬时心率抬高。平均后正负抵消掩盖了问题。解决瞬时心率曲线必须经过中值滤波再展示。用长度为 5 的中值窗口能干净利落地去掉单点异常同时回到 R 波检测环节把漏检的片段找出来单独调阈值。我在实际项目里发现很多所谓的“心率算法不稳定”问题根源不在心率计算而在 R 波检测。5.4 ATR 文件读进来后与自己的检测结果对不上现象用 WFDB 工具读200.atr得到的 R 波位置和自己检测出的位置总是差 5–10 个采样点。原因MIT-BIH 的注释位置标注的是 QRS 波群的某个参考点通常是最大斜率点或峰值附近不同版本的标注工具参考点定义不完全一致另一个常见原因是自己的信号预处理过程中用了np.diff没有注意到diff会使信号长度减 1后续所有索引整体偏移。解决先确认采样率是否一致——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%现象模型在训练集上表现优异独立测试集上性能明显下降。原因特征提取时使用了整个数据集的信息。最常见的泄漏是特征标准化用了全数据的均值和方差或者打乱训练集和测试集时没有按记录分组——同一条记录里相邻搏动的特征高度相关被分到训练集和测试集两侧时模型相当于“见过”了同类样本。解决按记录分组划分数据同一条记录的所有搏动只能出现在训练集或测试集中不能两边都出现。标注数据时也要注意是否存在时间上的重叠窗口。这条规则适用于心电、脑电、肌电等所有时间序列相关的分类任务。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 波峰和注释位置同时画在信号上用不同颜色区分真阳性、假阳性、漏检。一套完整的心电分析流程从滤波到房颤识别中间任何一环出问题都会在图上暴露出来。从那以后我每次跑新数据集都会强制走一遍这个流程——先把滤波前后波形叠画检查再跑 R 波检测最后对标签算指标缺一步都不安心。希望帮到你。本文还有配套的精品资源点击获取
