基于Python的心电信号处理：滤波、R波检测与心率计算

简介面向心电信号处理与算法研发人员这份Python工程围绕心电图分析全流程展开覆盖信号滤波、R波检测、心率计算、特征提取、心律失常与房颤/室颤室速识别等核心环节并提供伪差干扰处理思路与可视化模块适合生物医学工程、算法测试与临床数据挖掘场景。压缩包共91个文件以17个py源码脚本为主体辅以心电数据文件dat/hea/atr、7个xws文件、4个pdf文献、3个md说明及jpg/png可视化结果图整体约30.81MB目录按数据、算法、文档和测试组织便于直接运行与二次开发。目前已有672人学习下载资源内包含完整的测试工程和基础算法demo并内置MIT-BIH数据读取与转换脚本读者可对照标准数据库快速验证R波定位与心率计算效果并通过特征分布图理解分类依据是一份兼顾教学与工程落地的实用代码库。1. 使用 Python 编写心电算法的前提先分清信号、噪声和采样率原始心电记录在接入 Python 之前先要回答一个问题它要处理的是哪种噪声很多刚把采集电路跑通的人习惯把整段信号丢进单层低通滤波就找 R 波最终检测结果往往不稳。心电里真正干扰 R 波判断的是 50Hz 工频、肌电高频分量和呼吸造成的基线漂移这三类噪声混在一起时幅值一点不比 QRS 小。换句话说滤波、R 波识别和心率计算是一条完整流水线先用带通滤波将 QRS 所在的频段圈出来再用滑动窗口在时域确认峰值位置最后依靠 RR 间期的统计关系排除漏检和误检。下面从 numpy 数组和 scipy.signal 入手把三个环节串成可直接套用的代码适合刚把心电采集电路跑通、准备进入算法阶段的工程师。2. 心电滤波带通、滑动窗口和卡尔曼滤波怎样配合说到心电滤波先要建立一个认知不存在一个滤波器能同时解决基线漂移、工频和高频肌电这三类问题合理的做法是按频段拆开处理让每一级滤波只做一件事。2.1 心电噪声的频带分布与低通、高通边界心电各成分在频谱上有比较明确的分界。QRS 波群的主要能量集中在 5Hz 到 30HzT 波集中在 2Hz 以下P 波在 5Hz 左右。基线漂移由电极位移和呼吸引起频率基本在 0.5Hz 以下50Hz 工频来源于电网虽然硬件端会做屏蔽和右腿驱动数字端仍可能残留肌电噪声分布在几十到几百 Hz属于宽带随机干扰。正因为这些频带可分滤波的第一选择是带通。比较常用的通带是 0.5Hz 到 45Hz低端 0.5Hz 切掉基线漂移高端 45Hz 保留 QRS 形态同时压掉大部分肌电和高频噪声。有些心电设备只用 1Hz 高通那是针对监护场景希望尽量减少波形失真如果后续要做 P 波分析高端最好放到 30Hz 而不是 45Hz否则高频肌电会混进 P 波频段常见做法是保 0.5Hz 到 30Hz。这一段的关键是明确边界而不是追求滤波器阶数。心电信号采样率是 360Hz 时 Nyquist 频率为 180Hz带通到 45Hz 意味着归一化截止频点只有 0.25这个距离给滤波器阶数留了足够余量。采样率只有 125Hz 时则要注意45Hz 高端截止点会贴近 Nyquist此时应把高端改成 30Hz 左右避免滤波器边界效应扭曲波形。2.2 用 scipy.signal 实现巴特沃斯带通滤波带通滤波的 Python 实现并不复杂核心是 butter 生成系数、filtfilt 做零相位滤波这两个函数。filtfilt 会沿时间轴正向和反向各过滤一次抵消相位延迟滤完的波形与原始信号在时间位置上不产生偏移这对后面 R 波识别非常重要因为 R 峰坐标要对应到原始波形上的真实位置。下面的函数是一段可以直接抄进项目里的最小实现import numpy as np from scipy import signal def ecg_bandpass(data, fs, lowcut0.5, highcut45.0, order4): # 计算 Nyquist 频率并做归一化 nyq 0.5 * fs low lowcut / nyq high highcut / nyq # 生成巴特沃斯带通系数 b, a signal.butter(order, [low, high], btypeband) # filtfilt 做零相位滤波保持 R 波位置不变 return signal.filtfilt(b, a, data, axis0)参数含义要按实际数据来定。lowcut0.5 是基线漂移的分界线低于这个频率的缓慢变化会被滤掉highcut45.0 是高频截止点大于它的部分会被衰减。order4 表示四阶滤波器阶数越高过渡带越窄但同时数值稳定性会变差。对 500Hz 采样率的心电数据四阶巴特沃斯带通通常是够用的。调用时如果发现输出两端有明显震荡先检查数据里有没有 nan 或 inffiltfilt 对异常值非常敏感。还可以把数据先按片段切好保证片段长度不小于滤波器阶数的 3 倍否则边界填充区会占很大比例波形开头几十个点基本不可用。这类问题在做实时心电监测时最容易出现因为实时场景往往只给出一两秒的短窗口滤波边界的表现比长序列离线和差得多。2.3 滑动窗口滤波与卡尔曼滤波在实时场景里的取舍带通之后信号已经比较干净但很多人还想再加一道时域平滑这时滑动窗口滤波就被拿出来用了。滑动窗口滤波的本质是求窗口内平均值等价于一个低通滤波。窗口大小的选择直接影响 R 波形态窗口越大平滑效果越强但 R 波顶峰会被压扁时间坐标也会被拉偏。对心电来说窗口长度按 QRS 宽度来定比较合理取采样率对应的 100ms 到 150ms也就是 500Hz 下取 50 到 75 个点而不是随手取一个奇数。卡尔曼滤波则是另一个方向它把心电信号当作一个带有观测噪声的动态系统用状态方程递推估计真实波形。状态向量可以设为信号幅值和幅值变化率通过过程噪声协方差 Q 和测量噪声协方差 R 控制滤波平滑度。卡尔曼的优势在于逐点递归内存占用固定非常适合嵌入式设备上的实时处理但它的两个协方差矩阵需要经验调整调不好时容易把 R 波也当作噪声滤掉。实际工程项目里离线分析先跑带通在线实时处理则带通加滑动窗口更省算力卡尔曼滤波更适合嵌入式心电或者其他生理信号采集场景。三种方法的直观差异可以参考下表方法频响特性相位延迟对 R 波影响适用场景巴特沃斯带通带通过渡带陡零相位filtfilt基本不改变形态离线分析滑动窗口滤波低通窗口越大越平滑有滞后R 波变钝幅度降低实时流式处理卡尔曼滤波由 Q/R 决定有递归滞后可能压低 R 波峰值嵌入式在线估计我一般是这样分工的离线分析只做带通实时场景带通降阶数、加滑动窗口只有需要逐拍在线输出的嵌入式工程才引入卡尔曼。顺序上永远先带通再用滑动窗口卡尔曼放在最后考虑不要一上来就把三种滤波全叠上去叠多了波形会变得很“糊”R 波定位反而更难。3. R 波识别差分、平方、自适应阈值与不应期滤波完成后进入核心环节找出每个 QRS 波群里的 R 峰位置。直接对滤波后的波形找局部最大值很容易把 T 波当 R 波因为某些导联下 T 波幅值可以接近甚至超过 R 波。所以 R 波识别要先做波形变换把 QRS 特有的陡峭斜率放大再用阈值和不应期约束。3.1 QRS 波形变换一阶差分与平方包络QRS 波群在形态上最大的特点是斜率大T 波再高也是缓慢变化的。利用这一点对滤波后的信号做一阶差分R 波所在位置的差分值会明显高于 T 波差分同样会把噪声放大所以紧接着要做平方把差分结果变成正值再通过一个短滑动窗口积分成包络让 R 波位置对应一个平滑的凸起峰。这就是 Pan-Tompkins 方法的核心思路。具体步骤通常是这样一阶差分、逐点平方、窗口积分。窗口长度一般取 150ms 左右太短则包络上残留锯齿太长则两个相邻 R 波可能合并成一个峰。做平方之前先检查信号极性如果整段信号的 R 波朝下需要先乘 -1 再继续否则平方后位置会偏到 S 波上去。3.2 自适应阈值与不应期的设计阈值是 R 波识别里最容易翻车的参数。固定阈值在高噪声段会漏检在低噪声段会把 T 波误判为 R 波所以工程上普遍用自适应阈值。一种简单有效的做法是取当前窗口内包络最大值的 60% 作为阈值然后边移动窗口边更新阈值基准。心率突然加快时上一拍的最大值可能偏高导致下一拍被压住所以需要限制阈值的更新速度。不应期同样关键。生理上 R 波之后约 200ms 内不会出现另一个 QRS因此检测到一次 R 峰后至少在 200ms 内不再接受新的峰值。这段不应期由采样率换算成样本数例如 500Hz 采样率对应 100 个采样点。不应期的存在能有效消除同一个 R 波在包络上出现双峰造成的重复计数。3.3 R 波检测的 Python 实现与参数说明结合上面的思路可以写一个完整的检测函数。这里用 scipy.signal 的 find_peaks 替代手写阈值循环逻辑更清晰也方便调参from scipy import signal def detect_r_peaks(ecg, fs, smooth_ms120, refractory_ms200, threshold_ratio0.6, invertFalse): x -ecg if invert else ecg # 一阶差分突出 R 波陡峭斜率 diff_sig np.diff(x, prependx[0]) # 平方放大差分结果让能量集中 squared diff_sig ** 2 # 滑动窗口积分形成单峰包络 win_len max(1, int(fs * smooth_ms / 1000)) kernel np.ones(win_len) / win_len env np.convolve(squared, kernel, modesame) # 自适应阈值取包络峰值的比例 height threshold_ratio * np.max(env) # 不应期两个 R 峰之间的最小样本数 distance int(fs * refractory_ms / 1000) peaks, _ signal.find_peaks(env, heightheight, distancedistance) return peaks, env这段代码里prependx[0] 是为了让差分结果与原始信号长度一致避免错位。平方之后数据范围会变大所以 threshold_ratio 取 0.6 时阈值是包络最大值的 60%这个值需要根据噪声水平调整。find_peaks 的 distance 参数等于不应期的样本数它保证检测出的两个峰值之间至少间隔 200ms。调用这个函数时建议把原始波形和处理后的包络一起画出来对比。如果检测结果偏多多半是 threshold_ratio 设低了偏少则说明阈值太高或不应期太长。对一段 10 秒、采样率 500Hz 的数据检测到的 peaks 数组里每个元素是一个样本下标除以采样率就得到对应秒数。这里列的步骤其实就是平时调试时的标准动作先看包络形状再调阈值最后看不应期是否挡住重复峰。4. 心率计算RR 间期、BPM 与异常值剔除得到 R 峰位置之后心率计算看起来只是做个除法实际工程里麻烦的是漏检和误检带来的偏差。心率的正确单位是次/分英文简写 BPM它是通过 RR 间期换算出来的。4.1 RR 间期与平均 BPM 的关系相邻两个 R 峰之间的时间间隔叫 RR 间期单位是秒。heart 率的公式是 BPM 60 / RR。比如采样率 500Hz 时相邻 R 峰相隔 500 个采样点对应 1 秒那心率就是 60BPM。工程上更稳妥的做法是先算出所有 RR 间期再做平均而不是直接把单位时间内的峰数当作心率后者在心率变化较快时不准确。计算时需要特别注意采样率误差。很多便携采集设备的标称采样率与实际晶振频率之间有百分之零点几的偏差短时间看不出来但统计整段 RR 间期时误差会被放大。如果发现连续多段数据的心率都系统性偏离真实值先用一个已知频率的信号发生器标定采样率再回填到算法里。4.2 漏检和误检对心率统计的影响漏检时两个真实 R 峰之间只检测到一个峰RR 间期变成原来的两倍换算出的心率直接减半。误检则相反一个 R 波被记成两次RR 间期减半心率翻倍。所以心率计算前必须对 RR 间期序列做异常剔除这一步的价值不亚于 R 波识别本身。常用做法是取 RR 间期的中位数作为基准因为中位数比均值更抗异常值然后剔除偏离中位数超过 30% 的间期。这个比例覆盖了正常心率波动范围又能拦住典型的减半、翻倍错误。剔除后如果剩余间期少于原来的一半说明检测质量太差算法应该输出一个“信号质量不足”的标记而不是给一个看似正常的心率数字。心率区间对应 RR 间期常见误检表现40 BPM1.5s间期成倍增长心率显示减半60 BPM1.0s基线平稳时正常120 BPM0.5s间期减半心率显示翻倍200 BPM0.3s接近不应期检测不稳定4.3 心率计算的 Python 实现下面给出从 R 峰下标到心率输出的完整代码注意区分单个 RR 间期的心率和整体平均心率def compute_bpm(r_peaks, fs): # r_peaks 是 R 峰样本下标数组 if len(r_peaks) 2: return [], float(nan) # 换算成秒 rr_interval np.diff(r_peaks) / fs # 逐拍心率 bpm_series 60.0 / rr_interval # 用中位数剔除异常 RR 间期 med_rr np.median(rr_interval) keep np.abs(rr_interval - med_rr) / med_rr 0.3 rr_clean rr_interval[keep] if len(rr_clean) 0: return bpm_series, float(nan) # 平均心率由干净 RR 间期计算 avg_bpm 60.0 / np.mean(rr_clean) return bpm_series, avg_bpmbpm_series 是逐拍心率序列可以用于心率变异性分析avg_bpm 是整段平均心率对应监测界面上那个大数字。keep 条件里的 0.3 是剔除比例想要更严格的剔除就调到 0.2噪声高时可以适当放到 0.4。处理完的 bpm_series 如果长度明显小于 rr_interval说明原始检测存在较多误检应该回溯检查 R 波识别参数。长时监测场景下平均心率可以用 10 秒窗口滑动计算。每次只取最近 10 秒内的 R 峰位置算出该窗口的 RR 间期序列再做异常剔除和平均这种滑动窗口方式能更好反映心率随时间的变化趋势。5. 用验证脚本校准滤波与 R 波识别参数参数能不能用得靠一组带标注的 R 波位置来验证。所谓标注就是人工在波形上标出每个 R 峰对应的样本下标或者使用公开心电数据集中已经给出的参考结果。验证的核心是比较算法结果和标注位置是否在允许误差范围内。5.1 灵敏度与阳性预测值的计算评估一个 R 波检测器常用两个指标灵敏度和阳性预测值。灵敏度等于正确检出的 R 波数除以标注总数阳性预测值等于正确检出的 R 波数除以算法检出的总数。允许的误差一般取 20ms 到 150ms误差窗口越小对算法要求越严格。下面的评估函数可以直接用来对比def evaluate_detection(annotation, detection, fs500, tol_ms20): tol int(tol_ms / 1000 * fs) tp 0 for ann in annotation: if np.min(np.abs(detection - ann)) tol: tp 1 fp len(detection) - tp fn len(annotation) - tp sensitivity tp / len(annotation) if annotation else 0 ppv tp / len(detection) if detection else 0 return sensitivity, ppv, fp, fn这个函数把每个标注位置和最近的检测位置做差落在容差窗口内就计为正确。使用 20ms 容差来衡量 QRS 定位精度比较合理但初次调试可以用 150ms先确认找对位置再逐步收紧。如果灵敏度低说明漏检多优先调低 threshold_ratio如果阳性预测值低说明误检多优先检查不应期和滤波高端截止频率。5.2 调整参数时的三个优先顺序调整参数时我通常按三个方向依次排查。第一带通滤波的高端是否过高50Hz 工频有没有残留残留的工频会让包络上出现周期性小峰find_peaks 很容易把它们误认为 R 波。第二smooth_ms 是否匹配 QRS 宽度窗口太短包络上有锯齿窗口太长相邻 R 波合并。第三threshold_ratio 与实际信号幅值的关系心电信号幅值不稳定时用中位数代替最大值做分母会更稳定。还有一个容易被忽略的问题是滤波窗口与数据长度的关系。filtfilt 在做离线分析时表现很好但放到实时心电监测中它只能处理当前已采集的数据段每来一个新采样点就要重算一次这会引入明显延迟。实时场景更推荐用一阶低通加滑动窗口的组合牺牲一部分滤波效果换低延迟。参数调试不是一次性完成的。换一份数据集、换一个采样率、换一个导联位置原来能用的参数就要重新评估。把上面这个评估函数保存成一个独立的脚本每次改参数后跑一遍对比灵敏度和阳性预测值的变化比肉眼扫波形高效得多。等到灵敏度和阳性预测值都超过 95%再回到原始波形上看几处典型的 T 波位置确认没有因为调阈值带来新的误判这一轮参数调整才算真正收尾。本文还有配套的精品资源点击获取