简介:面向心电信号处理与算法研发人员,这份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 到 30Hz,T 波集中在 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, lowcut=0.5, highcut=45.0, order=4): # 计算 Nyquist 频率并做归一化 nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq # 生成巴特沃斯带通系数 b, a = signal.butter(order, [low, high], btype='band') # filtfilt 做零相位滤波,保持 R 波位置不变 return signal.filtfilt(b, a, data, axis=0)参数含义要按实际数据来定。lowcut=0.5 是基线漂移的分界线,低于这个频率的缓慢变化会被滤掉;highcut=45.0 是高频截止点,大于它的部分会被衰减。order=4 表示四阶滤波器,阶数越高过渡带越窄,但同时数值稳定性会变差。对 500Hz 采样率的心电数据,四阶巴特沃斯带通通常是够用的。
调用时如果发现输出两端有明显震荡,先检查数据里有没有 nan 或 inf,filtfilt 对异常值非常敏感。还可以把数据先按片段切好,保证片段长度不小于滤波器阶数的 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_ms=120, refractory_ms=200, threshold_ratio=0.6, invert=False): x = -ecg if invert else ecg # 一阶差分突出 R 波陡峭斜率 diff_sig = np.diff(x, prepend=x[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, mode='same') # 自适应阈值:取包络峰值的比例 height = threshold_ratio * np.max(env) # 不应期:两个 R 峰之间的最小样本数 distance = int(fs * refractory_ms / 1000) peaks, _ = signal.find_peaks(env, height=height, distance=distance) return peaks, env这段代码里,prepend=x[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 BPM | 1.5s | 间期成倍增长,心率显示减半 |
| 60 BPM | 1.0s | 基线平稳时正常 |
| 120 BPM | 0.5s | 间期减半,心率显示翻倍 |
| 200 BPM | 0.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, fs=500, tol_ms=20): 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 波位置,确认没有因为调阈值带来新的误判,这一轮参数调整才算真正收尾。
本文还有配套的精品资源,点击获取