简介这套基于Python的心电算法工程涵盖了滤波、R波检测、心率计算、特征提取、心律失常分类及房颤室速识别等关键环节从信号预处理到病理诊断形成完整链路适合生物医学工程研究者、心电信号处理初学者以及医疗AI开发人员学习参考。压缩包共91个文件主要包含17个Python脚本、16组心电记录数据dat/hea/atr、深度学习示例、可视化图片与PDF论文等整体大小约30.81MB目录结构按data、ecg、docs、tests等模块划分便于检索。该资源目前已有673人学习下载。除核心算法外还提供MIT-BIH数据库读取工具、小波变换与QRS检测等经典方法的演示代码、测试工程以及相关参考文献可帮助读者复现实验、理解原理并快速搭建自己的心电分析流程。此外资源还包含伪差干扰研究案例与可视化展示脚本覆盖了心电图分析从数据读到特征研究的常见环节是一份可直接用于教学和课题开发的综合性参考资料。1. 用Python处理心电信号R波检测为什么是第一道门槛拿到一段心电信号用Python写滤波、R波识别和心率计算这串算法第一道门槛几乎都在R波检测上。R波是QRS波群里幅度最大、形态最陡的波形也是后续心率、心律 variability 分析的地基R波位置错一个RR间期就错一次心率结果立刻失真。这篇文章把「滤波 → R波识别 → 心率计算」这条经典链路完整拆开给出可直接运行的代码、每个滤波器的参数设计理由以及我在实际数据上踩过的坑。适合刚接触心电信号处理、想用Python写一套能跑通的算法的开发者也适合被基线漂移和T波误检磨到没脾气的熟手。读完你至少能自己搭出一个可复现的检测流程而不是对着波形图发愣。2. 心电信号滤波噪声长什么样滤波器参数怎么设2.1 心电信号里的三类噪声先分清再动手心电信号是毫伏级的生物电信号采集过程中混入的噪声往往比信号本身还“精神”。我见过不少新手拿到数据后先上一个50Hz陷波发现波形更乱了原因就是没分清噪声类型。按频段分心电噪声主要就三类基线漂移、工频干扰、肌电干扰。基线漂移来自呼吸、电极与皮肤接触阻抗变化频率集中在0.05Hz到0.5Hz表现为整条波形在缓慢上下“游泳”。这种噪声会抬高或者压低R波的绝对幅度如果不处理固定阈值的R波检测器会在吸气段和呼气段表现出完全不同的灵敏度。工频干扰是50Hz部分地区60Hz的电源耦合噪声在波形上看是一条叠在ECG上的细密锯齿幅度可能达到R波的几分之一。肌电干扰来自肌肉收缩频谱很宽从几赫兹一直延伸到几百赫兹形态上是高频毛刺会直接影响差分运算的结果。搞清楚噪声频段之后滤波策略就很清晰了基线漂移用高通或中值滤波压掉工频干扰用窄带陷波肌电等高频噪声靠低通或带通限制。R波检测场景下最常用的是一个5Hz到25Hz左右的带通滤波器它一次性完成“去基线漂移 去掉大部分肌电”两件事。但如果你后续还要做ST段分析这个带通范围就不能用因为ST段能量主要落在0.5Hz到5Hz之间被滤掉就什么都看不到了。2.2 用scipy搭带通滤波器截止频率和阶数的取舍我用的是scipy.signal里的butter和filtfilt组合。先写一个通用的带通滤波函数import numpy as np from scipy import signal def butter_bandpass(lowcut, highcut, fs, order4): 设计Butterworth带通滤波器返回分子分母系数 nyq 0.5 * fs low lowcut / nyq high highcut / nyq b, a signal.butter(order, [low, high], btypeband) return b, a def bandpass_filter(data, lowcut, highcut, fs): 零相移带通滤波返回与输入等长的信号 b, a butter_bandpass(lowcut, highcut, fs) # padlen保证滤波起始段不被边缘效应吃掉 y signal.filtfilt(b, a, data, padlen3 * (len(b) - 1)) return y这里我特意用了filtfilt而不是lfilter。两者的差别很关键lfilter是因果滤波会给波形引入相移R波峰值位置会被推移几十毫秒filtfilt先正向后反向各滤一遍相位延迟相互抵消波形形态不会畸变R峰位置基本保真。做R波定位的滤波器零相移是硬要求。带通的范围怎么选我一般按采样率和用途分两类。做R波检测低端设5Hz、高端设20到25Hz这个频段既保留了QRS波群的主要能量又大幅度削弱T波和P波这些低频成分。做全波形分析或者ST段分析就放宽到0.5Hz到50Hz但这时候R波检测器如果直接跑在这个信号上T波误检率会明显上升。阶数方面4阶是性价比比较高的选择8阶的频带更陡峭但数值误差累积更明显滤波后波形可能在首尾出现更大的振铃。传统教材里常提到的0.5Hz到35Hz全频带滤波我实际做下来发现它更适合给人眼看波形而不是给算法做R波检测。算法需要的不是“看着舒服”而是“特征突出”所以检测专用的滤波器通常比显示用的滤波器更窄。这是新手最容易转不过来弯的地方。2.3 工频陷波和基线漂移去除一个都不能省带通滤波能解决大部分问题但工频干扰是窄带信号带通只能衰减它的一部分剩下来的残余还会影响差分步骤。更稳妥的做法是对50Hz单独做陷波。用iirnotch设计一个Q值可控的陷波器def notch_filter(data, fs, freq50.0, Q30.0): 50Hz工频陷波Q控制陷波带宽 b, a signal.iirnotch(freq, Q, fs) y signal.filtfilt(b, a, data) return yQ值是个需要调的参数我把它单独拿出来说。Q越大陷波带宽越窄对50Hz附近的QRS频率成分伤害越小但数字滤波器的频响特性会在窄带上变得更敏感略微的频偏就可能导致陷波失效。实测里Q取25到35之间的表现比较均衡——既能压住50Hz的锯齿又不会明显吃掉QRS波群的边沿。如果你发现陷波后波形出现“凹陷”大概率是Q值设得太低陷波带宽把49Hz到51Hz之外的有用成分一起削了。基线漂移的处理可以直接用高通滤波器也可以用中值滤波。中值滤波的优点是能适应非线性漂移。一个典型的做法是取1秒左右窗口的中值作为基线估计再从原始信号里减掉def remove_baseline(data, fs, window_sec1.0): 中值滤波估计基线并去除窗口必须为奇数 window max(3, int(fs * window_sec) | 1) # 保证奇数 baseline signal.medfilt(data, kernel_sizewindow) return data - baseline窗口大小直接影响处理效果窗口太大基线估计跟不上缓慢漂移窗口太小会把呼吸引起的低频成分当成信号保留下来。我一般在1秒到1.5秒之间试以滤波后整段波形没有明显上下浮动为准。整条预处理链路我通常这样串先去基线漂移再陷波最后带通。顺序有讲究——中值滤波会改变信号的基线状态如果放在带通后面会把滤波后的直流偏移又带回来陷波器对直流偏置不敏感但会放大数值误差放在中值之后更稳。最终的输出给到R波检测器之前顺手看一眼滤波后信号的峰值是否出现了毛刺这一步能省掉后面大量排查时间。3. R波识别从Pan-Tompkins原理到可运行的检测器3.1 为什么不直接用find_peaks很多人拿到滤波后的信号第一反应就是scipy.signal.find_peaks找局部极大值。这个思路在安静数据上确实能跑但换到带噪声的真实采集数据问题立刻暴露。find_peaks只能基于幅度阈值和最小间距工作它对T波和高耸的P波没有区分能力一旦T波的峰值超过了阈值就会产生误检。另一个问题是形态——R波的特征不仅是“高”更是“陡”。经典的Pan-Tompkins算法提供了一个更可靠的思路用带通滤波突出QRS频段用差分运算放大R波的斜率特征用平方运算让大斜率信号进一步突出最后用滑动窗口积分把分散的QRS能量聚合成一个单峰。这样处理后的信号R波对应位置是一个尖锐而宽阔的峰而T波和基线噪声被压得很低阈值分割变得容易得多。需要说明的是现在的实践中很少逐字照搬原始论文的参数因为它基于360Hz采样率设计直接套到250Hz或500Hz数据上会有偏差。更常见的做法是保留它“差分→平方→积分→自适应阈值”的骨架按自己的采样率重算窗口长度和阈值系数。这样做出来的检测器对多种采样率的数据都有稳定的表现。3.2 检测主流程滤波、差分、平方、滑动窗口积分下面这个函数是完整可运行的R波检测器输入原始ECG和采样率输出R峰位置索引。核心步骤都用注释标出import numpy as np from scipy import signal def preprocess_ecg(ecg, fs): 完整预处理基线漂移去除 - 带通 - 陷波 ecg remove_baseline(ecg, fs, window_sec1.2) ecg notch_filter(ecg, fs, freq50.0, Q30.0) ecg bandpass_filter(ecg, 5, 25, fs, order4) return ecg def detect_r_peaks(ecg, fs): 基于差分-积分法的R波检测返回R峰位置索引 x preprocess_ecg(ecg, fs) # 一阶差分近似斜率QRS段的斜率远大于T波和P波 diff np.diff(x) # 平方运算放大斜率差异小斜率成分进一步被压缩 squared diff ** 2 # 滑动窗口积分窗口长度约120ms覆盖一个典型QRS宽度 win_len max(1, int(0.12 * fs)) kernel np.ones(win_len) / win_len integrated np.convolve(squared, kernel, modesame) # 前2秒用于初始化阈值 init_len min(int(2 * fs), len(integrated)) init_block integrated[:init_len] signal_peak np.max(init_block) * 0.6 noise_peak np.mean(init_block) refractory int(0.3 * fs) # 不应期心率上限300bpm last_peak -refractory r_peaks [] for i in range(1, len(integrated) - 1): # 只检查局部极大值点避免在每个采样点上都做阈值判断 if not (integrated[i] integrated[i - 1] and integrated[i] integrated[i 1]): continue # 不应期内不重复检测 if i - last_peak refractory: continue if integrated[i] signal_peak: r_peaks.append(i) last_peak i signal_peak 0.7 * integrated[i] 0.3 * signal_peak else: noise_peak 0.7 * integrated[i] 0.3 * noise_peak # 阈值随信号和噪声估计动态更新 threshold noise_peak 0.25 * (signal_peak - noise_peak) # Python循环里直接访问全局变量比较慢这里只在每次局部极大值时做阈值更新 return np.array(r_peaks)这段代码有几个值得掰开讲的点。第一个是滑窗积分长度0.12秒乘以采样率对应QRS波群的典型宽度。窗口太短积分后的峰可能分裂成几个小峰窗口太长会把相邻的两个R波糊在一起。第二个是局部极大值检查它保证了每个峰只触发一次检测不会被积分波形上的次级波动重复触发。第三个是threshold的计算放在局部极大值出现之后这样阈值永远反映的是“当前最近的信号峰和噪声水平”。代码里的0.25是Pan-Tompkins风格的经验系数。0.25意味着阈值更靠近噪声估计灵敏度偏高如果想压低误检可以把0.25调成0.4甚至0.5。但调高阈值的同时会牺牲对低幅度R波的检出率这是灵敏度与精确度的经典取舍没有绝对正确的值。3.3 阈值自适应和不应期把误检率压下去固定阈值在实际数据上几乎撑不过一分钟。原因很简单心电信号的幅度会由于呼吸、体位变化、电极接触状态而缓慢变化一个在前10秒工作正常的固定阈值可能在30秒后因为信号整体变弱而一个R波都检不到。自适应阈值的本质就是让检测器持续跟踪信号幅度的变化。我用的更新公式很直接。被判定为R波的峰值用来更新signal_peak被判定为噪声的峰值用来更新noise_peak两边的更新系数一个是0.7对0.3形成“历史占主导、当前微调”的平滑效果。阈值取噪声估计与信号估计之间的差值加权。换句话说信号越强、噪声越弱阈值就越容易被拉开信号突然变弱时阈值也会在几个心跳周期内跟上。不应期是另一个关键环节。生理上人的心室不应期决定了两个正常R波之间至少隔约200毫秒换算成心率是300次/分钟。我不应期设0.3秒是因为剧烈运动或峰值心率场景下R波间距会缩短到300毫秒左右设成0.2秒会卡得太死设成0.3秒刚好在大多数场景下压制T波误检又能跟上高心率。如果你处理的是静息心电0.25秒就够处理运动负荷心电建议0.2秒。这里有一个常见的翻车点有些人把不应期直接用作find_peaks的distance参数但distance要求的是整数采样点数如果采样率是250Hz0.3秒对应75个点很多人忘掉换算设成distance0.3结果一个R波都检不出来。我吃过这个亏每次看到“检测结果为空”的报错第一反应就是去检查距离参数单位。4. 心率计算从RR间期到可展示的心率值4.1 从R峰位置到瞬时心率和平均心率R峰位置找到之后心率计算本身很简单但单位换算和异常处理才是真正拉开差距的地方。相邻两个R波之间的采样点数除以采样率得到的是以秒为单位的RR间期60除以它就是这一跳的瞬时心率。平均心率可以有多种定义有人用所有RR间期的均值有人用中位数还有人用总时长除以总心跳数。用哪种取决于你的场景。静息状态下RR间期相对稳定均值和中位数差别不大。但只要有漏检或误检均值会被极端值拉偏中位数则稳健得多。我的习惯是先算RR间期序列剔除异常值后再取均值这样既保留了中位数对离群值的容忍度又不会丢掉正常心跳间的细微差异。总时长除以总心跳数的做法适合长程统计但在短时间片段里容易把边缘的半个间期也算进去我不太推荐。直接看计算代码def compute_heart_rate(r_peaks, fs): 输入R峰索引数组和采样率返回 rr_intervals: 所有RR间期秒 inst_hr: 每个心跳的瞬时心率bpm avg_hr: 过滤异常间期后的平均心率bpm if len(r_peaks) 3: raise ValueError(R峰太少无法计算心率) rr np.diff(r_peaks) / fs # 秒 inst_hr 60.0 / rr # bpm # 用中位数做异常过滤而不是用均值 rr_median np.median(rr) lower 0.5 * rr_median upper 1.8 * rr_median valid (rr lower) (rr upper) avg_rr np.mean(rr[valid]) avg_hr 60.0 / avg_rr return rr, inst_hr, avg_hr过滤范围0.5到1.8倍中位数是怎么来的正常心律下相邻两次心跳的间期波动一般在±30%以内运动状态下波动会大一些。0.5倍以下的间期出现几乎可以断定是误检把T波或噪声当成R波1.8倍以上的间期出现基本都是漏检中间少跳了一个R波。这个范围对绝大多数心率计算场景够用而且实现极其简单不用引入复杂的统计模型。还有一个容易被忽略的点60.0 / rr在rr接近0时会得到巨大的瞬时心率。比如某个误检峰跟前面的R波只隔了0.1秒瞬时心率会显示600bpm。这显然不是真实生理值计算平均心率前不加过滤一个这样的离群点就能把均值拉飞。4.2 漏检和误检怎么暴露RR间期分布是照妖镜心率算完别急着画结论先把RR间期序列打出来看一遍。我每次调检测器参数都会先打印RR间期序列它比任何指标都直观。正常心律的RR间期序列应该是一条在均值附近小幅波动的横线如果中间出现了一个接近2倍的大值说明这里漏检了一次如果出现一个接近0.5倍的小值说明这里多检了一个波。漏检和误检对应的修正方式不一样。漏检时常见做法是把这个大间期等分成两个间期相当于在中间补齐一个R波。但这个操作要谨慎——如果漏检是因为T波幅度压过了阈值补进去的可能正好是把T波当成R波。误检时直接删除这个异常间期同时把对应的R峰索引从数组里去掉。修正完之后重新计算平均心率这个值才是能对外展示的。我习惯于把上述判断逻辑做成一个简单的表格方便在调试时对照RR间期特征最可能的原因建议处理接近正常值的0.5倍R峰误检多为T波或噪声删除对应R峰接近正常值的2倍R峰漏检结合波形确认后补齐连续多个间期忽长忽短检测器阈值或滤波参数不匹配回到第3章重新调参间期整体偏大且波动小采样率与R峰索引单位不一致检查fs是否与数据匹配单位问题看起来低级实际频率很高。有次某个开发者拿了360Hz的数据代码里却写fs250R峰索引除以250后所有间期被拉长了44%心率显示40多以为出了室性心律不齐最后发现只是采样率常量写错。R峰位置是采样点序号它离不开采样率这个参照系任何一步换算错后面全错。5. 心电算法避坑指南5个真实翻车场景与排查方法5.1 滤波后波形很干净但一个R波都检不到现象是滤波后的波形看起来平滑规整肉眼能清楚看到R波但检测器返回空数组。原因通常是带通滤波器的高通截止频率设太高了。有人为了压基线漂移把低端截止频率直接设到5Hz以上甚至10Hz结果QRS波群的能量被削掉一大半集成后的峰值比T波还低阈值初始化时把噪声当成了信号。解决方法是回看带通范围。R波检测的低端截止频率建议设在3到5Hz之间别超过8Hz。另外检查一下是不是用了lfilter——因果滤波的相位延迟会把R峰位置移动几十毫秒虽然不会导致完全检测不到但会让标记位置看起来明显偏斜。5.2 陷波器用完后QRS波形出现“凹陷”现象是陷波后的心电图上每个QRS波群的位置多了一个向内凹的缺口看起来像被咬了一小口。原因是陷波滤波器的Q值设太小陷波带宽太宽把50Hz附近的低频分量连着QRS的边沿一起衰减了。解决办法是把Q值从10调高到30再试。Q10时陷波带宽大约5Hz覆盖49Hz到51Hz但实际数字滤波的过渡带会向外延伸QRS波群中含有的40Hz到60Hz成分被误伤。Q30时带宽收窄到1.5Hz左右对QRS影响小得多。如果调整后锯齿纹仍明显优先检查数据里的工频是否为50Hz有些采集设备标注50Hz实际干扰在60Hz。5.3 T波被反复当成R波现象是检测结果里每隔两次正常R峰中间就多了一个峰位置刚好落在T波上。出现这种情况说明T波在滤波后的幅度接近甚至超过了R波。原因有两类一是带通范围太宽保留了T波的能量二是不应期设置得太短给了T波可乘之机。处理时先看带通高端把highcut从25Hz降回20Hz对T波有一定的压制作用。如果T波幅度依然很大则把不应期从0.3秒稍微延长比如0.35秒这会以损失极高心率检测能力为代价但对静息数据很有效。还有一个被低估的细节如果数据里T波本身比R波还高这个数据可能是单极导联或者导联接触异常先检查原始波形别急着调参数。5.4 心率结果忽高忽低平均心率没法看现象是打印出的瞬时心率在60、120、180之间跳动。看一眼RR间期序列就能定位问题——出现0.5倍中位数左右的短间期说明某个噪声峰或T波被算成了R波出现2倍左右的长间期说明某次R波没被检出来。解决时先确认滤波参数再检查阈值更新系数。如果阈值系数0.25太灵敏试着改成0.35或0.4误检会减少。漏检往往是因为初始阈值设置过高把初始的signal_peak从0.6倍调成0.5倍或者把noise_peak初始化方式从均值改为更小的分位数有助于把低幅度的R波也捞回来。5.5 filtfilt报padlen错误或首尾段信号异常现象是运行滤波时报错ValueError: The length of the input vector x must be greater than padlen或者滤波后波形两头明显摆动。原因是filtfilt需要对信号首尾做边缘填充输入信号太短或者padlen设置超过了信号长度。解决办法是检查信号长度是否足够至少要几百个采样点然后显式指定合理的padlen。注意padlen最小要求是3 * (len(b) - 1)我通常在此基础上加一点余量。如果数据本身只有一秒不到建议先做拼接或放弃这一段的滤波而不是强行缩短padlen。6. 进阶用评价指标验证你的心电检测器而不是看几张图6.1 用灵敏度和阳性预测值量化检测效果肉眼看图容易产生“看起来不错”的错觉。要判断检测器能不能上线需要两个量化指标灵敏度Sensitivity和阳性预测值PPV。灵敏度衡量的是“真实R波里被检出多少”阳性预测值衡量的是“检出的峰里有多少是真的”。两者都高才说明检测器既没漏检也没误检。def evaluate_detector(true_peaks, detected_peaks, fs, tolerance0.15): 用±150ms容差判定TP/FP/FN true_peaks、detected_peaks: R峰索引数组 tolerance: 容差秒 tol int(tolerance * fs) detected_set set() tp 0 for p in true_peaks: # 在容差范围内找最近的检测峰 near detected_peaks[(detected_peaks p - tol) (detected_peaks p tol)] if len(near) 0: tp 1 detected_set.add(int(near[0])) # 每个检测峰只匹配一次 fp len(set(detected_peaks) - detected_set) fn len(true_peaks) - tp sensitivity tp / (tp fn) if (tp fn) 0 else 0 ppv tp / (tp fp) if (tp fp) 0 else 0 return sensitivity, ppv, tp, fp, fn判定规则里容差取150毫秒是因为R峰目视标注本身存在几十毫秒的人为误差卡得太紧会把人工标注的偏差也算成算法错误。你会需要一份标注好的R峰位置作为参考标准。实际做的时候可以先在公开标注数据集上验证参数再去跑自己的数据。一个检测器如果灵敏度高但PPV低说明它“检得多但检不精”常见于阈值设得偏低。反之灵敏度低而PPV高说明检出质量高但漏检明显是阈值偏高的典型表现。理想的调整路径是先把灵敏度拉到95%以上再逐步提高阈值看着PPV曲线上升直到灵敏度掉到不可接受为止。6.2 可视化验证的三个实用技巧量化指标之外可视化能快速定位问题区间。我每次调完参数都会做三张图。第一张是原始ECG叠加检测到的R峰标记直接在信号上画圆点一眼就能看出漏检和误检集中在哪个波段。第二张是滤波前后的波形对比放大图截取三到五个心跳周期的片段确认滤波后的R波形态是否保留完整。第三张是RR间期序列图横轴为心跳序号纵轴为间期毫秒。这张图的诊断价值最高正常的RR间期应该在均值附近窄幅波动出现单点跳变说明该处有漏检或误检出现阶梯式变化说明信号质量在某段时间内发生了改变。我最早做这套检测器时花了两天排查一个R峰偶尔丢的问题最后发现就是带通滤波器的低端截止频率设到了1Hz呼吸波的低频能量漏了进来阈值被慢慢抬高。把截止频率改回5Hz后问题立即消失。现在的习惯是每次拿到心电数据先画RR间期序列再动手调滤波参数顺序反了就是白调。希望帮到你。本文还有配套的精品资源点击获取
