简介面向心电图算法研究与开发的Python工程包系统覆盖滤波、R波识别、心率计算、特征提取、心律失常分类、房颤与室颤室速识别以及伪差干扰研究并配有可视化模块和测试工程。压缩包共91个文件、约30.81MB包含可直接运行的Python源码、MIT-BIH心电数据文件、深度学习相关脚本、PDF专业论文与图片资料脚本涵盖数据读取与格式转换、信号去噪、峰值检测、分类模型训练等环节便于按模块复现算法流程。已有673人浏览学习。读者可借助完整项目掌握从心电数据读取、信号预处理到模型评估的完整链路对于需要实现心电滤波、R波定位或心律失常分类的开发者能够直接参考巴特沃兹滤波、小波变换、SVM/随机森林/神经网络等实现并验证效果节省从零搭建的时间。1. 心电算法为什么滤波、R波识别、心率计算要放在同一条流水线里把心电算法拆成滤波、R波识别、心率计算三个环节看起来是一条教科书式的标准流水线但真拿Python把它跑通时翻车点一个接一个采回来的信号带着基线漂移50Hz工频混在QRS波群里滤波参数没选好R波直接被磨没阈值设得不对T波又被当成R波心率瞬间翻倍。这篇文章就沿着这条流水线往下走先讲预处理阶段到底滤什么再给出R波识别从特征信号构建到自适应阈值的完整Python实现最后落到心率计算和实时窗口设计。适合正在写心电信号处理算法的开发者、做可穿戴设备算法验证的工程师也适合准备把离线算法迁到实时系统的人。篇幅里的示例基于numpy和scipy默认采样率250Hz。2. 心电预处理要滤什么三种噪声的频段划分与带通参数取舍2.1 三种噪声的频段决定了滤波策略心电信号的有效能量集中在哪里决定了带通滤波器的两个端点。体表心电的主要成分在0.5~35Hz之间QRS波群的斜率能量集中在5~20Hz这个区间保留得越完整后面的差分检测越有利。麻烦来自三路噪声呼吸和肢体运动造成的基线漂移能量集中在0.05~1Hz供电系统的工频干扰在50Hz或60Hz的单一频点及谐波肌电干扰的频带很宽从20Hz一直延伸到数百Hz。三种噪声频段基本错开所以预处理通常不是“一个低通走天下”而是先用带通把能量限制在心电主频段内再根据工频的实际残留决定要不要陷波。噪声类型常见频段不处理时对R波识别的干扰基线漂移0.05~1Hz把幅值阈值整体抬高或压低R波“忽大忽小”工频干扰50Hz或60Hz在QRS波群上叠细碎抖动差分后被放大成假峰肌电干扰20Hz以上高频毛刺容易被误检成额外的R波这个表格的用途不是让你背参数而是排查时可以对照如果某类错误突增先想对应频段有没有滤干净。比如R波识别突然多出一串间隔规律的“R峰”优先怀疑工频残留波形整体上下飘优先怀疑基线漂移还没有被压住。下面我把两次滤波的代码和参数选择说清楚。2.2 Butterworth带通滤波的Python实现与参数设置常见做法是先用Butterworth带通把信号框到0.5~35Hz。之所以选这个范围是因为低于0.5Hz的成分对R波识别基本没有贡献只会干扰阈值判断高于35Hz的成分里QRS能量已经很少保留太多反而让肌电噪声在差分阶段占到优势。阶数取4是比较稳的折中2阶的过渡带太宽工频和基线漂移压制得不彻底8阶虽然更陡但相位失真更明显用零相位滤波时边缘效应也被放大。import numpy as np from scipy import signal fs 250 # 采样率单位Hz f_low, f_high 0.5, 35 # 带通截止频率单位Hz order 4 # Butterworth阶数 b, a signal.butter(order, [f_low, f_high], btypebandpass, fsfs) ecg_filt signal.filtfilt(b, a, ecg_raw)代码里有两个位置不能照抄完事。第一是截止频率我一般会先用0.5~35Hz跑通再根据数据表现微调儿童或运动负荷数据可以放宽到40Hz但必须同步观察T波误检率如果工频干扰很重不要急着把上限抬高而是先加陷波。第二是filtfilt它做的是零相位滤波正向和反向各滤一遍幅度和相位都保得很好适合离线分析和算法验证代价是滞后处理不能直接用于实时。你可以在调试阶段固定用filtfilt最后做实时化时再换成lfilter并处理相位延迟两段代码不要混用。另外一个容易忽略的点滤波后要检查R波形态是否被削弱。心电算法的经验是带通滤波对QRS幅值的影响很小但如果截止频率取成30Hz甚至25Hz部分窄QRS的高频成分会被削掉识别的峰位会往后拖。所以代码里我习惯在滤波后打印一次np.max(np.abs(ecg_filt))与原始信号对比超过30%的衰减就要检查截止频率是不是卡得太紧。注意不要在一个工程里混用filtfilt和lfilter做结果对比相位特性不同改换滤波方式后必须重新验证R峰定位精度。2.3 50Hz/60Hz工频干扰什么时候必须加陷波带通滤波在50Hz或60Hz处已经有一定衰减但衰减量常常不够。尤其是采集环境里有强供电干扰时差分信号里会残留周期性的小波动被滑动积分聚成一条条“假包络”心率直接翻倍。判断要不要加陷波不要凭感觉我一般会在滤波后对差分信号做一次快速傅里叶变换看工频频点是否仍然凸出或者直接观察滤波后的波形QRS波群上是不是还有密集的锯齿。# 50Hz工频残留明显时用陷波器打掉Q值选30左右 b_notch, a_notch signal.iirnotch(50, 30, fs) ecg_notch signal.filtfilt(b_notch, a_notch, ecg_filt)陷波器的Q值决定凹口的宽窄。Q值越大凹口越窄对QRS波群的高频成分伤害越小但凹口太窄时如果工频实际在50.2Hz陷波点会偏掉削弱效果大打折扣。我一般先用Q30如果效果不干净再降Q值如果降Q后QRS形态明显变钝就要回到采集端去查接地和导联屏蔽而不是在后端继续堆滤波器。注意存在起搏器的患者数据里窄起搏脉冲含有丰富的宽频带能量工频陷波对它的影响不大真正要防的是第二道形态校验这一点放到第5章再说。3. R波识别用差分、平方和滑动积分把QRS波群从干扰里找出来3.1 幅值阈值为什么会在T波和基线漂移面前翻车拿到滤波后的信号最容易想到的检测方案是设一个幅值阈值超过就是R波。这个方案在静息、信号干净、R波明显高于T波的时候能用但真实采集数据里几乎处处是反例T波在某些导联上幅度可以接近R波基线漂移没滤干净时整段波形被抬起来幅值阈值会同时把T波和噪声包络判成心搏更麻烦的是同一个阈值换一个受试者、换一个导联位置就要重调无法形成稳定算法。这不是“阈值不聪明”而是幅值特征本身不够稳定。因此主流的做法是把R波最独特的斜率特征先提取出来再在特征信号上做阈值判断这也是Pan-Tompkins那一类方法的核心思路。不要把它当黑匣子接下来我把特征构建的每一步拆开。3.2 构建特征信号差分、平方、滑动积分R波与P波、T波相比最大的区别是斜率陡、变化快。第一步用一阶差分把斜率转成幅值变化QRS对应的差分值会远大于T波和P波第二步做平方既把差分变成正值又放大高斜率段的差异第三步做滑动积分在120~150ms窗口内把平方信号聚合成包络避免单个噪声尖峰直接触发阈值。这个三步变换做完QRS的特征从“一段陡峭的波形”变成“一个高而宽的包络”T波和基线残留会被压到很低的水平。# 差分突出斜率QRS波群的差分能量远高于T波和P波 diff_ecg np.diff(ecg_filt) # 平方统一方向同时放大高斜率段的差异 squared diff_ecg ** 2 # 滑动积分150ms窗口把局部能量聚合成包络 win int(0.12 * fs) kernel np.ones(win) / win integrated np.convolve(squared, kernel, modesame)这段代码里要留意两个参数。窗口长度win是0.12秒的采样点数对应临床上QRS宽度加一小段余量窗口太短会保留过多单点噪声太长会把相邻的T波也糊进来modesame让积分输出和输入等长但数组两端会有边缘误差精确定位时不要依赖前50ms和最后50ms的峰值。另外差分输出比滤波信号短一个点integrated的索引和原始信号存在一个采样点的偏移这个偏移在后续回原始信号精调时会被覆盖掉但如果直接在积分信号上取峰误差会留在结果里。3.3 自适应阈值与不应期R峰定位的完整实现特征信号构建好之后剩下的工作是两件定阈值、加不应期。定阈值不能只用一个全局最大值乘以系数因为信号质量会随时间变化早搏、导联松动都会让峰值幅度突变我一般会用一个缓慢更新的自适应阈值。不应期则是利用心脏的生理特性一次QRS之后200ms内不可能再有一次完整的QRS检测器在这段时间内跳过避免把同一个波群的次级峰误认为新心搏。# 用前1秒信号初始化阈值 threshold 0.3 * np.max(integrated[:fs]) refractory int(0.2 * fs) # 200ms生理不应期 r_peaks [] k 0 while k len(integrated): if integrated[k] threshold: # 在不应期内找特征信号最大值 end min(k refractory, len(integrated)) peak k np.argmax(integrated[k:end]) r_peaks.append(peak) k peak refractory # 自适应更新新阈值由旧阈值和当前峰共同决定 threshold 0.7 * threshold 0.3 * (0.3 * integrated[peak]) else: k 1 # 把候选峰映射回滤波信号在30ms窗口内找绝对最大值得到R峰精确位置 fine_peaks [] for p in r_peaks: lo max(0, p - int(0.03 * fs)) hi min(len(ecg_filt), p int(0.03 * fs)) seg ecg_filt[lo:hi] fine_peaks.append(lo np.argmax(np.abs(seg)))这个实现是完整可跑的但有两个地方需要按数据调。第一初始化阈值里的0.3倍来自“R波包络远高于噪声包络”的经验值如果数据噪声偏大可以先取0.4~0.5宁可漏检一些也不要让噪声触发一连串假峰。第二阈值更新式里0.7和0.3是权重新阈值不会因为一个异常峰大幅跳动如果要更快适应长时间漂移可以把0.3换成0.4但相应的误检风险也会上升。最后一步精调很重要差分和积分会使峰位产生一个采样点的偏移所以要在原始滤波信号附近重新找极值取绝对值策略是为了兼容R波向上和R波倒置两种情况。4. 心率计算从RR间期到8秒滚动窗口实时心率怎么做4.1 RR间期与瞬时心率的换算R峰位置序列确定了心率计算就只是换算问题。相邻两个R峰之间是RR间期单位取秒瞬时心率用60除以RR间期得到。代码很短但概念上容易出错第一个R峰没有前置间期所以n个R峰只能算出n-1个RR间期如果直接把60除以“相邻采样点差”而没有除以采样率算出来的数值会大fs倍。rr_samples np.diff(fine_peaks) # 相邻R峰的采样点间隔 rr_sec rr_samples / fs # 换算成秒 bpm_instant 60.0 / rr_sec # 瞬时心率瞬时心率的物理意义是“当前这一拍的心率”在屏幕上显示时会出现明显跳动。比如一次早搏之后短RR算出180bpm下一秒长RR又算出60bpm这种跳变会让体验很差。所以实际产品里几乎不会直接展示瞬时心率而是把它作为中间量用于心律失常分析和后续的滑动窗口统计。4.2 滑动窗口里的稳健心率用中位数而不是均值展示给用户的平均心率通常取一个滑动窗口内的RR间期统计值。窗口长度常见的是8秒或10秒窗口越短响应越快但越容易被单次异常拍带偏窗口越长越平稳但运动场景下心率变化的响应会显得迟钝。我的习惯是先定8秒再根据产品对响应时间的要求微调。window_sec 8 # 8秒滚动窗口 r_peak_buffer [] def push_rpeak(r_idx): r_peak_buffer.append(r_idx) # 丢掉窗口外的旧R峰 while r_peak_buffer and r_peak_buffer[0] r_idx - window_sec * fs: r_peak_buffer.pop(0) if len(r_peak_buffer) 2: rr np.diff(r_peak_buffer) / fs return 60.0 / np.median(rr) return None这段代码把每个新R峰的位置传进来函数返回当前窗口内的稳健心率。关键选择是用中位数而不是均值两个相邻RR间期如果夹了一个误检峰均值会被这个异常值明显拉高中位数在窗口内有多个正常间期时基本不受影响。如果窗口内只有两个RR间期中位数退化成均值所以至少要等积累到3个以上R峰再输出心率否则宁可显示“信号测量中”。4.3 漏检与误检的修复先过滤再回原始信号补峰漏检一个R波会让某个RR间期接近两倍心率看起来砍半误检一个假峰则相反心率直接翻倍。处理办法分三层。第一层是生理范围过滤RR间期落在0.3~2.0秒之间才保留对应心率30~200bpm超出这个范围的间期标记为可疑第二层用窗口内中位数兜底单点异常不会明显影响输出第三层才是真正的修复一般只在离线后处理里做。median_rr np.median(rr_sec) for i in range(1, len(fine_peaks)): if rr_sec[i - 1] 1.6 * median_rr: # 怀疑在相邻两峰之间漏检了一个R峰回去找 lo fine_peaks[i - 1] int(0.25 * fs) hi fine_peaks[i] - int(0.25 * fs) if hi - lo int(0.05 * fs): local_idx lo np.argmax(np.abs(ecg_filt[lo:hi])) # 补检前再比一次能量防止把T波补进来 if np.max(np.abs(ecg_filt[lo:hi])) 0.5 * median_rr_amplitude: fine_peaks.insert(i, local_idx)这段代码里的lo和hi刻意让开了前后各250ms避免重复检测到已经识别的R波补检后还要做能量阈值校验。这种修复在离线分析里有价值但在实时系统里要非常克制你无法确认“这个空白间隔里到底有没有一个被噪声淹没的心搏”所以更好的做法是标记异常段交给信号质量模块去提示用户而不是在界面上强制补出一个心率值。5. 心电算法踩坑清单导联脱落、滤波相位延迟、单位与采样率5.1 导联脱落现象是R波连续误检根因是信号质量失效现象某一段数据里R波检测器输出一串连续的“R峰”间隔极短心率疯狂跳动但对应波形里根本看不到正常的QRS形态。原因导联脱落或接触松动后输入不再是心电信号而是一段饱和平线加毛刺固定阈值被毛刺反复触发。解决在检测之前加信号质量判断质量不合格的片段直接不参与心率计算。可以用一个简单函数先粗筛正常心电信号的RMS在一个合理范围内过零率远低于纯噪声两者越界就判为坏段。真正上线时还要配合长时间无峰检测连续2秒没有新R峰就重置状态机避免把噪声尾巴当持续心搏处理。def check_signal_quality(seg): # 心电信号质量的粗筛阈值按设备量纲调整 rms np.sqrt(np.mean(seg ** 2)) zc np.sum(np.diff(np.sign(seg)) ! 0) / len(seg) return (0.02 rms 1.5) and (zc 0.4)5.2 滤波相位延迟R峰位置整体偏移几十毫秒现象把检测算法的R峰定位结果和心电标注对齐时发现整体偏移几十到一百毫秒心率能算准但R峰位置用来对齐心拍就不行。原因用lfilter做实时滤波时IIR滤波器会给不同频率分量引入不同相位延迟QRS主峰被整体平移。解决离线分析直接换成filtfilt零相位实时处理则用lfilter加固定延迟补偿或者在候选峰附近开一个30~50ms的搜索窗在原始信号上重新找极值。后一种做法对系统性偏移的抵抗很好因为搜索窗把误差框住了即使滤波器延迟略有变化也能覆盖。5.3 幅值单位mV和μV混用阈值一个都不对现象同一套分析流程换了一台采集设备后一个R波都检不出来或者阈值乱跳。原因有的设备输出mV有的输出μV两者差1000倍固定阈值、固定RMS门限在不同的量纲下完全是两回事。解决在数据入口统一做量纲归一化。我的做法是先打印信号的最大绝对值根据量级判断单位如果整段信号最大值在几百以上通常意味着μV在1左右则意味着mV归一化之后所有阈值和RMS门限才具有跨设备的可比性。这个单位问题很容易被忽略因为滤波代码本身不会报错只会让检测结果默默变坏。5.4 采样率变了窗口长度却没有跟着变现象采样率从250Hz换成500Hz后同一个检测器出现大量重复检测。原因代码里把不应期写成了固定采样点数比如refractory52在250Hz下是208ms在500Hz下只剩下104ms生理不应期被腰斩同样的问题也会出在滑动积分窗口上。解决所有时间相关参数统一用int(0.2 * fs)、int(0.12 * fs)这种写法只写时间不写点数。换采样率时拿着旧参数清单从头对一遍凡是出现“500个点”这种写法的窗口都要怀疑不要只改fs变量就以为万事大吉。5.5 起搏器脉冲窄而尖的假R波靠形态校验过滤现象带起搏器的患者数据里起搏脉冲幅值很高、宽度极窄被识别成额外心搏心率数值翻倍。原因起搏脉冲持续不到1ms在差分和平方变换后形成极强的窄峰积分窗口反而把它放大。解决在峰位定位后加形态校验统计R峰附近30ms内的信号上升时间和峰宽正常QRS至少要几十毫秒宽而起搏脉冲是瞬时的尖刺也可以比较候选峰与相邻R峰幅值超过正常比例太多就删除。形态校验要放在RR间期过滤之前因为起搏规律可能导致异常RR也“看起来很有规律”。6. 最后一步用标准数据库验证算法再把离线流程改成在线流式6.1 用公开数据库算灵敏度与阳性预测值算法写完后第一步不是继续调参而是找一份带有R峰标注的公开心电数据库把检测结果和标注逐点比较。常见的比对规则是检测点与标注点相差±150ms以内算正确。然后算两个指标灵敏度反映漏检阳性预测值反映误检。两个指标必须一起看只盯一个会掩盖问题。比如灵敏度很高但阳性预测值很低说明检测器靠“多报”提高了召回这种算法在临床应用里同样不达标。指标公式关注点Se敏感度TP / (TP FN)漏检PPV阳性预测值TP / (TP FP)误检平均定位误差检测点与标注点的差值统计R峰位置精度6.2 把离线流程改成在线流式的三个关键改动离线代码跑通只是第一步多数真实场景要求实时输出。在线化时我一般只改三个地方其他逻辑尽量不动。第一滤波从filtfilt换成lfilter因果滤波无法避免相位延迟要么做延迟补偿要么用候选峰回检抵消误差。第二阈值从“全段计算”改成在线更新检测器只维护信号峰和噪声峰两个估计值每检到一个候选峰就按小步长更新。第三检测循环改成状态机空闲态检测超阈值检测到候选峰后进入不应期不应期结束回到空闲态避免在线滑动里对同一个波群重复计数。我自己的血泪教训是所有离线调好的参数切到在线处理前一定要在新采样率、新导联位置的数据上重跑一遍不要因为波形看着干净就跳过验证。心电算法是输入质量决定上限的方向滤波和后处理能救回一部分坏数据但救不了一切。把验证流程固定下来每次改动都做一次灵敏度回归算法才能真正稳定交付。希望帮到你。本文还有配套的精品资源点击获取