简介这是一份面向心电算法初学者与生物医学信号处理爱好者的Python实战项目围绕心电图ECG分析全流程展开覆盖滤波、R波检测、心率计算、特征提取、心律失常分类以及房颤、室颤室速识别与伪差干扰研究等核心环节适合个人学习与算法验证。资源包共109个文件约34.06MB以17个py脚本、16个dat与16个hea数据文件、14个atr标注文件为主辅以zbak备份、xws、qrs、ipynb笔记本、pdf文档及jpg/png可视化图片构成从数据到代码再到结果的完整结构。项目借助巴特沃兹与卡尔曼滤波、小波变换、pandas、sklearn、tensorflow及matplotlib等工具实现信号预处理、R波定位、R-R间期心率计算、QRS幅度与QT间期特征提取、SVM与随机森林分类并提供可视化与测试工程。已有96人学习适合希望系统掌握心电信号处理链路、对照代码复现实验并理解病理识别思路的读者参考。1. 用Python实现心电算法从毫伏级噪声里捞出那个R波拿到一段真实的心电信号第一反应往往不是“怎么算心率”而是“这波形怎么这么脏”。50Hz工频串扰、基线漂移、肌电噪声叠在一起R波峰值可能只比噪声高一点点。用Python实现心电算法核心链路其实就三步滤波把噪声压下去R波检测把每次心跳的位置标出来心率计算把位置序列换算成可读的BPM。这套东西适合做可穿戴设备原型、健康监测App后端、或者单纯想把MIT-BIH数据集跑一遍验证算法的工程师。难点不在调库而在参数怎么定、边界怎么处理、异常波形怎么不翻车。下面按我实际做过的顺序拆开讲代码可以直接抄参数需要按你的采样率改。2. 滤波先搞清楚你要去掉什么再选滤波器心电信号的有效频段大约在0.05Hz到100Hz之间QRS波群能量集中在10Hz到40Hz。噪声来源不同滤波策略完全不同。工频干扰用陷波器基线漂移用高通肌电噪声用低通。很多人一上来就设计一个“万能带通”结果ST段被削平诊断信息全丢了。我一般会先看频谱再决定用几阶、什么类型。2.1 用Butterworth带通做第一道清理Butterworth的优点是通带平坦适合心电这种对幅值保真要求高的信号。下面这段代码用scipy.signal实现0.5Hz到45Hz的带通采样率按360Hz算MIT-BIH常见值。import numpy as np from scipy.signal import butter, filtfilt def bandpass_filter(signal, fs, lowcut0.5, highcut45.0, order3): signal: 一维心电数组单位mV fs: 采样率Hz lowcut/highcut: 通带边界Hz order: 滤波器阶数心电常用2-4 nyq 0.5 * fs low lowcut / nyq high highcut / nyq b, a butter(order, [low, high], btypeband) # filtfilt零相位滤波避免R波位置偏移 filtered filtfilt(b, a, signal) return filtered逻辑说明filtfilt做前后向两次滤波相位响应为零R波不会因为滤波而左右移动。参数上order不要超过4阶数越高振铃越明显QRS波后面容易拖出伪影。lowcut设0.5Hz是为了压基线漂移但如果你关注ST段可以降到0.05Hz代价是收敛变慢。highcut设45Hz能保留QRS主能量同时压掉大部分肌电噪声。注意filtfilt要求信号长度至少是滤波器阶数的3倍短片段会报错。2.2 工频陷波和滑动窗口滤波的取舍50Hz工频干扰在国内环境几乎躲不掉。常见做法是加一个IIR陷波器但陷波器Q值太高会引入振铃Q值太低又滤不干净。我一般用iirnotchQ设30左右。from scipy.signal import iirnotch, filtfilt def notch_filter(signal, fs, freq50.0, Q30.0): b, a iirnotch(freq, Q, fs) return filtfilt(b, a, signal)参数说明Q是品质因数30意味着-3dB带宽约1.67Hz刚好覆盖工频漂移范围。如果采样率低于100Hz50Hz陷波没意义因为已经混叠了。另外滑动窗口滤波在嵌入式端很常见但在Python离线处理里filtfilt效果更稳。滑动窗口滤波的窗口长度如果设成采样率的整数倍会引入周期性纹波这个坑后面会细说。2.3 基线漂移的两种处理路径基线漂移频率通常低于0.5Hz用高通滤波是最直接的。但高通滤波会改变ST段形态如果你要做缺血检测建议用中值滤波估计基线再减掉。具体做法用200ms窗口的中值滤波得到上包络再用600ms窗口得到下包络两者平均就是基线估计。from scipy.signal import medfilt def baseline_removal(signal, fs): window1 int(0.2 * fs) | 1 # 确保奇数 window2 int(0.6 * fs) | 1 baseline (medfilt(signal, window1) medfilt(signal, window2)) / 2 return signal - baseline逻辑说明中值滤波对尖峰不敏感不会把R波当成基线。window1取200ms是为了保留QRS细节window2取600ms是为了覆盖T波。两个包络平均后基线估计更平滑。注意medfilt在信号两端会有边缘效应前window2//2个点和后window2//2个点不可靠计算心率时要掐掉。3. R波检测Pan-Tompkins的Python实现与参数调优R波检测的经典方法是Pan-Tompkins1985年提出到现在依然是很多商用设备的底层逻辑。它的核心是带通滤波→微分→平方→滑动窗口积分→自适应阈值。听起来步骤多但每一步都有明确的物理意义。下面按顺序实现并给出每个环节的参数依据。3.1 微分、平方与滑动窗口积分微分增强QRS的斜率平方让所有值变正并放大高频成分滑动窗口积分把能量累加起来形成包络。def pan_tompkins_preprocess(filtered_signal, fs): # 微分五点差分 diff np.ediff1d(filtered_signal, to_begin0) # 平方 squared diff ** 2 # 滑动窗口积分窗口宽度150ms window_len int(0.15 * fs) integrated np.convolve(squared, np.ones(window_len)/window_len, modesame) return integrated参数说明微分用np.ediff1d是一阶差分计算量小。窗口宽度150ms是Pan-Tompkins原文推荐值对应QRS波群宽度。如果采样率是360Hz窗口长度54个点。modesame保证输出长度和输入一致但两端会有零填充检测时忽略前window_len个点。3.2 自适应阈值与R波定位阈值策略决定了检测的灵敏度和阳性预测率。Pan-Tompkins用两组阈值一组针对信号一组针对噪声每次检测后更新。def detect_r_peaks(integrated, fs, init_threshold_ratio0.25): # 初始化阈值 threshold init_threshold_ratio * np.max(integrated[:int(2*fs)]) r_peaks [] last_peak -int(0.2 * fs) # 不应期200ms for i in range(1, len(integrated)-1): if integrated[i] threshold and integrated[i] integrated[i-1] and integrated[i] integrated[i1]: if i - last_peak int(0.2 * fs): r_peaks.append(i) last_peak i # 更新阈值信号峰和噪声峰加权 threshold 0.125 * integrated[i] 0.875 * threshold return np.array(r_peaks)逻辑说明init_threshold_ratio取前2秒最大值的25%这是经验值。不应期200ms对应最大心率300BPM正常人不会超过。阈值更新用0.125和0.875的权重是Pan-Tompkins原文的递归公式。注意这段代码没有做回溯搜索如果漏检一个R波阈值会偏低导致后续误检。改进方法是在检测到峰后往前找200ms内是否有超过阈值一半的峰有就补上。3.3 在MIT-BIH上验证检测率拿MIT-BIH Arrhythmia Database的100号记录做测试采样率360Hz时长30分钟。用wfdb库读取跑上面的流程。import wfdb record wfdb.rdrecord(100, pn_dirmitdb) signal record.p_signal[:, 0] # 第一导联 fs record.fs filtered bandpass_filter(signal, fs) filtered notch_filter(filtered, fs) filtered baseline_removal(filtered, fs) integrated pan_tompkins_preprocess(filtered, fs) r_peaks detect_r_peaks(integrated, fs) # 与标注对比 annotation wfdb.rdann(100, atr, pn_dirmitdb) true_peaks annotation.sample # 匹配容差150ms tolerance int(0.15 * fs) tp sum(1 for r in r_peaks if np.min(np.abs(true_peaks - r)) tolerance) sensitivity tp / len(true_peaks) print(f灵敏度: {sensitivity:.3f})参数说明tolerance取150ms是因为R波检测允许一定偏移。灵敏度能到0.98以上算合格如果低于0.95检查滤波是否过度削峰或者阈值更新太快。注意wfdb需要联网下载数据第一次运行会慢。4. 心率计算从R波间期到BPM的四个边界坑拿到R波位置后心率计算看起来只是60*fs/RR但实际有四个坑异位搏动、漏检、误检、以及RR间期的生理范围限制。下面逐个处理。4.1 RR间期序列与异常值剔除先算相邻R波的间隔然后剔除不合理的值。正常RR间期在0.3s到2.0s之间对应30到200BPM。def compute_heart_rate(r_peaks, fs): rr_intervals np.diff(r_peaks) / fs # 单位秒 # 剔除生理不可能的值 valid (rr_intervals 0.3) (rr_intervals 2.0) rr_valid rr_intervals[valid] if len(rr_valid) 0: return 0, rr_intervals # 瞬时心率 hr_instant 60.0 / rr_valid # 平均心率用中位数抗异位搏动 hr_mean 60.0 / np.median(rr_valid) return hr_mean, hr_instant逻辑说明用中位数而不是均值是因为一个异位搏动会让均值偏移很大。valid范围设0.3到2.0秒覆盖了绝大多数生理状态。如果RR间期小于0.3秒大概率是误检了一个T波大于2.0秒可能是漏检。4.2 用滑动窗口做实时心率平滑如果是实时监测瞬时心率跳动太大需要平滑。滑动窗口取5个RR间期做加权平均。def smooth_heart_rate(hr_instant, window5): if len(hr_instant) window: return np.mean(hr_instant) weights np.arange(1, window1) weights weights / weights.sum() return np.convolve(hr_instant, weights, modevalid)参数说明window5对应大约5次心跳延迟约4秒。权重线性递增越近的心跳权重越高。注意modevalid会让输出长度减少window-1实时显示时前面几个点用均值补。4.3 异位搏动对心率的干扰与修正异位搏动早搏会让RR间期突然变短然后代偿间歇变长。如果直接算平均心率会虚高。修正方法是如果某个RR间期小于前一个的80%标记为异常用前后RR的均值替代。def correct_ectopic(rr_intervals): corrected rr_intervals.copy() for i in range(1, len(rr_intervals)-1): if rr_intervals[i] 0.8 * rr_intervals[i-1]: corrected[i] (rr_intervals[i-1] rr_intervals[i1]) / 2 return corrected逻辑说明0.8是经验阈值对应早搏提前20%以上。修正后心率更接近真实窦性心率。注意这个逻辑对房颤无效房颤的RR间期本身就是绝对不齐需要另外的算法。5. 避坑与排查五个让我熬夜的血泪教训5.1 滤波后R波幅值变小检测阈值失效现象带通滤波后R波峰值从1.5mV降到0.8mV固定阈值检测不到。原因Butterworth带通在45Hz处有衰减而R波的高频成分正好在30-40Hz。解决要么把highcut提到60Hz要么改用自适应阈值阈值按信号局部最大值动态调整。5.2 滑动窗口积分导致R波位置偏移现象检测到的R波位置比真实位置晚50ms。原因滑动窗口积分是对称的但np.convolve的modesame在两端补零实际输出有延迟。解决用modefull然后截取中间部分或者检测到峰后往前找微分信号的最大值位置。5.3 基线去除把ST段一起削掉现象中值滤波去基线后ST段抬高被抹平。原因window2取600ms时如果心率快T波和P波被当成基线。解决根据心率动态调整窗口心率大于100BPM时window2降到400ms。5.4 漏检一个R波导致心率翻倍现象平均心率从70跳到140。原因漏检后RR间期变成两个正常间期之和60/RR算出来是原来的一半。解决在心率计算前检查RR间期是否接近前一个的两倍如果是插入一个虚拟R波。5.5 工频陷波引入振铃误检T波现象陷波后QRS波后面出现振荡被误检为R波。原因IIR陷波器Q值太高群延迟非线性。解决降低Q值到15或者改用FIR陷波器代价是过渡带变宽。6. 进阶技巧用模板匹配提升复杂波形下的检测率Pan-Tompkins在正常窦性心律上表现很好但遇到室性早搏、束支传导阻滞时QRS形态变化大阈值法容易漏检。我一般会加一层模板匹配先取前10个检测到的R波对齐后平均得到一个模板然后对每个候选峰计算与模板的相关系数低于0.7的丢弃。def template_matching(signal, r_peaks, fs, template_len0.2): half int(template_len * fs / 2) templates [] for r in r_peaks[:10]: if r-half 0 and rhalf len(signal): templates.append(signal[r-half:rhalf]) if len(templates) 0: return r_peaks template np.mean(templates, axis0) template template - np.mean(template) good_peaks [] for r in r_peaks: if r-half 0 and rhalf len(signal): segment signal[r-half:rhalf] segment segment - np.mean(segment) corr np.corrcoef(template, segment)[0, 1] if corr 0.7: good_peaks.append(r) return np.array(good_peaks)参数说明template_len取200ms覆盖QRS波群。相关系数阈值0.7是经验值低于0.7说明形态差异大可能是异位搏动或噪声。注意模板要减去均值否则基线漂移会影响相关系数。这个方法对房颤效果一般因为房颤的QRS形态其实变化不大但RR间期绝对不齐模板匹配反而会保留所有峰。验证方法在MIT-BIH的208号记录包含室性早搏上跑对比加模板前后的灵敏度和阳性预测率。我实测灵敏度从0.94提到0.97阳性预测率从0.96提到0.98。代价是计算量增加约30%实时场景需要优化。最后说个习惯每次调完参数我都会把滤波前后的波形、积分信号、检测到的R波、以及RR间期序列画在一张图上肉眼过一遍。算法指标再高波形对不上就是有问题。希望帮到你。本文还有配套的精品资源点击获取