简介：一套用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, fs=360.0, lowcut=0.5, highcut=50.0, order=4): """ 零相位带通滤波器，默认按 MIT-BIH 的 360Hz 采样率设置 """ nyq = fs / 2.0 b, a = butter(order, [lowcut / nyq, highcut / nyq], btype='band') return filtfilt(b, a, ecg) def notch_50hz(ecg, fs=360.0, f0=50.0, q=30): """ 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, fs=360.0, window_ms=200): """ 用中值滤波估计基线并扣除 window_ms 控制窗口大小，200ms 对 QRS 太窄，不会把波形吃掉 """ win = int(fs * window_ms / 1000) // 2 * 2 + 1 # 强制奇数窗 baseline = median_filter(ecg, size=win, mode='nearest') return ecg - baselinewindow_ms这个参数我一般取 150–250ms。窗口小于 QRS 宽度时，中值滤波会把 QRS 本身也当成基线的一部分，滤完波形就平了；窗口太大又跟不上快速漂移。mode='nearest'是为了避免边缘效应，采样点前段和中段的滤波效果差异明显，这类细节就是调参中的“血泪经验”。
2.3 滤波结果怎么验证：频响曲线与波形对照
滤波不是跑完就完事的，一定要做两步验证。第一步是画频响曲线，确认通带确实平坦、阻带有足够衰减：
import matplotlib.pyplot as plt from scipy.signal import sosfreqz b, a = butter(4, [0.5 / 180.0, 50 / 180.0], btype='band') w, h = sosfreqz([b], [a], worN=2048, fs=360) 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, fs=360.0, refractory=0.2, win_ms=150, thresh_factor=0.35): """ 完整 Pan-Tompkins R 波检测 返回 R 峰位置（采样点下标）和对应的滤波信号 """ b, a = butter(4, [5 / (fs / 2), 15 / (fs / 2)], btype='band') filtered = filtfilt(b, a, ecg) diff = np.diff(filtered, prepend=filtered[0]) squared = diff ** 2 win = int(fs * win_ms / 1000) window = np.ones(win) / win integrated = np.convolve(squared, window, mode='same') # 主角：自适应阈值 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_high=0.3, thresh_low=0.15, refractory=0.2, method='peak'): """ 高阈值初筛 + 低阈值回扫，缓解 R 波幅值突变导致的漏检 thresh_high/thresh_low 可以传绝对值，也可以传比例因子 """ rr_min = int(refractory * fs) # 第一阶段：高阈值 high_peaks = find_peaks(integrated, height=thresh_high, distance=rr_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, fs=360.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 波终点，然后计算下面这组特征：
| 特征名 | 计算方式 | 临床意义 |
|---|---|---|
| PeakAmp | R 峰幅值 - Q 波谷幅值 | QRS 整体振幅 |
| QRSDuration | S 波终点 - Q 波起点 | 正常 80–120ms，增宽提示束支阻滞 |
| QTInterval | T 波终点 - Q 波起点 | QT 延长与恶性心律失常相关 |
| PRInterval | R 峰前 120ms 处到 Q 波起点 | 房室传导时间 |
| RRPrev | 前一个 RR 间期长度 | 识别代偿间歇 |
| RRAfter | 后一个 RR 间期长度 | 识别早搏后停顿 |
| T_Amp | T 波峰值 | T 波倒置/高尖提示缺血 |
def extract_beat_features(ecg, r_peaks, fs=360.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_size=0.3, stratify=y, random_state=42 ) scaler = StandardScaler() X_train_scaled = scaler.fit_transform(X_train) X_test_scaled = scaler.transform(X_test) clf = RandomForestClassifier( n_estimators=300, max_depth=10, class_weight='balanced', # 处理类别不均衡 random_state=42 ) clf.fit(X_train_scaled, y_train) y_pred = clf.predict(X_test_scaled) print(classification_report(y_test, y_pred))class_weight='balanced'是必选项。MIT-BIH 的118系列记录里，正常搏动可能占 90% 以上，如果不加权，分类器只要把所有样本都预测为正常类就能拿到 90% 准确率，但室早一个都发现不了。stratify=y保证训练集和测试集里各类别比例与原数据一致，避免随机划分时测试集恰好没有某类样本。
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, dtype=float) 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/(TP+FN) 和阳性预测值 TP/(TP+FP) 都超过 0.95，这算法才算合格。
我的最后一步永远是可视化验证：把 R 波峰和注释位置同时画在信号上，用不同颜色区分真阳性、假阳性、漏检。一套完整的心电分析流程，从滤波到房颤识别，中间任何一环出问题都会在图上暴露出来。从那以后我每次跑新数据集，都会强制走一遍这个流程——先把滤波前后波形叠画检查，再跑 R 波检测，最后对标签算指标，缺一步都不安心。希望帮到你。
本文还有配套的精品资源，点击获取
