简介面向心电信号处理学习者和算法研究人员这份资源提供了一套完整的Python心电分析实现覆盖数字滤波、R波检测、心率计算、特征提取、心律失常分类以及房颤/室颤识别并包含可视化与测试工程。压缩包共109个文件以17个Python脚本、16组心电数据dat/hea/atr为主辅以PDF/Markdown说明、Jupyter Notebook示例和结果图片整体约34MB结构清晰。目前已有96人学习适合具备一定Python基础、希望系统掌握ECG处理全流程的开发者。资源源自网络分享仅用于学习交流使用时请注意版权声明。1. 一段裸心电波形到屏幕上的心率数之间发生了什么拿到一份几十秒的心电波形屏幕上跳动的那个心率数是怎么来的很多人以为心率是「数波峰」数出来的实际上一路走下来你会发现真正难的不是数数而是让每个R波都出现在它该在的位置。滤波、R波检测和心率计算这三个步骤几乎是所有心电分析算法从手环到Holter的公共地基自己用Python写一遍Pan-Tompkins链路比用现成库调接口更能建立对信号处理的直觉。如果你刚想入门心电信号处理、做穿戴设备算法或者只是被期末项目逼到Python面前下面的内容会把每一步的原理、参数和踩过的坑一起讲清楚看完你就能自己跑出一份可靠的心率曲线。2. 滤波选型基线漂移、工频干扰和肌电噪声分别用哪一招对付2.1 三种干扰和一段心电图里的频段分工心电信号本身只有毫伏级环境里的噪声比它大得多。先认识三个最常见的对手基线漂移、工频干扰和肌电噪声。基线漂移的频率在0.15到0.5赫兹之间来源是呼吸和电极片的微动表现就是整条波形像在水面上浮动如果不处理后续差分运算会把这种缓慢起伏当成大斜率变化阈值检测会一片混乱。工频干扰在50Hz左右放大器屏蔽没做好的时候波形上会叠一层细密花纹肉眼看起来就是一条加粗的线这种干扰幅度不大但会抬高信号能量基线。肌电噪声频率范围最宽从几十赫兹到几百赫兹都有肌肉紧张时波形像一团毛刺。数字滤波的做法本质上和硬件里RC滤波电路、π型滤波解决的是同一个问题只是把电阻电容换成了滤波器系数。做心电滤波前要先分清楚你的目标如果只是检测R波算心率QRS复合波的能量集中在5到20Hz取5-15Hz窄带就够了如果要做ST段分析之类的诊断任务低频必须保留到0.5Hz否则ST段抬高的临床特征会被一起滤掉。这个语义差别决定了两套完全不同的滤波参数先想明白再动手。2.2 中值滤波去基线漂移窗口按秒换算强制奇数去基线漂移我一般不用多项式拟合因为漂移的形态不是标准正弦波拟合阶数选不好反而把真实的ST段当漂移抠掉了。中值滤波更省心它在一个滑动窗口里取中位数对QRS这种短时突刺不敏感窗口内的中位数大概率取自基线附近于是剩下的baseline就是漂移本身原始信号减掉它即可。from scipy.signal import medfilt def remove_baseline(ecg, fs): # 窗口取0.5秒对应的采样点数并强制为奇数 kernel fs // 2 if kernel % 2 0: kernel 1 baseline medfilt(ecg, kernel_sizekernel) ecg_detrend ecg - baseline return ecg_detrend逻辑说明medfilt是排序滤波器输出每个位置邻域的中位数。QRS波群本身的宽度大约80到120毫秒在一个0.5秒的窗口里只占大约1/5剩余时间都是平缓的基线所以中位数能落在基线上而不是被R波拖走。参数说明kernel必须按采样率换算而不是写死整数。数据集可能是125Hz、250Hz、360Hz、500Hz同样0.5秒对应的采样点数完全不同。窗口越小跟踪漂移越灵敏比如患者说话动作大时漂移频率升高可以把窗口缩短到0.2秒但窗口小于0.2秒后中位数会开始包含QRS的骨架基线估计里带上心跳轮廓这部分要反复实验。medfilt要求kernel_size为奇数fs//2是偶数时记得加1这是最常见的低级翻车点。2.3 带通滤波与SG滤波一个用于检测一个用于打磨波形去完基线之后如果目标就是R波检测带通5-15Hz是最直接的路线。这个频带把P波、T波这些低频成分压掉保留R波的主要能量还能顺便压掉50Hz工频连专门的陷波器都不用做。如果目标是出一张好看的心电波形图还要保留P波T波带通范围放宽到0.5-40Hz再叠一层SG平滑。from scipy.signal import butter, filtfilt def bandpass_filter(ecg, fs, low5, high15, order2): nyq fs / 2.0 low low / nyq high high / nyq b, a butter(order, [low, high], btypeband) # filtfilt是零相位滤波R波位置不偏移 filtered filtfilt(b, a, ecg) return filtered逻辑说明butter返回的是数字滤波器的传递函数系数b和afiltfilt把信号正向滤一遍、反向滤一遍相位延迟互相抵消输出波形与原始信号在时间轴上一一对齐。这一点对R波检测很关键如果用lfilter之类的普通滤波器输出会整体延迟几十毫秒R波位置系统性偏移心率算出来就是带偏差的。参数说明order2对心电足够阶数越高过渡带越窄但瞬态响应变长、信号两端畸变更明显。low和高参数就是截止频率单位Hz注意代码里除以nyq做了归一化。检测场景5-15Hz是经验起点如果只做可视化把low改成0.5、high改成40即可。带通之后波形已经比较干净但带内还残留一些毛刺这时可以上SG滤波打磨一下from scipy.signal import savgol_filter def smooth_ecg(ecg_filt, fs): # 取100ms窗口配合R波宽度 win fs // 10 if win % 2 0: win 1 # polyorder3三次多项式局部拟合 ecg_smooth savgol_filter(ecg_filt, window_lengthwin, polyorder3) return ecg_smooth逻辑说明SG滤波在滑动窗口内做局部多项式拟合它的特点是在平滑噪声的同时保留波形的尖锐转折不像均值滤波会把R波峰顶磨圆。window_length取100毫秒基本贴合QRS宽度polyorder3是保守选择polyorder越高越保留细节但过拟合噪声的风险也越高。噪声类型频率范围处理手段关键参数基线漂移0.15-0.5 Hz中值滤波窗口0.2-0.5秒按fs换算取奇数工频干扰50 Hz国内带通压掉或陷波陷波带宽1-2Hz肌电噪声20 Hz以上带通SG平滑带通上限15-40HzSG窗口100msR波主能量约5-20 Hz带通5-15Hz检测用order23. R波检测从差分平方到滑动积分复现Pan-Tompkins的完整链路3.1 为什么直接阈值找峰不靠谱很多新手第一步写的是「np.max定位」或者「把信号大于某阈值的位置当R波」结果阈值调了一下午高一点漏检低一点把T波当R波。原因在于R波幅度个体差异太大0.5到2.5毫伏都有T波正常情况下只有R波的四分之一到三分之一但某些导联、某些患者的高尖T波幅度可以逼近R波。单纯看幅度根本分不清R波和T波必须借助一个更本质的特征——斜率变化率。R波的上升沿是整个心电周期中斜率最陡的部分这就是Pan-Tompkins算法的切入点。3.2 差分、平方、滑动窗口积分把R波变成脉冲Pan-Tompkins的思路是先对滤波后的信号做一阶差分把斜率转换成幅值再对差分结果做平方让负向的S波也变成正值同时放大R波和高频噪声之间的差距最后做滑动窗口积分把QRS周围的多个差分峰融合成一个平滑平台供阈值检测使用。import numpy as np def moving_integral(ecg_filt, fs): # 一阶差分R波斜率最陡差分值最大 diff_x np.diff(ecg_filt, prependecg_filt[0]) # 平方S波方向翻正同时拉开与T波的差距 squared diff_x ** 2 # 150ms积分窗口恰好覆盖QRS宽度 win int(0.15 * fs) window np.ones(win) / win # 卷积实现滑动平均 mw np.convolve(squared, window, modesame) return mw逻辑说明三个运算各司其职。差分后R波附近数值最大T波因为斜率平缓差分值小了一个量级平方进一步放大这种差距比如原本差6倍平方之后差36倍滑动积分的目的是把QRS内部因R波上升、下降产生的多个差分峰合并成一个单峰避免检测器在一个心动周期里被抖动触发多次。这里的卷积就是经典的滑动窗口滤波模型和烟雾传感器里用的滑动平均滤波算法是同一个运算。参数说明win取150毫秒对应生理上QRS波群的宽度这是Pan-Tompkins原始论文里反复验证过的值。如果采样率是360Hzwin就是54个采样点如果是500Hzwin就是75个采样点。卷积模式用same保证输出长度和输入一致否则后面索引对齐会乱。窗口过短会保留多个小峰过长会把下一个P波都融进来这两端都会导致检测位置偏移。3.3 自适应阈值与不应期动态变化的心律才测得住积分后的信号里R波对应一个大平台平台内部和外部都有噪声。固定阈值的问题是心率变化时信号能量跟着变运动状态下R波幅度猛增静息状态下低平固定阈值要么漏检要么误检。自适应阈值的做法是持续跟踪一段时间的信号水平每次检测到新峰就更新一次阈值让阈值跟随信号幅度走。不应期则是从生理上约定的最短心跳间隔正常情况下心室完成一次除极后需要一定时间才能再次激动200到300毫秒内不可能出现第二个真正的R波处于这个时间窗内的峰直接当作噪声处理。def detect_r_peaks(mw, fs, refractory0.25, th_ratio0.3): n len(mw) peaks [] # 用前8秒的99分位数初始化信号水平避免被单次突刺拉高 init_len min(int(8 * fs), n) signal_level np.percentile(mw[:init_len], 99) threshold th_ratio * signal_level refractory_n int(refractory * fs) last_peak -refractory_n pending -1 for i in range(n): # 超过1.5秒没检测到峰阈值放宽10%防漏检 if i - last_peak int(1.5 * fs): threshold * 0.9 if mw[i] threshold: # 进入峰值区间记录区间内最大值的位置 if pending -1 or mw[i] mw[pending]: pending i continue if pending ! -1: # 峰值区间结束判断这次峰值是否有效 if pending - last_peak refractory_n: # 有效R波更新信号水平并接受该峰 signal_level 0.875 * signal_level 0.125 * mw[pending] threshold th_ratio * signal_level peaks.append(pending) last_peak pending else: # 不应期内的峰当作噪声但也参与更新 signal_level 0.875 * signal_level 0.125 * mw[pending] threshold th_ratio * signal_level pending -1 return np.array(peaks)逻辑说明信号超过阈值时暂存为pending退出峰值区间时再判断是否为有效R波。信号水平和噪声水平共用一套递归更新公式新阈值的权重是0.875倍的旧水平加上0.125倍的当前峰值幅度这是Pan-Tompkins原始论文中的递归平均思路让阈值能平滑适应信号幅度的缓慢变化。如果1.5秒内没有产生任何有效峰大概率是阈值飘高了每次把阈值乘0.9做紧急下拉。参数说明refractory0.25秒是静息心电检测的默认值心率上限240bpm之内不会冲突如果专门做心律失常分析比如需要捕捉早搏需要把refractory降到0.15秒否则提前收缩会被不应期吞掉。th_ratio0.3表示阈值为信号水平的30%信号平稳时可以调到0.2提升灵敏度噪声偏大时调到0.4更稳妥。3.4 完整调用链路从原始信号到R波位置把前面的函数串起来就可以对一段原始心电完成整个R波检测流程。数据源常见做法是使用公开的MIT-BIH心律失常常用数据集通过wfdb库读取本地的dat文件import wfdb from scipy.signal import butter, filtfilt, medfilt # 读取一段心电数据sampto只取前100秒 record wfdb.rdrecord(105, sampto100 * record_info.fs) fs record.fs ecg record.p_signal[:, 0] # 三步预处理 kernel fs // 2 if kernel % 2 0: kernel 1 ecg_detrend ecg - medfilt(ecg, kernel_sizekernel) nyq fs / 2.0 b, a butter(2, [5 / nyq, 15 / nyq], btypeband) ecg_filt filtfilt(b, a, ecg_detrend) # R波检测 mw moving_integral(ecg_filt, fs) peaks detect_r_peaks(mw, fs) print(f采样率: {fs} Hz, 检测到 R 波: {len(peaks)} 个)逻辑说明wfdb.rdrecord是本课题最常用的数据读取入口record.p_signal取第一导联record.fs是数据集的采样率。这里用一个伪变量record_info.fs示意时间长度实际使用时直接替换成你确认过的采样率即可。预处理链路固定为「先去基线、再带通」的顺序不要反过来先带通会把基线漂移的低频成分截断再做中值滤波时基线估计反而更容易混入QRS能量。4. 心率计算RR间期、瞬时心率与平均心率三个口径别混用4.1 从R波位置到RR间期序列先除以采样率再做离群点滤波R波检测完成后得到的是一个个采样点索引转成时间间隔很简单相邻索引求差再除以采样率单位就是秒。但这里有个常见的脏数据问题漏检会让一个RR间期变成正常值的两倍误检会让间期缩短到正常值一半以下。直接把这些异常值丢进平均计算心率会明显偏离真实值。所以要先做一轮离群点滤波。def compute_rr(peaks, fs, min_rr0.4, max_rr1.5): # RR间期单位秒 rr np.diff(peaks) / fs # 过滤生理上不可能的间期400ms 或 1500ms 大概率是误检/漏检 valid_mask (rr min_rr) (rr max_rr) rr_valid rr[valid_mask] return rr_valid逻辑说明400毫秒对应150bpm心率如果用0.25秒不应期检测还能检到这个级别的间期说明已经到了剧烈运动场景1500毫秒对应40bpm低于这个心率通常需要特别谨慎地确认是否有漏检。这里的min_rr和max_rr是硬过滤不是数值修整过滤后的RR序列用于后续所有计算。参数说明如果你做的是新生儿心电正常心率上限可以到180bpm以上min_rr需要降到0.33秒左右如果是重症监护里的心动过缓患者max_rr可以放宽到2秒。参数跟着人群走不要一套打天下。4.2 瞬时心率与平均心率两个口径两个答案瞬时心率指每个RR间期对应的瞬时速率的倒数换算即60除以RR秒数单位bpm。注意心律不齐时瞬时心率天然会跳来跳去比如呼吸性窦性心律不齐可以让相邻两个瞬时心率差出十几这不是算法错了是生理事实。平均心率严格说有三种算口径第一种是60除以平均RR间期临床上最常用也叫平均心室率第二种是滑动窗口内数R波个数再折算即rate 窗口内R波数 * 60 / 窗口长度运动设备里更常用因为它反映的是一段时间内的平均速率而非单纯间期均值第三种是把瞬时心率序列做滑动平均常用于画趋势曲线。三种口径在心律绝对整齐时结果几乎相同心律不齐时差异明显输出时要说清楚你用的是哪种。def compute_hr(peaks, fs, window_secNone): rr np.diff(peaks) / fs rr rr[(rr 0.4) (rr 1.5)] # 口径一平均RR间期的倒数 hr_avg_beat 60.0 / np.mean(rr) if window_sec is None: # 口径二窗口内R波个数折算 duration (peaks[-1] - peaks[0]) / fs hr_avg_count len(peaks) * 60.0 / duration else: # 口径三滑动窗口均值 win_n int(window_sec * fs) centers np.arange(0, len(ecg), win_n // 2) # 步进为半窗 hr_avg_count [] for c in centers: left max(0, c - win_n // 2) right min(len(peaks), c win_n // 2) n_r np.sum((peaks left) (peaks right)) hr_avg_count.append(n_r * 60.0 / window_sec) hr_avg_count np.asarray(hr_avg_count) return hr_avg_beat, hr_avg_count逻辑说明function里的窗口滑动计算逻辑用了peaks作为索引严格说是按采样点滑动的center遍历的是np.arange(0, len(ecg), ...)其中ecg需要作为外部变量传入。我把这个实现拆得更清晰一点把原始信号长度作为参数传进来避免闭包引用外部变量导致不可复现的尴尬。改进后的版本在最后一段给出。参数说明window_sec一般取10秒对应临床上常用的心率观察窗口步进取半窗会产生重叠使得心率曲线更平滑。如果你只关心整段数据的平均心率直接用口径一最稳定。4.3 一个完整的可复用pipeline从原始信号到心率曲线把前面所有模块收拢成一个函数输入原始信号和采样率输出R波位置、RR间期序列和两类心率值。这是我在实际脚本里反复用到的骨架每次拿到新数据集先跑一遍它确认波形和R波标注没问题再做后续分析。def extract_hr(ecg, fs): # 1. 中值滤波去基线 kernel fs // 2 if kernel % 2 0: kernel 1 ecg_d ecg - medfilt(ecg, kernel_sizekernel) # 2. 带通5-15Hz nyq fs / 2 b, a butter(2, [5 / nyq, 15 / nyq], btypeband) ecg_f filtfilt(b, a, ecg_d) # 3. 差分-平方-滑动积分 diff_x np.diff(ecg_f, prependecg_f[0]) sq diff_x ** 2 win int(0.15 * fs) mw np.convolve(sq, np.ones(win) / win, modesame) # 4. R波检测 peaks detect_r_peaks(mw, fs) # 5. RR间期与平均心率 rr np.diff(peaks) / fs rr_valid rr[(rr 0.4) (rr 1.5)] hr 60.0 / np.mean(rr_valid) # 6. 瞬时心率曲线 hr_inst 60.0 / rr_valid return peaks, rr_valid, hr, hr_inst逻辑说明步骤1到4是固定的信号链路步骤5到6是心率的两个口径。注意hr_inst与rr_valid是一一对应的画图时横轴用对应R波的时间位置而不是均匀时间轴。拿到结果后至少做一次可视化验证这是整个流程中最容易被跳过的环节却往往能一秒钟暴露问题。python数据分析与可视化的画图能力在这里直接派上用场import matplotlib.pyplot as plt t np.arange(len(ecg)) / fs fig, ax plt.subplots(2, 1, figsize(12, 6), sharexTrue) # 上滤波后的心电R波标注 ax[0].plot(t, ecg_f, linewidth0.8) ax[0].plot(np.array(peaks) / fs, ecg_f[peaks], ro, markersize3) ax[0].set_title(ECG filtered with R-peak markers) # 下瞬时心率 ax[1].plot(np.array(peaks[1:len(hr_inst)1]) / fs, hr_inst, o-, markersize2) ax[1].set_ylabel(Instant HR (bpm)) plt.tight_layout() plt.show()逻辑说明上子图把R波位置画在滤波后心电波形上一眼就能看出每个R波是否被正确标出、有没有误标在T波或噪声上。下子图画瞬时心率序列如果曲线在某个点突然跳到200bpm以上说明那个位置大概率漏检或误检。这两张图画完算法的真实水平就心中有数了。5. 避坑指南五个让检测结果翻车的经典问题与解法5.1 中值滤波窗口偏小ST段被整段吃成直线现象预处理后的信号里R波和ST段看起来都变成了平直线心跳轮廓消失了。 原因中值滤波窗口设置得过小比如直接写死成50或100个采样点没有按采样率换算。窗口小于0.15秒时中值滤波把QRS当作局部毛刺全部滤除baseline估计里就包含了心跳本身的轮廓相减之后有用的信号反被抹掉了。 解决窗口按采样率换算为0.5秒对应的点数即fs//2保证窗口内有足够多的基线样本来计算中位数。修改参数后重新画图对比QRS恢复到正常形态即可。5.2 带通范围设反把R波主能量滤干净了现象经过带通滤波后的信号一片平滑R波几乎看不见检测结果为零或者极少量。 原因低通上限设得太低比如把high设成3Hz甚至更低。R波主要能量分布在5到20Hz3Hz的低通上限相当于把R波当成高频噪声滤掉了。这个错误往往发生在「直接把网上找的心电诊断滤波参数0.5-40Hz改成了检测用」思路没转过弯。 解决检测场景统一切到5-15Hz先不用管波形好不好看等R波检测稳定后再单独评估形态。如果上限超过20Hz肌电噪声会大量混入误检率直线上升。5.3 高尖T波被当成R波一个心动周期检出两个峰现象明明心电图是正常的窦性心律检测结果却显示心率是实际值的一倍瞬时心率曲线出现大量200bpm以上的点。 原因有些导联的T波较高且陡差分平方积分之后形成的局部峰值超过了当前阈值且间隔恰好大于不应期。正常T波确实不会被误检但高尖T波是心律失常患者的常见表现这才真正考验算法。 解决把不应期从0.25秒加长到0.3到0.35秒。代价是漏掉真实早搏的风险增加因此做心律失常分析时应期要单独调低这个场景只能做权衡。我一般同时输出平均RR间期如果一个心动周期内意外出现两个间隔低于0.3秒的峰先用人工标注确认T波形态再决定参数。5.4 初始阈值被前几秒的突刺抬高起步就漏检一大段现象算法跑完后开头几秒的R波一个都没有中段开始正常前部的瞬时心率空缺。 原因初始阈值用了整段信号的最大值或很长一段窗口的最大值。如果前8秒内有一个幅度极大的基线漂移或电极干扰尖峰max值被瞬间拉高阈值超过正常R波幅度导致起步区域漏检。 解决初始信号水平用前8秒窗口的99分位数而不是最大值。百分位数对单次突刺有天然免疫力8秒窗口保证了采样量的同时不至于太长。5.5 用了lfilter做带通R波位置整体偏移几十毫秒现象把自动检测的R波位置和原始波形叠加后发现标注点永远落在R波下降沿之后一点整体右移了固定几个采样点。平均心率差别不大但如果后续做HRV分析RR间期的系统误差会污染全部指标。 原因lfilter是非线性相位滤波输出信号的每个频率分量都有不同的相位延迟波峰位置发生系统性偏移。这里没有玄学就是滤波器本身的性质。 解决需要零相位滤波时一定用filtfilt正向反向各滤一遍相位互相抵消。代价是filtfilt不能用于实时流式处理因为它需要整段数据才能做反向滤波实时场景要换IIR滤波并在检测结果上补偿群延迟或者直接用纯因果的FIR带通加固定延迟修正。6. 进阶把算法实时化并用公开数据给检测结果打分如果要把这套流程跑在可穿戴设备或实时监护上离线整段滤波的做法不能原样照搬。常见做法是把信号按1秒窗口切成块步进0.5秒中值滤波和带通只针对当前块加历史缓冲做R波检测函数维护一个跨块更新的阈值和不应期计数。这样单次处理的延迟控制在几百毫秒内心率更新频率可以做到每秒2次RNN跑滑动窗口实时心率不是瓶颈瓶颈始终在R波检测的稳定性上。验证算法好坏的硬指标是灵敏度和阳性预测值。把检测结果和公开数据集的专家标注对齐在标注点50毫秒范围内存在检测点记真阳性TP检测到了但标注点不在范围内记假阳性FP漏掉了标注点记假阴性FN。灵敏度是TP除以TP加FN阳性预测值是TP除以TP加FP两个指标都大于95%才说明算法基本可靠。我第一次实现时只看了平均心率对得上就没做逐点匹配评估后来发现室性早搏患者的数据漏检惨烈平均心率却偏差不大从那以后每次跑新数据都先算这两个指标再谈别的。验证时手头没有现成标注集也可以用最笨的土办法把检测点画在波形图上人工抽查前30秒和最后30秒数一遍漏检误检数。这一看就能暴露出参数问题比任何自动化分数都直观。希望帮到你。本文还有配套的精品资源点击获取