简介面向心电图信号处理学习者与科研人员的Python脚本包聚焦原始心电信号的噪声抑制与心跳检测。压缩包共7个文件含核心算法脚本ECG.py、CSV示例原始数据、4张结果可视化图时域、频域、心率及最终结果以及Markdown说明文档其中README还说明了环境依赖与运行方式整包仅211KB。脚本基于numpy与scipy实现基线漂移去除、巴特沃斯低通滤波等预处理流程有效抑制噪声随后通过峰值检测识别QRS波群计算RR间期并估算每分钟心跳数同时生成对比图直观展示处理前后时域与频域差异。代码模块化程度高既有数据读取与预处理步骤也有可视化输出方便研究者针对自身数据修改参数并复用适合医学信号分析与生物医学工程初学者作为入门练习也可为心率变异性等后续研究提供参考。目前已有273人学习浏览整体结构清晰便于快速上手。1. 从noise.csv到心拍序列一个Python脚本能做什么拿到一个名为ECG-Signal-Processing-master的Python工程目录里面有noise.csv、ECG.py和几张结果图时多半意味着你需要对原始ECG信号做噪声抑制、峰值检测和心跳估算。这个脚本要解决的核心问题是手里只有一段未标注的心电数据却要快速得到每个心拍的位置和心率。它不依赖专业心电软件只靠numpy、scipy和matplotlib就能跑通所以非常适合入门信号处理的开发者、刚接触生物医学数据的算法工程师以及需要验证硬件采集质量的嵌入式工程师。接下来我会把这个脚本的常见实现拆开从预处理、R峰定位到心率计算给出可以直接运行的代码片段和参数调整依据让你既能读懂原包也能改成自己的工具。这类脚本的输出通常是一张标注了R峰的result.png一份记录了每搏间隔的heart-rate.PNG以及一个可以直接打印出平均心率的终端结果。对于医生或研究者来说这些数据能帮助他们快速判断信号质量对开发者而言这个流程则是后续心率变异性分析、心律失常分类的基础。在正式开始前你只需要熟悉Python基础、numpy数组操作以及一点信号处理中的频率概念——这三点足够让你把下面每一段代码都落地跑起来。2. ECG信号预处理基线漂移与高频噪声的抑制方法2.1 读取noise.csv并构造时间轴ECG.py要处理的第一件事是读取noise.csv。常见做法是用numpy的loadtxt或pandas的read_csv但我建议先检查一下文件的前三行确认分隔符、表头行数和列顺序。下面这段代码假设第一列是时间秒第二列是电压mV并用np.median计算采样率这样能避免偶发的时间戳跳变把采样率算歪import numpy as np data np.loadtxt(noise.csv, delimiter,, skiprows1) time data[:, 0] ecg_raw data[:, 1] fs 1.0 / np.median(np.diff(time)) print(f采样率: {fs:.2f} Hz)如果你拿到的noise.csv只有一列数据没有时间列那么直接让ecg_raw data[:, 0]并手动指定采样率类似fs 360。很多公开数据集用的是360Hz而可穿戴设备常用250Hz或500Hz。识别采样率很关键因为后面所有滤波器的截止频率、峰值检测的最小间隔参数都依赖它。脚本如果因为读取失败而报错先检查分隔符是否是英文逗号以及skiprows是否正好跳过了表头。最容易被忽略的是时间轴不均匀问题。某些采集设备在蓝牙传输时丢包导致时间戳之间的间隔不是常数。如果你发现np.median和np.mean算出的采样率差异超过1%建议先做插值重采样到均匀网格否则后续滤波器的系数会与实际频率对不上峰值检测的distance参数也会失去意义。2.2 低通滤波与工频陷波用巴特沃斯消除高频噪声原始ECG信号中QRS波群的能量集中在520Hz但肌肉电噪声和电磁干扰会把频谱拉高到几十甚至几百Hz。常见处理是设计一个低通巴特沃斯滤波器截止频率设在4060Hz。这里需要提醒不要为了“干净”把截止频率设在30Hz以下那样会削掉QRS波的高频分量R峰变钝后续定位误差会增大。下面是一个标准的零相位低通滤波实现from scipy.signal import butter, filtfilt def butter_lowpass(cutoff, fs, order4): nyq 0.5 * fs normal_cutoff cutoff / nyq b, a butter(order, normal_cutoff, btypelow, analogFalse) return b, a b, a butter_lowpass(50, fs, order4) ecg_filt filtfilt(b, a, ecg_raw)这里使用filtfilt而不是lfilter是因为filtfilt会正向和反向各滤一遍实现零相位偏移。ECG峰值定位对时间位置很敏感如果使用lfilter信号会产生与频率相关的延迟R峰位置会整体偏移几十毫秒最终计算出的心率虽然可能差别不大但心拍标注会错位。order4是折中值阶数越高过渡带越窄但数值稳定性变差二阶滤波器幅频响应太缓六阶以上容易引入振铃。对于50Hz或60Hz工频干扰可以再加一个带阻滤波器。iirnotch是scipy提供的陷波工具quality参数控制陷波宽度值越大陷波越窄对邻近频段伤害越小但太大容易数值不稳定。下面是处理60Hz干扰的写法from scipy.signal import iirnotch notch_freq 60 # 根据当地电网频率调整中国大部分地区用50Hz quality 30 b_notch, a_notch iirnotch(notch_freq, quality, fs) ecg_notch filtfilt(b_notch, a_notch, ecg_filt)如果你不确定自己的数据里有没有工频干扰先用matplotlib画出频谱找到那个像针一样的尖峰。如果尖峰在50Hz附近就设50在60Hz附近就设60如果两个都有就级联两个陷波器。注意陷波器会稍微扭曲QRS波的高频边缘所以在做心跳检测时不要对信号反复多次滤波。2.3 去基线漂移移动平均与多项式拟合基线漂移主要来自电极移动、呼吸和皮肤阻抗变化频率通常在0.5Hz以下。抑制基线漂移最直接的办法是估计出这条缓慢变化的基线然后从原信号里减掉。这里我推荐使用savgol_filter它本身是平滑滤波器但配一个较大的窗口后可以很好地近似慢变基线。下面的代码用1.5秒窗口和3阶多项式from scipy.signal import savgol_filter window int(fs * 1.5) # 1.5秒窗口覆盖一到两个心搏 baseline savgol_filter(ecg_notch, window_lengthwindow, polyorder3) ecg_clean ecg_notch - baseline这里的逻辑是窗口越长估计出的基线越平滑但也越容易忽略较快的基线起伏polyorder越高拟合越灵活但过高会把QRS波本身的一部分也当成基线。经验上window取12秒polyorder取2或3。如果你发现减完后QRS波幅度明显变小说明窗口太短或polyorder太高把R波包络也吸进基线里了。如果基线仍有缓慢波浪就把窗口加大到3秒。下表列出了不同采样率下常用的滤波参数组合供你对照初始值采样率低通截止频率陷波频率基线窗口polyorder250Hz40~50Hz50/60Hz1.2~2s2~3500Hz45~60Hz50/60Hz1.5~2.5s31000Hz60~80Hz50/60Hz2~3s2~3在实际项目中我一般先用默认参数跑一遍再根据freq-d.PNG和time-d.PNG的对照结果微调。如果频谱中0.5Hz以下仍有明显能量就加大基线窗口如果50Hz附近还有毛刺就降低quality值或增加陷波次数。不要一次性把多个参数同时调大那样很难定位是哪个环节造成的波形异常。3. QRS峰值检测与心跳估算的实现3.1 R峰定位的阈值法预处理完成后中心任务就是找到每个心拍中的R峰。找峰的核心思想是设定一个阈值只保留那些明显高于周围波形的高幅值点。但ECG信号幅值会随呼吸、运动而变化固定阈值很容易漏检或误检。常见做法是先用信号包络动态估计阈值再用find_peaks在原始信号上精确定位。下面这段代码展示了最简洁的实现from scipy.signal import find_peaks abs_ecg np.abs(ecg_clean) window int(0.08 * fs) # 80ms窗口约等于一个QRS波宽度 kernel np.ones(window) / window envelope np.convolve(abs_ecg, kernel, modesame) threshold np.median(envelope) * 2.5 peaks, props find_peaks(ecg_clean, heightthreshold, distanceint(0.25*fs)) print(f检测到 {len(peaks)} 个R峰)这里的关键点有两个一是envelope用80ms滑动平均把QRS波包络提取出来二是阈值取包络中位数的2.5倍。中位数比均值稳健即使有大幅噪声尖峰也不会把阈值拉得过高。distance0.25*fs表示两个R峰之间至少间隔250ms对应最高心率240次/分这在绝大多数场景下足够。如果你发现检测出的峰数量明显多于实际心跳往往是T波被当成了R峰。T波的幅度通常约为R波的三分之一到一半但在运动状态下T波会变大。这时可以把height阈值提高或者增加一个形态约束检查候选峰前后20ms内的斜率是否足够陡峭因为R波比T波陡得多。我的习惯是先用默认参数跑一次把peaks打印到文本文件里对照原始波形人工抽查10秒数据再决定调整方向。3.2 自适应阈值与不应期修正固定阈值在信号平稳时够用但当某段噪声突然增大时阈值可能失效。更稳健的做法是把信号分成若干固定长度窗口在每个窗口内单独计算中位数阈值这样阈值能跟随信号幅值的变化。下面是一个窗口化自适应检测的实现def adaptive_rpeak_detection(ecg, fs, win_sec5.0, factor2.0): win_len int(win_sec * fs) n len(ecg) peaks [] rms_env np.sqrt(np.convolve(ecg**2, np.ones(int(0.03*fs))/int(0.03*fs), modesame)) for start in range(0, n, win_len): end min(start win_len, n) segment rms_env[start:end] if len(segment) 10: continue thresh np.median(segment) * factor local_peaks, _ find_peaks(ecg[start:end], heightthresh, distanceint(0.25*fs)) peaks.extend(start local_peaks) return np.unique(peaks)这里用RMS包络代替绝对值包络RMS对噪声更敏感但能更好地突出QRS波的高能量区域。factor默认为2.0信号质量好时可以用1.8噪声大时调到2.53.0。注意这里虽然用RMS包络算阈值但真正找峰还是在原始预处理信号上这样做不会因为包络的平滑而改变R峰的位置。另一个重要概念是“不应期”。在生理上心肌细胞在一个心搏后有一段时间对刺激不反应体现在心电图上就是QRS波后的一段时间内不会再有另一个R波。distance参数其实就是一个简单的数字不应期。对于严重漏检的数据还可以在检测到峰之后在原始信号局部最大值位置附近细化校正即在候选峰前后±30ms内搜索最大幅值点并替换为精确位置。3.3 从峰值序列计算瞬时心率找到所有R峰位置后心跳间隔RR interval就是相邻峰值索引差除以采样率瞬时心率就是60除以RR间隔。但这中间必须做异常值过滤否则漏检和误检会直接让心率输出变得离谱。常见的过滤规则是以RR间隔的中位数为基准剔除偏离30%以上的点。中位数比均值更抗异常即使有少数几个错误检测也不会把基准拉跑。rr_intervals np.diff(peaks) / fs median_rr np.median(rr_intervals) valid (rr_intervals 0.7 * median_rr) (rr_intervals 1.3 * median_rr) rr_valid rr_intervals[valid] heart_rate 60.0 / rr_valid print(f平均心率: {heart_rate.mean():.1f} bpm)下面这张表总结了不同场景下的异常过滤阈值场景下限上限说明静息状态0.8 * median1.2 * median心率平稳范围可收紧日常活动0.7 * median1.3 * median允许呼吸性波动运动状态0.6 * median1.4 * median心率变化快防止误删过滤完之后如果你需要每个心拍的时间戳可以这样输出r_peak_times peaks / fs r_peak_times_valid r_peak_times[valid]这里的valid数组长度比rr_intervals短一个因为rr_intervals是通过np.diff得到的。如果你把valid直接索引到r_peak_times上长度对不上。所以正确写法是先用valid过滤rr_intervals再决定保留哪些峰。常见错误是直接用错误长度的布尔数组索引结果导致越界或错位。另外脚本里通常会把R峰位置标记在信号图上生成类似heart-rate.PNG的图片这时需要把时间轴和检测到的peak索引一起传入绘图函数。4. 参数调优与常见坑让脚本在不同设备数据上可用4.1 滤波器阶数与截止频率的选择很多新手拿到脚本后直接在自带数据上跑通但换了自己的采集设备就一团糟。问题往往出在滤波器参数与采样率不匹配。巴特沃斯低通滤波器的阶数越高下降沿越陡但阶数过高会导致filtfilt出现振铃尤其在QRS波这样陡峭的波峰附近振铃会被误检为额外的心搏。我的选择标准是采样率低于250Hz用2阶250500Hz用4阶超过1000Hz才考虑6阶。如果你的信号来自胸贴式单导联设备采样率通常在125250Hz2阶滤波配合40Hz截止频率足够如果是医院12导联设备采样率在500Hz以上可以设到60Hz甚至80Hz保留更多心电细节。观察脚本自带的freq-d.PNG能帮你判断截止频率是否合适。频域图中QRS主能量应该完整保留而高于截止频率的部分应明显衰减。如果你发现滤波后的波形在每一个QRS波前后出现小幅振荡那就是阶数过高或截止频率过低导致的振铃。这时候应该降低阶数或提高截止频率而不是继续加强滤波。4.2 阈值和窗口的调节策略find_peaks中的height和distance是一对需要联调的参数。height过高会漏检小的QRS过低则会把T波甚至把滤波后的残余噪声尖峰当作R峰。一个实用的方法是先画出信号的直方图观察幅值分布。正常的ECG信号幅值分布中QRS波会拖出一个长尾T波和噪声则集中在低幅值区。阈值设在高斯分布尾部与QRS长尾交界处附近通常在中位数的23倍之间。如果检测结果中出现“心率是实际两倍”的现象基本可以判定是T波被误检。T波在II导联中与R波方向一致在运动状态下幅值会升高有时能达到R波的70%。此时除了调高height还可以利用生理知识R波与T波之间的间隔通常不小于200ms。增大distance到0.3秒可以强制忽略T波。反过来如果检测结果出现“心率是实际一半”说明有R波被漏检常见原因是某个心拍的幅值突然变小比如呼吸性基线波动未完全去除。这时可以降低factor到1.51.8或缩短RMS包络窗口到20ms让阈值更快跟上幅值变化。另一种提升稳健性的做法是模板匹配。从信号中裁剪几个清晰的QRS波作为模板然后通过滑窗计算相关系数。相关系数高于0.8的位置视为R峰。这种方法对小幅度QRS更鲁棒但计算量大而且模板需要人工挑选不适合全自动流水线。我一般只在离线分析中使用在线场景还是用自适应阈值更实用。4.3 处理noise.csv格式不一致的常见错误自带的noise.csv是标准逗号分隔文件但你自己的数据可能完全不同。常见情况包括分隔符是分号、第一列不是时间而是采样序号、电压单位是μV而脚本预期是mV、文件带BOM头导致第一列名带特殊字符。读取时建议使用pandas它处理混乱格式的能力明显更强import pandas as pd df pd.read_csv(noise.csv, delimiter,, enginepython, encodingutf-8-sig) t df.iloc[:, 0].values v df.iloc[:, 1].values.astype(float)如果时间轴不均匀比如蓝牙丢包导致某些时间戳缺失你需要做插值。最简单的做法是np.interp把原始时间戳插值到[0, duration]上均匀网格target_fs 250 target_t np.arange(0, t[-1], 1/target_fs) v_uniform np.interp(target_t, t, v) fs target_fs注意np.interp要求原始t单调递增如果数据中有重复时间戳要先做去重。还有一种情况是电压幅值漂移特别大比如某一段信号整体抬高了2mV这其实属于基线漂移应该回到第2.3节的去基线步骤重新处理而不是在峰值检测阶段强行用阈值去适应。下面是一张快速排查表能帮你定位大多数检测异常现象可能原因调整方向心率约为实际一半阈值过高或distance过大降低height减小distance心率约为实际两倍T波误检提高height增大distance到0.3sR波变钝峰位偏移低通截止频率太低提高到60Hz以上低频漂移仍然明显基线窗口太短增大到2~3s某段信号全部漏检局部噪声过大改用滑动窗口自适应阈值5. 从result.png到heart-rate.PNG验证你的处理链路5.1 对比时域和频域图脚本包里的time-d.PNG展示的是预处理前后信号叠加freq-d.PNG是频谱对比。拿到这两张图第一步不是数峰而是看频谱。正常结果中50Hz和60Hz附近的尖峰应该消失1Hz以下的低频能量被明显压制但1020Hz区间的QRS主峰应该保留。如果主峰没了那滤波过度了。第二步看时域QRS波应当尖锐T波没有被放大基线在0mV附近波动波动幅度不超过0.2mV。5.2 人工核对R峰标注在信号上叠加标注点用matplotlib逐屏检查。下面这段代码把检测到的R峰画成红色圆点并列出每个峰的时间import matplotlib.pyplot as plt r_times peaks / fs plt.figure(figsize(12, 4)) plt.plot(time, ecg_clean, b-, alpha0.7, labelECG) plt.plot(r_times, ecg_clean[peaks], ro, markersize3, labelR peaks) plt.xlabel(time (s)) plt.ylabel(amplitude (mV)) plt.legend() plt.show()核对时注意两点一是红点是否都在R波顶点附近偏差超过20ms说明滤波产生了相位偏移或峰值查找窗口设置不当二是相邻红点的间隔是否规律若某个间隔突然变成前后间隔的两倍说明中间漏检了一个心拍若出现极短间隔则多半是误检。5.3 用滑动窗口输出实时心率如果你想把脚本扩展成流式计算可以每次读入2秒数据检测峰值后只记录当前窗口最后一个峰值的位置与下一窗口拼接。拼接时注意去重如果后一个窗口检测到的第一个峰与前一窗口最后一个峰的时间差小于0.3秒则认为是同一个峰。输出心率时使用指数平滑能让数值更稳定hr_smooth alpha * hr_curr (1-alpha) * hr_smoothalpha取0.3左右。这个技巧对可穿戴设备的实时显示很实用因为瞬时心率容易受单次误检影响而平滑后的数值既保留了变化趋势又不会剧烈跳动。本文还有配套的精品资源点击获取