Python心电算法实战:从滤波到R波检测与心率计算

发布时间:2026/10/18 12:30:51
Python心电算法实战:从滤波到R波检测与心率计算 简介面向心电信号处理与算法研究的Python项目覆盖滤波、R波识别、心率计算、特征提取、心律失常分类、房颤与室颤室速识别等核心环节并内置了可视化模块与独立测试工程。项目借助SciPy、PyWavelets、scikit-learn等库完成数字滤波、小波变换R波定位以及机器学习模型训练适合生物医学工程学习者、算法工程师和需要快速验证ECG算法的研究者。资源包共91个文件主要包括Python脚本、MIT-BIH格式心电数据文件.dat/.hea/.atr、算法说明文档与相关论文PDF压缩包大小30.81MB目录结构清晰便于按模块查找。脚本涵盖信号预处理、R波检测、特征选择与分类预测数据文件可配合脚本在标准数据库上复现实验PDF与markdown文档则提供理论推导和实现细节。目前已有673人学习下载既可当作教学案例逐行研读也能直接作为心电分析原型的开发基础。1. 心电算法看着像玄学本质是一条可复现的信号处理链路第一次用 Python 接触心电算法的人十有八九会先打开一个深度学习教程把原始波形直接丢进模型。但真实工程里滤波、R波识别、心率计算这三件事90% 是靠成熟的信号处理方法完成的。把一条 ECG 原始记录变成一份能看心率数值的处理流程通常只需要 NumPy/SciPy几百行代码就能跑通。这个方向适合三类人一是在做可穿戴设备数据处理的开发者二是刚进入医疗信号领域的算法工程师三是想把手动看波形的分析过程自动化的研究人员。心电算法并不神秘前提是你先把信号读懂。2. 读懂心电信号采样率、导联和滤波的边界2.1 心电信号里有哪些成分QRS、基线漂移和工频干扰心电信号不是一条干净的曲线它由不同频率成分叠加而成。P 波和 T 波属于慢成分能量主要集中在 0.5~5 Hz 附近QRS 波群是心电里最陡峭也最明显的部分主要能量在 5~20 Hz。基线漂移是身体呼吸、电极移动造成的低频起伏通常低于 0.5 Hz工频干扰来自供电线路在中国大陆一般是 50 Hz部分设备环境里还有 60 Hz 及其谐波肌肉抖动带来的肌电干扰频带更宽从 20 Hz 一直延伸到 100 Hz 以上。R波识别之所以很少直接找一个全局最大值就是因为信号里混杂着这些干扰。基线漂移会让整段波形上下缓慢变化工频干扰会在 QRS 上叠加一层细小的锯齿肌电干扰会让波峰毛刺变大。滤波的目标不是把所有非 QRS 成分都清掉而是在压制干扰的同时保留 P 波、T 波和 QRS 的形态。如果滤波过头R 波会被削平后面做 R波识别时反而更难。2.2 用 Python 读取原始心电数据先确认采样率写任何心电算法之前第一件事是确认采样率。采样率不是抄过来的而是从数据里核对出来的。很多心电数据以 CSV 格式保存一列是时间一列是幅值。用 NumPy 读入后先看时间轴间隔不要直接写死 250 Hz。import numpy as np # raw.csv 里第一列是时间(s)第二列是心电幅值(mV) data np.loadtxt(raw.csv, delimiter,, skiprows1) time data[:, 0] ecg data[:, 1] # 用时间间隔的中位数估计采样率避免个别时间戳异常影响结果 dt np.median(np.diff(time)) fs int(round(1.0 / dt)) print(采样率:, fs, Hz) print(信号长度:, len(ecg), 个采样点)这段代码里np.diff(time)算的是相邻两个时间点之间的间隔np.median比np.mean更抗个别跳变。如果时间轴是等间距的fs就是数据真实的采样率。如果数据源本身给出了采样率也要用上面的方式校验一遍因为很多采集设备实际输出频率和配置并不完全一致。采样率一旦算错后边所有滤波器的截止频率都会偏移。2.3 滤波前必须先确认的三个参数滤波不是随便套一个低通或高通就完事。每次拿到新数据我都会先确认三个参数采样率、滤波频带、滤波器阶数。这三个参数决定了你能保留什么、滤掉什么。参数常见做法判断依据采样率 fs250 Hz 或 500 Hz以数据时间轴实际间隔为准高通截止频率0.5~1 Hz低于此频率主要是基线漂移陷波频率50 Hz 或 60 Hz看供电环境或频谱尖峰滤波器阶数高通 2~4 阶陷波 Q 值 20~30阶数太高容易产生振铃和相位畸变心电滤波有一个反直觉的原则滤波不是越干净越好。能看清 QRS 波、保留 P 波和 T 波才是合适的滤波结果。如果为了压基线漂移把高通截止频率调到 5 HzR 波虽然还在P 波和 T 波会被明显削弱后续如果要做心电特征分析数据就不够用了。注意滤波器的截止频率和采样率有关系。同样是 0.5 Hz 高通在 250 Hz 和 500 Hz 采样率下传给滤波器的归一化截止频率完全不同参数不能直接搬。3. 用 Python 把滤波写稳基线漂移、工频干扰和参数选择3.1 去掉基线漂移高通滤波与中值滤波两种写法基线漂移是最常见的心电干扰。呼吸造成的胸部起伏会让整条波形以 0.3~0.5 Hz 的频率上下摆动如果不去除R 波检测阈值会被抬到一个无效的位置。处理基线漂移有两种常用写法一种是高通滤波一种是中值滤波。高通滤波适合基线漂移频率相对稳定的场景计算快参数直观。我一般用二阶 Butterworth 高通截止频率 0.5 Hz再用sosfiltfilt做零相位滤波。零相位的意思是滤波后波形没有整体延迟R 波位置不会发生移动这对后面 R波识别非常关键。from scipy.signal import butter, sosfiltfilt def highpass_filter(ecg, fs, cutoff0.5, order2): # 归一化截止频率 截止频率 / (采样率/2) sos butter(order, cutoff / (fs / 2), btypehighpass, outputsos) return sosfiltfilt(sos, ecg) ecg_hp highpass_filter(ecg, fs, cutoff0.5, order2)这里order2是 Butterworth 滤波器的阶数sosfiltfilt会正反各滤一次等效阶数更高但相位延迟为零。阶数不建议超过 4因为阶数越高滤波器在截止频率附近的振铃越明显R 波前后会出现人为的起伏。中值滤波的思路不一样它不是把低频成分滤掉而是直接从原始信号里把一条基线估计出来再用原始信号减去这条基线。中值滤波对形态不规则、偶发突变的基线漂移更稳因为高通滤波在信号边缘容易产生过冲。from scipy.ndimage import median_filter def remove_baseline_median(ecg, fs, window_seconds0.2): # 窗口长度取一个 QRS 宽度以上才能保住 QRS window int(fs * window_seconds) if window % 2 0: window 1 baseline median_filter(ecg, sizewindow, modenearest) return ecg - baseline中值窗口的window_seconds一般取 0.2 秒左右。如果窗口设得太短比如 0.05 秒中值会把 QRS 波本身也当成高频细节剔除R 波幅度会明显下降。如果设得太长估计出来的基线跟不上漂移速度去不干净。modenearest是为了处理边界避免波形两端出现空值。3.2 抑制工频干扰陷波器与加权移动平均工频干扰表现为波形上叠加的细密高频纹波频率在 50 Hz 附近。判断它是否存在最直接的办法是看原始波形上有没有一条稳定的细锯齿或者对信号做频谱分析看 50 Hz 处是否有明显尖峰。确认之后用陷波器是最常见的做法。from scipy.signal import iirnotch, filtfilt def notch_ecg(ecg, fs, freq50.0, q30.0): b, a iirnotch(freq, q, fs) return filtfilt(b, a, ecg) ecg_clean notch_ecg(ecg_hp, fs, freq50.0, q30.0)q是品质因子决定陷波带宽。Q 值越大陷波频带越窄对相邻频率的影响越小Q 值太小会把 50 Hz 附近的心电成分也削掉。一般取 20~30 比较合适。如果现场供电环境是 60 Hz把freq改成 60.0 就行其他参数不用动。还有一种老式做法是用移动平均做陷波窗口长度取一个工频周期内的采样点数。比如 50 Hz、250 Hz 采样率下一个周期是 5 个点移动平均窗口可以取 5。这个方法实现简单但会同时模糊 QRS 波形而且对谐波压不干净。我一般只在采样率很低、陷波器数值不稳定的情况下才用移动平均能上陷波器就优先陷波器。3.3 滤波效果怎么验证滤波做完不能只看一眼感觉“顺眼”就够了。我的验证习惯是时域和频域各看一次。时域看滤波后的 R 波是否还尖锐、P 波和 T 波是否还在频域看基线漂移频段和工频处的能量是否真的降下来了。from scipy.signal import welch f, psd welch(ecg_clean, fsfs, npersegmin(4 * fs, len(ecg_clean))) # 检查 0.5 Hz 以下的能量是否被压下去 lowband np.sum(psd[(f 0.01) (f 0.5)]) # 检查 50 Hz 附近的能量是否明显低于原始值 power_band np.sum(psd[(f 48) (f 52)]) print(低频能量:, lowband, 工频带能量:, power_band)welch返回频谱的横轴f和功率谱密度psdnperseg决定每个分段的长度取 4 秒左右的窗口比较稳定。只看数值还不够最好把原始波形和滤波后波形叠在一起对比重点看 R 波峰值位置是不是没动。如果滤波后 R 波位置和原始波形上的尖峰有明显偏移说明滤波器相位特性不合格。滤波顺序也很重要。常见做法是先去基线漂移再陷波工频最后做一个很轻的低通平滑。反过来先陷波往往会让基线漂移段产生额外振荡后续高通处理又有可能把这些振荡放大。先处理低频再处理高频逻辑上更顺。4. R波识别与心率计算从阈值检测到RR间期4.1 为什么只找最大值会翻车很多人一开始实现 R波识别就是np.argmax(ecg)找最大的峰。这个写法在安静、波形规范的数据上也许能用但只要换一段数据就很容易翻车。R 波不一定都是正向的某些导联或电极贴反时R 波是负向尖峰找全局最大值会直接找到高耸的 T 波或基线毛刺上去。即使 R 波是正向的心电幅值也会随呼吸、运动、导联接触状态变化固定阈值对幅度变化非常不敏感。R波识别真正依赖的是 QRS 波的局部特征它比 P 波、T 波更陡持续时间更短且在正常心电周期里有一个生理上的不应期。R 波识别算法可以很简单但不能只靠一个全局最大值至少要考虑局部窗口、阈值基准和不应期三个因素。4.2 一个能跑的R波检测流程带方向参数、局部窗口和不应期下面这段是我的常见做法先确定 R 波极性再用局部窗口找极大值最后用不应期跳过相邻的误检。代码不做任何机器学习只依赖信号滤波后的形态。import numpy as np def detect_r_peaks(ecg, fs, refractory0.2, threshold_ratio0.4, r_polarity1): # r_polarity1 表示R波向上-1表示R波向下 x ecg if r_polarity 0 else -ecg x x - np.median(x) # 用99分位数作为幅值基准避免单个噪声点把阈值抬得太高 threshold threshold_ratio * np.percentile(x, 99) window int(fs * 0.05) # 50ms局部窗口 peaks [] i 0 while i len(x): left max(0, i - window) right min(len(x), i window 1) local_max np.argmax(x[left:right]) left if x[local_max] threshold: peaks.append(local_max) # R波出现后进入不应期跳过生理上不可能出现下一次QRS的时间段 i local_max int(fs * refractory) else: i 1 return np.array(peaks, dtypeint) peaks detect_r_peaks(ecg_clean, fs, refractory0.2, threshold_ratio0.4)这段代码的关键参数有三个。threshold_ratio决定检测灵敏度取 0.4 表示阈值为信号 99 分位幅值的 40%如果波形质量差可以降到 0.3但要接受更多误检的风险。refractory是不应期默认 0.2 秒对应人体生理上每分钟最多 300 次心跳的上限可以在 0.18~0.25 秒之间调。window取 50 ms用来确定局部最值点太短会把一个 QRS 波分裂成多个候选峰太长会把 T 波并进来。r_polarity参数很容易被忽略但它很值得单独提一下。处理 II 导联时 R 波通常向上r_polarity1处理某些小导联或反接数据时把整个信号取负再走同一套最大值检测逻辑就避免了负向 R 波检测不到的问题。4.3 从R波位置到心率计算RR间期、瞬时心率和平均心率R 波位置确定后心率计算只剩简单算术。相邻两个 R 波的时间间隔叫 RR 间期单位通常用秒。瞬时心率就是 60 除以 RR 间期。平均心率不应该直接对瞬时心率取平均而应该先对 RR 间期取中位数再换算成每分钟心跳数。这样做能减少个别异常 RR 间期对平均值的冲击。rr np.diff(peaks) / fs # 单位秒 instant_hr 60.0 / rr # 每个RR间期对应的瞬时心率 # 用生理边界过滤明显异常值静息数据一般取0.4~1.5秒 valid_mask (rr 0.4) (rr 1.5) valid_rr rr[valid_mask] mean_rr np.median(valid_rr) avg_hr 60.0 / mean_rr sdnn np.std(valid_rr) rmssd np.sqrt(np.mean(np.diff(valid_rr) ** 2)) print(检测到R波数量:, len(peaks)) print(平均心率: {:.1f} bpm.format(avg_hr)) print(SDNN: {:.3f} s, RMSSD: {:.3f} s.format(sdnn, rmssd))valid_mask过滤的是明显不符合生理规律的 RR 间期比如误检造成的 0.1 秒或者漏检造成的 2 秒。0.4~1.5 秒对应 40~150 bpm这是一个偏保守的静息范围。如果你处理的是运动负荷或儿童数据需要把范围放宽比如 0.3~1.5 秒否则算法会把 200 bpm 附近的心率直接过滤掉。SDNN是所有正常 RR 间期的标准差RMSSD是相邻正常 RR 间期差值的均方根。这两个指标是心率变异性分析里最常用的基础值。如果valid_rr里留下的点数太少这两个指标几乎不可信至少要有 30 段有效 RR 间期再谈 HRV。5. 心电算法避坑指南5个常见的翻车点与排查方法5.1 采样率算错滤波器参数全对结果全乱现象滤波后 R 波形态变得很钝或者原本 50 Hz 的工频干扰依然存在R波识别结果时好时坏。原因心电数据的时间列不是严格的等差数列有人用相邻时间差的均值算采样率某个时间戳异常跳动直接把fs拉偏更多时候是直接从配置里抄了一个没核对过的采样率。滤波器的截止频率是相对采样率归一化计算的fs一旦有偏差真实的滤波频带就跟着错了。解决读数据后一定用np.diff(time)核对时间间隔取中位数不要取均值。如果数据没有时间列要问清楚采集设备实际输出频率不能根据文件名猜。5.2 滤波顺序反了陷波器先跑波形出现振铃现象原始波形本来只是基线漂移加一点工频滤波后 QRS 波前后出现明显的来回振荡R 波峰值位置也变了。原因先做 50 Hz 陷波再做高通去基线高通会把陷波器在信号突变处产生的振铃进一步放大。陷波器本质是窄带滤波在 R 波这种陡峭波形上容易产生瞬态振荡。解决把滤波顺序固定为先去基线漂移再陷波工频。如果滤波后还在某些到访处出现振铃试着把高通阶数从 2 降到 1或者把陷波 Q 值从 30 调到 20。不要同时叠加太多滤波器每加一级就是一次信号变形。5.3 固定阈值碰上幅度变化安静段正常活动段大量漏检现象算法在静息数据上准确率很高一换到走路、握拳、深呼吸的数据R 波就漏掉一大半心率显示频繁跳变。原因阈值是按整段信号的固定百分比设的。静息段幅值大运动段基线起伏和噪声让有效 R 波幅值变小固定阈值落在 R 波峰值之上自然检不到。心电信号是非平稳的幅度会持续变化。解决把阈值计算从“整段一次”改成“滑动窗口更新”。常见做法是每次检测到有效 R 波后用最近若干秒的最大幅值重新计算阈值或者把 R 波幅值做一个滑动中位数作为基准。同时对静息、深呼吸、轻微运动三段数据分别做验证不要让一段数据的表现代表整体。5.4 负向R波和T波过高最大值检测把T波当R波现象R 波位置检测结果里有一半落在 T 波上相邻 RR 间期忽长忽短心率计算结果完全不规律。原因没有考虑导联极性。某些导联的 T 波幅值比 R 波还高或者 R 波本身是负向峰算法仍以正向最大值作为检测目标。T 波虽然比 R 波慢但局部幅值已经超过了阈值。解决先用肉眼或自动方式确认 R 波方向r_polarity参数设为 -1 就对信号取负再检测。同时把不应期设到 0.2 秒左右T 波落在 R 波后约 0.3 秒处不应期天然挡掉一部分 T 波误检。如果 T 波实在太高可以考虑对信号做一阶差分利用 R 波的陡峭边界定位而不是直接看幅值。5.5 心率计算直接对瞬时心率取平均一个错检把结果拉到 180现象明明只有两三个地方误检了打印出来的平均心率却 160、180和真实心率完全对不上。原因一个 0.3 秒的误检 RR 间期换算成瞬时心率就是 200 bpm对这个值取平均会把整体均值拉高。更常见的是漏检产生的超长 RR 间期算出的瞬时心率低到 20直接把平均值拉低。解决先过滤异常 RR 间期再取中位数不要取平均值。中位数对异常值不敏感这是心率计算里最重要的一步。另外打印结果时顺手打印 R 波数量和有效 RR 数量如果有效 RR 数量远小于 R 波数量说明滤波或检测环节有系统性问题不是心率算法的问题。6. 从“能跑”变成“能信”用对照标注做R波检测验证与参数微调R 波检测写出来很简单真正难的是证明它可靠。我的习惯是用一小段人工标注的 R 波位置做验证跑一个灵敏度sensitivity和阳性预测率positive predictive value再决定要不要改参数。没有量化指标所有“效果还行”都是主观感受。def evaluate_detector(ref_peaks, det_peaks, fs, tolerance0.15): tol int(tolerance * fs) # 先用广播计算每个检测峰和所有参考峰的距离 diff np.abs(det_peaks[:, None] - ref_peaks[None, :]) matched_ref np.zeros(len(ref_peaks), dtypebool) matched_det np.zeros(len(det_peaks), dtypebool) for i in range(len(det_peaks)): hit diff[i] tol if np.any(hit): # 只把最近的那个参考峰标记为匹配 j np.argmax(diff[i]) if False else np.where(hit)[0][0] if not matched_ref[j]: matched_ref[j] True matched_det[i] True tp matched_ref.sum() fn len(ref_peaks) - tp fp len(det_peaks) - matched_det.sum() se tp / (tp fn) ppv tp / (tp fp) return se, ppv se, ppv evaluate_detector(ref_peaks, peaks, fs, tolerance0.15) print(灵敏度: {:.3f}, 阳性预测率: {:.3f}.format(se, ppv))这里的tolerance0.15表示检测峰和参考峰落在 ±150 ms 内就算匹配成功。这个容差偏宽松适合验证整体检测流程如果你想做更严格的逐拍对比可以缩到 ±50 ms。matched_ref保证一个参考峰只被匹配一次避免一个真实 R 波同时被多个候选峰反复匹配。拿到指标后参数调整就不是拍脑袋了。灵敏度偏低说明漏检多优先把threshold_ratio降下来或者把局部窗口从 50 ms 调大一点。阳性预测率偏低说明误检多优先把不应期调大或者把阈值提上去。如果两个指标都很差先回头看滤波不要死磕检测参数。我还有一个习惯验证数据至少分三段一段安静平躺一段深呼吸一段轻微运动。心电算法最怕的就是只在“最好的一段数据”上表现好换个场景就翻车。我现在拿到任何新数据第一件事永远是把采样率、滤波参数、R波检测参数写成一个独立配置区任何改动都重新跑一遍三段数据保存前后指标做对比。这个习惯让我少改了很多轮。希望帮到你。本文还有配套的精品资源点击获取

关于本文作者

来自尧图内容编辑团队

尧图内容编辑团队 内容团队

尧图内容编辑团队

本文由尧图网络内容编辑团队执笔。团队由资深项目经理、前端工程师与设计师组成,所有内容均来自亲手交付的真实项目,先讲清问题、再给出可落地的解法。尧图深耕北京网站建设十年,服务过京华建材集团、智造科技等各行业客户,把一线经验沉淀为可复用的行业观察。

  • 十年建站经验,覆盖建材、制造、服务、文创等
  • 项目经理把关选题与事实准确性
  • 工程师与设计师联合撰写专业细节
  • 统一编辑规范,保证文风与排版一致
  • 每月复盘转化数据,迭代选题方向

延伸阅读

相关资讯与近期热门内容

深度阅读推荐

建站决策前值得细读的三篇

网站改版的5个关键决策
2024-08-12

网站改版的5个关键决策

什么时候该改版、改到什么程度、如何避免流量掉光,京华建材集团改版复盘给出答案。

获取专属建站方案

看完文章,把您的行业与预算告诉我们,免费获取一份量身定制的官网建设方案与报价。

立即免费咨询