Python心电信号处理工程：滤波、R波识别与心率计算流水线

简介面向心电信号处理研究者的Python算法工程系统覆盖滤波去噪、R波检测、心率计算、特征提取、心率失常分类以及房颤、室颤室速等病理状态识别并包含伪差干扰研究与可视化模块配套测试工程便于验证。项目按功能拆分代码清晰展示从原始ECG数据读取、scipy.signal滤波设计、小波/短时傅立叶变换定位R波到基于R-R间期计算心率、构建分类模型识别异常心律的完整链路能帮助医疗AI或生物医学工程开发者快速理解心电分析算法落地方案。资源共91个文件压缩包约30.81MB主要包含Python源码17个py、MIT-BIH心电数据dat/hea/atr格式、PDF论文与Markdown笔记、结果展示图片及Jupyter Notebook等目录结构规整便于按模块检索与二次开发。已有673人下载学习适合想复现心电算法流程、研究心律失常分类或作为课程设计参考的读者。1. 把滤波、R波识别、心率计算串成流水线的 Python 心电算法工程拿到一段心电图数据第一眼看到的往往不是心动而是满屏的干扰基线上下漂工频叠在上面肌电毛刺把 R 波边缘搅得乱七八糟。滤波参数调错一档后面 R 波识别、心率计算、心律失常分类全部跟着错位。这套工程把这些环节串成了一条可以直接跑通的处理链覆盖滤波、R 波识别、心率计算、特征提取以及房颤、室颤室速等异常心律的识别流程还带了可视化模块和测试工程。适合刚接触心电信号处理的新手做整体复现也适合已经写过简单算法的人对照目录快速定位到自己缺的那一环。它解决的核心问题是让原始心电数据从波形走向可解释的指标而不是卡在“R 波到底在哪”的手工找峰里。2. 数据层优先mit-bih-arrhythmia、mit-bih-noise、aha、paf 四份目录的读取与标注边界解压工程后先看到的不是算法模块而是一排数据目录。很多人习惯性跳过数据层直接打开signal_process.py我在这类项目上踩过的坑是先搞清四份数据分别干什么、标注怎么读、噪声叠加逻辑在哪后面调算法的效率能高一倍。这一章把read_data.py、wfdb_read.py和数据目录的关系讲清楚。2.1 四份数据目录的分工逐搏标注、噪声压力测试、房颤片段与恶性心律工程把原始数据包留在data下并保留了四个来源目录mit-bih-arrhythmia、mit-bih-noise、aha、paf。mit-bih-arrhythmia是公开心电数据库的原始记录双导联动态心电每条记录约 30 分钟带逐搏标注适合做 R 波检测和心律失常识别的基准测试。mit-bih-noise是噪声压力测试集把电极运动伪迹、肌电噪声、基线漂移三类真实干扰叠加到干净心电上专门用来考验滤波器和伪差识别能力。aha目录里是早期室速/室颤竞赛用的短片段波形比较原始适合验证急救场景判断逻辑。paf对应阵发房颤记录提供了房颤发作与不发作的分段区间。目录名内容对算法验证的价值mit-bih-arrhythmia双导联动态心电逐搏标注R 波检测、心律失常分类基准mit-bih-noise三类真实噪声叠加干净心电滤波压力测试、伪差干扰研究aha室速室颤短片段恶性心律失常识别场景验证paf阵发房颤记录房颤发作/不发作的分段验证这四份数据的原始格式并不一致aha早期文件的头格式与 WFDB 不兼容paf的记录时长和标注粒度也不同所以read_data.py的核心作用是把异构来源统一成同一个结构信号数组、采样率、节拍标注数组。不理解这一点直接拿wfdb_read.py去读所有文件会很快碰壁。2.2 wfdb_read.py 的读取套路rdsamp 取信号rdann 取标注wfdb_read.py负责最底层的 WFDB 格式解析。读取一段记录的固定动作是rdsamp拿信号、rdann拿标注。工程数据在本地时优先用本地相对路径避免网络拉取超时把流程打断。import wfdb record_path ./data/100 # 记录名即文件前缀不含 .dat/.hea record wfdb.rdsamp(record_path) # 返回 (信号数组, 头信息字典) ann wfdb.rdann(record_path, atr) # atr 是节拍标注类型 sig record[0] # shape (采样点数, 导联数) fs record[1][fs] # 采样率常见 250 或 360 labels ann.symbol # 每个节拍的标注符号N/V/F/S/A 等 samples ann.sample # 标注对应的全局采样点索引rdsamp返回的二维数组第一列是第一导联第二列是第二导联fs必须显式保存后面滤波器阶数、滑动窗口长度、R 波不应期全部依赖这个值。rdann里的atr表示注释类型有的记录用qrs标注或st事件标注读取前可以先看头文件里的注释列表。samples的单位是采样点而不是毫秒转成秒要除以fs这是新手复现时最容易忽略的地方。2.3 标注质量决定检测基准用噪声库做信噪比压力测试R 波检测算法在干净信号上准确率再高也不算数只有拿到有标注的数据才能算出漏检率和误检率。标注符号里N是正常窦性搏动V是室性早搏S是室上性早搏F是融合波A是房性早搏。如果检测器预期只挑主 QRS把V当成 R 峰通常可以接受但特征提取时要明确标签映射否则异位搏动会污染心率变异性指标。mit-bih-noise的正确用法是把它做成压力测试把干净信号与噪声按不同信噪比叠加重复跑滤波和 R 波检测看检出率随 SNR 下降的曲线。import numpy as np def add_noise(clean_sig, noise_sig, snr_db): # 按目标信噪比叠加噪声snr_db 越低噪声越大 p_signal np.mean(clean_sig ** 2) p_noise p_signal / (10 ** (snr_db / 10)) factor np.sqrt(p_noise / max(np.mean(noise_sig ** 2), 1e-12)) return clean_sig factor * noise_sigSNR 的定义是10 * log10(Ps / Pn)。20 dB 时是轻微噪声5 dB 时信号几乎被噪声淹没。对不同 SNR 重复跑预处理加 R 波检测把准确率随 SNR 的变化画出来就能定量评价滤波的抗噪能力伪差干扰研究也就有了可复现的对照而不是拿偶然的一段抖动数据碰运气。3. 滤波链路signal_process.py 的高通、陷波与低通及参数边界滤波看起来简单实际翻车最多。工程里的signal_process.py把心电预处理做成了可复用函数。这一章按顺序拆解为什么先高通、再陷波、最后低通给出一版可直接替换进工程的函数并说明换采样率后哪些参数必须重算。3.1 滤波顺序比滤波器类型更关键频段划分与 QRS 形态保留心电信号里的噪声有明确频段。基线漂移受呼吸和电极接触影响集中在 0.05 Hz 到 2 Hz工频干扰在 50 Hz部分地区 60 Hz及其谐波肌电噪声和高频毛刺覆盖 10 Hz 以上。QRS 波群的主要能量集中在 5 Hz 到 20 HzT 波在 0.5 Hz 到 2 Hz 附近。所以滤波顺序应该是高通滤掉极低频基线陷波除掉工频低通压制高频在保留 QRS 形态的前提下最大程度去噪。为什么要先做高通基线漂移能量大如果带着漂移去做斜率计算基线上跳会被当成 QRS 候选点误检率直接上升。陷波放在高通之后是因为信号幅值归一化后陷波参数更稳定。低通放最后是因为后续差分或窗口积分会放大高频噪声前置低通能减少衍生干扰。三个滤波器串联时相位响应会叠加工程里统一用零相位滤波避免 QRS 峰值在时间轴上被平移。3.2 可直接替换的预处理函数butter 与 filtfilt 的组合用法从signal_process.py里抽出的核心逻辑整理成下面这个版本import numpy as np from scipy import signal def preprocess_ecg(sig, fs250, highcut40.0, notch_freq50.0): 心电预处理链路高通、陷波、低通全部零相位滤波。 sig: 单导联心电1D 数组 fs : 采样率影响所有归一化截止频率 # 1) 高通 0.5Hz去掉基线漂移和呼吸扰动 b, a signal.butter(2, 0.5 / (fs / 2), btypehigh) sig signal.filtfilt(b, a, sig) # 2) 50Hz 工频陷波Q 值取 30 平衡带宽与形态保持 b, a signal.iirnotch(notch_freq, Q30.0, fsfs) sig signal.filtfilt(b, a, sig) # 3) 低通 40Hz压制肌电和高频噪声 b, a signal.butter(2, highcut / (fs / 2), btypelow) return signal.filtfilt(b, a, sig)butter第二个参数是截止频率归一化到奈奎斯特频率fs / 2所以fs变化时必须重算不能直接拿固定分母的代码换数据跑。阶数取 2 是因为 QRS 频带较宽阶数太高会引入纹波并拉长计算延迟阶数太低过渡带太宽压不住噪声。filtfilt做零相位滤波相当于正向反向各滤一次不产生相位偏移R 波位置不会滑动。iirnotch的 Q 值取 30Q 越大陷波越窄、保留的 QRS 能量越多但工频频率偏移时残留也越多Q 太小会把 QRS 的高频成分削掉一般取 20 到 35 之间。注意如果信号里基线漂移幅度不大高通截止频率可以降到 0.5Hz运动伪迹很重时升到 1Hz代价是压低部分 ST 段的低频斜坡。这是滤波和保留诊断信息之间的折中不是玄学。3.3 换采样率必须重算参数250、360、500 场景下的常见翻车很多复现者从网上抄到一套fs250的函数直接拿到fs360的数据上跑滤波效果大变。根因是归一化截止频率的分母变了0.5 / (fs/2)在 250 与 360 下取值不同。正确做法是把fs作为参数从文件头读出所有与时间相关的窗口都用int(fs * 秒数)计算而不是写死点数。三个具体翻车点第一陷波参数没跟上采样率fs360时仍用 250 的设计参数陷波中心偏移工频残留解决是重算iirnotch。第二移动窗口积分写成固定点数比如窗口写死 15 个点在fs500时窗口只有 30ms同一个 QRS 被拆成两半出现重复检测。第三低通截止频率设到 100Hz肌电残留多误检率上升设到 10HzQRS 被削钝低幅 R 波漏检。40~50Hz 是常见均衡值具体以当前数据的噪声构成为准。4. R波识别与心率计算从能量整形到RR间期附排查记录R 波检测是心率计算和特征提取的前置条件。直接拿原始波形找峰T 波高一点就会翻车。工程里的算法基础是先把 R 波整形为尖峰再做阈值和不应期过滤。4.1 用能量整形把 R 波变成尖峰差分平方加滑动平均QRS 波群里 R 波是斜率最陡的位置。先把信号做一阶差分再取平方放大高频分量最后用 150ms 滑动窗口做平均R 波就变成明显的凸峰T 波被压下去这是常见 QRS 检测论文里的基础结构。import numpy as np def qrs_energy(sig, fs): 将心电转成 QRS 能量曲线供后续找峰使用 diff np.diff(sig, prepend0) sq diff ** 2 win int(0.15 * fs) kernel np.ones(win) / win return np.convolve(sq, kernel, modesame) def detect_r_peaks(sig, fs, thr_factor0.35, refractory0.25): energy qrs_energy(sig, fs) thr thr_factor * np.mean(energy) r_peaks [] last_pos -9999 for i in range(len(energy)): if energy[i] thr: continue # 不应期过滤避免同一 QRS 被重复检出 if i - last_pos int(refractory * fs): last_pos i continue # 在候选点附近找幅值最高点作为 R 峰 lo max(0, i - int(0.1 * fs)) hi min(len(sig), i int(0.1 * fs)) pos lo int(np.argmax(sig[lo:hi])) r_peaks.append(pos) last_pos i return np.array(r_peaks)diff计算斜率平方放大高频滑动窗口把同一 QRS 的多重峰合二为一。thr_factor默认 0.35意思是阈值取平均能量的 35%噪声多时提高尖峰清晰时降低。refractory是 0.25 秒的绝对不应期对应 240bpm 的极限心率超过这个频率的候选峰大概率是噪声或 T 波。4.2 从 RR 间期算心率与心率变异性中位数比平均值稳拿到 R 峰索引后计算 RR 间隔很直接但单次误检会让平均值心率跳变。工程实践里通常用中位数抗干扰同时输出常用心率变异性指标。rr np.diff(r_peaks) / fs if len(rr) 2: return None hr_instant 60.0 / rr # 每搏瞬时心率bpm hr_median float(np.median(hr_instant)) # 抗单次误检的心率估计 sdnn float(np.std(rr)) # RR 间期标准差 rmssd float(np.sqrt(np.mean(np.diff(rr) ** 2))) # 相邻 RR 差均方根 pnn50 float(np.mean(np.abs(np.diff(rr)) 0.05) * 100) # 相邻差超50ms占比rr的单位是秒瞬时心率单位是 bpm。sdnn反映长程变异性rmssd反映短程变化pnn50是相邻 RR 差超过 50ms 的比例。房颤时rmssd与pnn50通常显著升高这两个特征可以直接喂给第 5 章的分类器。4.3 排查误检、漏检、心率跳变的四类常见翻车现象一同一 QRS 被检出两个 R 峰心率变成实际的两倍。原因是能量窗口太短QRS 的上升沿和下降沿被拆成两个峰或者不应期小于 0.2 秒。解决是把qrs_energy的窗口设在 120~200ms不应期保持不小于 0.25 秒。如果用scipy.signal.find_peaks给distanceint(0.25*fs)即可。现象二高 T 波被当成 R 峰心率忽高忽低。原因是阈值基于能量均值T 波形态较高时能量与 R 波接近。解决是提高幅值阈值或者改用 R 波上升沿的陡峭特征加形态宽度判断候选峰宽度与 QRS 接近时才保留。现象三低幅 R 波漏检房颤和室速片段里尤其明显。原因是整段信号的均值阈值偏高低幅 R 波被淹没在噪声里。解决是改成滑动窗口自适应阈值用最近 1s 能量的滑动中位数替代整段均值。现象四运动伪差或电极摆动被当成 R 波心率瞬间飙上去。原因是伪差能量与 QRS 频谱重叠单导联无法分辨。解决是双导联投票两导联各自检测后只保留位置一致的峰或者联合局部 SNR 判断低于阈值的片段直接丢弃。这正是mit-bih-noise存在的意义跑一遍不同 SNR 的压力测试这类问题能提前暴露。5. 特征提取与心律失常分类从 RR 间期到房颤、室颤室速与伪差识别特征提取和分类是工程的价值终点但上手顺序要稳妥先搭特征矩阵再选传统方法还是深度学习最后必须搞清楚房颤、室颤、室速和伪差的边界否则分类结果只是数字游戏。5.1 特征矩阵怎么搭时域、形态与非线性三层特征单靠 RR 规律性不足以区分类别因为运动伪差也会让 RR 变乱。特征矩阵分三层。特征类别代表特征计算来源主要判别方向时域中位 RR、SDNN、RMSSD、pNN50RR 间期序列房颤、心动过速/过缓形态QRS 宽度、QRS 幅度、T/R 比值滤波后信号切片室速、室早非线性样本熵、RR 间期熵RR 序列房颤、室颤样本熵计算时m取 2r取0.15 * std(RR)是工程里常见的配置。形态特征要从每个 R 峰切窗提取窗口左边界取峰前 50ms右边界取峰后 100ms。def extract_features(sig, fs, r_peaks): rr np.diff(r_peaks) / fs if len(rr) 2: return None rmssd float(np.sqrt(np.mean(np.diff(rr) ** 2))) pnn50 float(np.mean(np.abs(np.diff(rr)) 0.05) * 100) qrs_widths [] for p in r_peaks: lo max(0, p - int(0.05 * fs)) hi min(len(sig), p int(0.1 * fs)) seg sig[lo:hi] if len(seg) 0: continue # 用局部幅值超过0.5倍峰值的点数估算宽度 qrs_widths.append(np.sum(np.abs(seg) 0.5 * np.max(np.abs(seg))) / fs) return { mean_rr: float(np.mean(rr)), rmssd: rmssd, pnn50: pnn50, mean_qrs_width: float(np.mean(qrs_widths)) if qrs_widths else 0.0, }特征提取完成后每个 30 秒窗口可以用一条特征向量表示。宽度估算只用于相对比较不追求精确到毫秒因为分类器更依赖特征之间的相对差异。5.2 分类链路basic_alg 的传统方法与 deep_learning 目录的取舍工程同时提供了basic_alg和deep_learning两套路径。传统方法先用固定时间窗拼特征向量再喂给 SVM 做交叉验证结果可解释且对样本量要求低。from sklearn.svm import SVC from sklearn.model_selection import cross_val_score X feature_matrix # shape (n_segments, n_features) y segment_labels # 0正常, 1房颤, 2室颤, 3室速 clf SVC(kernelrbf, C10.0, gammascale) scores cross_val_score(clf, X, y, cv5, scoringf1_macro) print(fmean F1: {scores.mean():.3f})RBF 核处理非线性边界C10控制误分类惩罚gammascale自动按特征方差缩放。交叉验证用 5 折即可但要注意按记录分组而不是按段分组否则同一条记录里相邻的窗口会泄漏到训练集。deep_learning目录里的网络更适合端到端学习输入原始波形或小波时频图不再依赖人工特征小样本下很容易过拟合所以工程实践里通常先跑basic_alg建立基准再看深度模型有没有实际增益。5.3 房颤、室颤室速与伪差的判定边界先排除伪差再下结论房颤的典型特征是完全无规律的 RR 间期、P 波消失并被细小 F 波替代可以用 30 秒窗口的rmssd和样本熵描述。室颤是心室的无序电活动没有完整 QRS 形态体表记录为波浪状粗 f 波频率约 150~250 次/分。室速表现为连续宽大畸形的 QRS 波群心率常在 100~220 次/分。伪差里的电极脱落是直线或瞬间大幅跳变最容易识别运动伪差的波形有时看起来非常像室颤单靠幅值阈值一定会误报。def heuristic_triage(feat): 研究演示用启发式规则不能作为临床诊断依据 if feat.get(rmssd, 0) 0.12 and feat.get(entropy, 0) 1.5: return suspicious AF / highly irregular if feat.get(mean_hr, 0) 150 and feat.get(qrs_width, 0) 0.12: return suspicious VT if feat.get(flatline_ratio, 0) 0.5: return artifact / lead off return unknown这些阈值都是演示性质实际使用必须结合带标注的数据库调参。工程里把伪差干扰研究样本纳入测试集目的是让模型学会输出“不确定”而不是强行分到某一类。在临床落地场景伪差误报成室颤是非常危险的错误方向分类链路里保留 artifact 分支比追求单类准确率更重要。6. 可视化走读与复现验证用叠加绘图把每一步结果钉在图上6.1 ecg_display.py 的画法一张图同时看到原始、滤波和 R 波位置ecg_display.py提供的核心能力是把预处理前后波形放在同一张图里对比并把检测到的 R 峰位置直接画在信号上。调试时能看到滤波有没有削掉幅值、R 峰对齐是否准确而不是对着数组猜。import matplotlib.pyplot as plt import numpy as np def plot_ecg_pipeline(sig_raw, sig_filtered, fs, r_peaksNone, titleECG Pipeline): t np.arange(len(sig_filtered)) / fs fig, ax plt.subplots(figsize(12, 5)) # 原始信号浅灰滤波后信号黑色直观对比形态保留程度 ax.plot(t, sig_raw, colorlightgray, lw0.8, labelraw) ax.plot(t, sig_filtered, colorblack, lw1.2, labelfiltered) if r_peaks is not None: ax.scatter(t[r_peaks], sig_filtered[r_peaks], s50, colorred, zorder3, labelR peaks) ax.set_xlabel(time (s)) ax.set_ylabel(mV) ax.legend(locupper right) ax.set_title(title) fig.tight_layout() return fig使用要点原始信号用浅灰色细线滤波后信号用黑色粗线R 峰用红色散点。三个图层叠加后能确认三件事滤波后 QRS 形态没有明显变形R 峰与真实 QRS 对齐T 波高耸的段有没有被误标。这个函数可以直接替换工程里ecg_display.py的同类绘图入口。6.2 复现验证先对拍标注再谈算法指标工程包里的tests目录已经提供了基本测试工程但自己动手验证一套流程更稳妥。固定做法是读取标注 → 跑预处理 → 跑 R 波检测 → 与标注样本对拍统计 TP、FP、FN再算灵敏度和阳性预测值。def evaluate_detection(r_peaks, ann_samples, fs, tol_ms150): tol int(tol_ms / 1000 * fs) ann_set set(ann_samples) tp sum(1 for p in r_peaks if any(abs(p - a) tol for a in ann_samples)) fp len(r_peaks) - tp fn sum(1 for a in ann_samples if all(abs(a - p) tol for p in r_peaks)) se tp / max(tp fn, 1) ppv tp / max(tp fp, 1) return {sensitivity: se, positive_predictive_value: ppv}容差窗口取 ±150ms这是心电算法验证里常见的对拍口径。从那以后我每次拿到新的数据源上来第一件事就是强制跑一遍“读数据 → 滤波 → R 峰检测 → 对标注”只看两个数字灵敏度和阳性预测值。这个习惯帮我挡掉了至少三次格式不兼容、标注错位和采样率写死带来的假实验结果。希望帮到你。本文还有配套的精品资源点击获取