企业官网建设流程全解析

简介：这是一套面向生物医学信号处理开发者与算法学习者的Python实现心电图R峰检测算法源码包，重点针对MIT-BIH等心电数据库的心跳定位、HRV计算与心律失常指标评估，适合有Python基础并希望深入ECG分析的读者。包内共31个文件，约600KB，包含10个Python脚本和7个C程序：Python侧覆盖Pan-Tompkins、Jafarnejadsani等经典检测器及测试、绘图流程，C程序负责滤波与心率统计，另有CSV/TSV结果文件和不同采样率模板，可直接复现基准测试。整个工程自带安装脚本、Makefile、示例数据和运行入口，从环境配置到结果输出层次清晰，能帮助读者快速跑通检测流程，并对照源码理解阈值自适应、峰谷判定与噪声抑制等细节。已有2007人学习，适合正在开展可穿戴监测、远程心电预警或生物信号课程设计的开发者作为参考。

1. 拿到心电图数据之后，第一件事为什么是找 R 波

手里有一份心电图数据，几万甚至几十万个采样点，真正能用的信息其实只有一小段：每一次心跳对应一个 QRS 波群，其中最陡、最高、最醒目的那个尖峰叫 R 波。所谓 Python 心电图心跳检测算法，核心就是把这个尖峰的位置找出来，剩下的心率计算、心律不齐分析、疲劳监测，全都建立在 R 波定位之上。很多人一开始就想着上深度学习，但我个人建议先把经典算法跑通，因为这类算法不依赖标注数据也能调试，参数透明，出了问题一眼能看出是滤波还是阈值的问题，非常适合低成本落地和做技术预研。

这篇文章要讲的是我从零搭起一组用 Python 实现的心跳检测流程：先做带通滤波把噪声压下去，再用斜率、幅度和滑动窗口积分把 R 波从波形里抠出来，最后用公开数据库核对精度，并整理出我在实际数据上踩过的几个坑。如果你刚装好 Python 和 NumPy 还没想好从哪个方向入手，或者已经用阈值法做过检测但总觉得误检很多，这篇应该能给你一条可以直接照着走下去的路。

2. 预处理做对一半：为什么先带通滤波再检测，而不是直接找峰值

2.1 心电信号里的噪声都藏在哪些频段

直接对原始心电信号做峰值检测，结果通常惨不忍睹。原因很简单：ECG 是一种低频信号，但采集过程中会叠加上各种噪声。基线漂移来自电极移动和呼吸，频率在 0.5 Hz 以下；肌电干扰来自肌肉紧张，频率在 20 Hz 以上；工频干扰在 50 或 60 Hz，幅度可能比 QRS 波本身还大。而 QRS 波群的主要能量集中在 5 到 20 Hz 之间。这一条频带信息就是整个预处理设计的出发点。

常见做法是先设计一个带通滤波器，只保留 5 到 20 Hz 这个频段，把低频漂移和高频肌电一起切掉。我用的是 Butterworth 滤波器，因为它的通带内最平坦，参数少且行为可预期。很多人第一次写滤波代码时直接调用 scipy.signal.butter 拿到系数就上 lfilter，结果发现波形整体往后平移了，R 波位置全部偏移几十毫秒。这是因为普通滤波器都有相位延迟，而心跳检测对位置敏感，所以我这里用 filtfilt 做零相位滤波，前向和反向各过一次，相位延迟相互抵消。

2.2 用 scipy 搭一版零相位带通滤波

下面这段代码是我常用的预处理函数，输入原始信号和采样率，输出滤波后的信号。

import numpy as np from scipy.signal import butter, filtfilt def bandpass_ecg(signal, fs, lowcut=5.0, highcut=20.0, order=2): """ 对心电信号做零相位带通滤波。 signal: 一维数组，原始 ECG 波形 fs: 采样率，单位 Hz lowcut/highcut: 通带上下限 order: 滤波器阶数，2 阶通常够用 """ nyquist = 0.5 * fs low = lowcut / nyquist high = highcut / nyquist # 注意 butter 返回的是传递函数的分子分母系数 b, a = butter(order, [low, high], btype='band') # filtfilt 做零相位滤波，R 波位置不会漂移 filtered = filtfilt(b, a, signal) return filtered

逻辑说明：butter 函数把频率参数归一化到 Nyquist 频率之后，返回 IIR 滤波器的分子分母系数。filtfilt 先正向滤波一次，再把结果反转后反向滤波一次，抵消相位偏移，代价是计算量翻倍，但对于离线处理的数据完全可接受。order 取 2 是因为更高阶的滤波器虽然过渡带更陡，但数值稳定性和波形畸变风险也会增加，对心跳检测来说 2 到 4 阶足够。

参数说明：lowcut 建议设在 3 到 5 Hz，太高会把 QRS 波的低频成分切掉导致波形变尖，太低则挡不住基线漂移。highcut 建议设在 15 到 25 Hz，T 波和 P 波频率较低，设到 20 Hz 时它们会被明显抑制，对 R 波影响不大。fs 必须跟采集设备一致，200 Hz 和 1000 Hz 的归一化计算完全不同，这一步错了后面全是错的。

2.3 滤波前后到底该看什么

怎么确认滤波参数是对的？我的习惯是把原始信号和滤波后信号画在同一张图上对比，横坐标按秒而不是按采样点。画图的时候注意一个坑：采样点太多横坐标会挤成一团，根本看不出节奏，所以要自己换算成秒再画。

判断标准有三个。第一，滤波后的基线应该基本是一条平线，没有大幅上下起伏；第二，QRS 波群仍然保留明显的尖峰形态，没有被磨圆；第三，T 波幅度应该明显小于 R 波，如果 T 波还跟 R 波一样高，说明高通截止频率设得太低了。

3. 核心检测逻辑：从滤波信号里把 R 波的位置抠出来

3.1 斜率、平方、滑动窗口积分：把尖峰变成一块突出的区域

滤波之后的信号仍然是一个连续波形，直接在原始幅度上做阈值检测会遇到一个问题：QRS 波的幅度在不同导联、不同人身上差异很大，有的 R 波 0.5 mV，有的高达 2 mV，一个固定阈值没法通用。经典 Pan-Tompkins 算法给出了一套不依赖绝对幅度的处理流程，核心分三步：求导、平方、滑动窗口积分。

求导是为了突出斜率，R 波的上升沿是心电信号中变化最快的地方，求导之后在对应位置会产生一个明显的正峰值，而 T 波和 P 波因为变化缓慢，导数很小。平方操作把所有负值变成正值，同时进一步放大高频成分的差异。最后做滑动窗口积分，把一段短时间内的能量累加起来，QRS 波群所在的位置会形成一个凸起的平台，这个平台远比原始尖峰容易检测。

下面这段代码把三步连在一起实现，注意移动窗口积分我用了累计和做差分，避免反复求和导致计算量爆炸。

def qrs_energy_transform(filtered, fs, win_sec=0.15): """ 对滤波后的信号做差分、平方、滑动窗口积分。 返回变换后的能量曲线和窗口长度。 """ # 一阶差分近似求导，注意长度减少 1 diff = np.diff(filtered) # 平方突出高频成分，转为正值 squared = diff ** 2 # 窗口长度按采样率换算，150ms 是常见值 win_len = max(1, int(win_sec * fs)) # 用 cumsum 高效算滑动窗口和 cum = np.cumsum(np.concatenate(([0.0], squared))) integrated = np.zeros_like(squared) for i in range(len(squared) - win_len + 1): integrated[i] = cum[i + win_len] - cum[i] return integrated, win_len

逻辑说明：np.diff 得到的导数序列长度比原信号少 1，所以后面所有操作的索引都会整体偏移一点，定位 R 波时要补偿回来。cumsum 的用法是关键，它先算累计和，然后窗口和用两次索引相减得到，这里看起来是 for 循环，但实际只访问 O(n) 次，比每次 sum 一个窗口快非常多。

参数说明：win_sec 是滑动窗口的物理时长，我这里取 150 ms。这个值的物理意义是 QRS 波群本身的持续时间约 80 到 120 ms，窗口太小会把一个 QRS 波切成两段产生双峰，窗口太大又会让两个相邻的心跳融合在一起，150 ms 是经验上最稳的选择。如果采样率是 250 Hz，win_len 就是 38 个采样点；如果是 1000 Hz，win_len 就是 150 个采样点，务必按 fs 换算而不是用固定点数。

3.2 动态阈值和不应期：两个最容易翻车的设计点

能量变换之后，R 波位置会变成一处处明显的凸起，但这还不够，你仍然需要回答一个问题：多高的凸起算一次心跳？如果设固定阈值，信号幅值漂移时就会失效。我常用自适应动态阈值，思路是维护两个量：信号峰值水平 signal_peak 和噪声水平 noise_peak。每次检测到新峰值，就按比例更新这两个水平，阈值取两者的加权值。

此外还有一个必须处理的问题：R 波本身在积分后有一定宽度，如果不做抑制，同一个心跳会被连报好几次。解决办法是不应期：检测到一次 R 波后，在其后 200 ms 内不再接受新的检测，同时对这个时间段内的最大峰值做一次比较，保留最高峰作为真实位置。这就是把精度从“次”提升到“位置”的关键一步。

def detect_r_peaks(integrated, fs, win_len, threshold_ratio=0.4, refactory_sec=0.2): """ 在能量曲线上做动态阈值检测，返回 R 波位置索引。 integrated: 上一步滑动窗口积分的输出 fs: 采样率 win_len: 滑动窗口长度，用于补偿差分导致的位置偏移 """ peaks = [] signal_peak = 0.0 noise_peak = 0.0 refractory_samples = int(refactory_sec * fs) last_peak_idx = -refractory_samples # 初始化为足够远 for i in range(len(integrated)): val = integrated[i] # 更新噪声/信号峰值估计 if val > signal_peak: signal_peak = 0.8 * signal_peak + 0.2 * val else: noise_peak = 0.8 * noise_peak + 0.2 * val threshold = noise_peak + threshold_ratio * (signal_peak - noise_peak) if val > threshold and (i - last_peak_idx) > refractory_samples: # 进入不应期：在这个时间窗内找局部最大值 search_start = i search_end = min(len(integrated), i + refractory_samples) local_max = search_start + int(np.argmax(integrated[search_start:search_end])) peaks.append(local_max) last_peak_idx = local_max # 把搜索范围跳过去，直接以局部峰为下一次起点 i = local_max # 差分导致波形整体往前偏移，补偿回来得到原始信号上的位置 peaks = [p + win_len // 2 + 1 for p in peaks] return peaks

逻辑说明：signal_peak 和 noise_peak 都用了指数滑动平均，0.8 和 0.2 的权重决定了跟踪速度。权重越大跟踪越慢，抗干扰能力越强，但对幅值突变反应迟钝。阈值本身是噪声水平加上信号与噪声差值的 40%，这个比例意味信号峰只要比噪声高出一定比例就会被接受，不会因为整体幅度变大变小而失效。

参数说明：threshold_ratio 的常见取值在 0.3 到 0.5 之间，取 0.4 是比较均衡的起点。refactory_sec 设 200 ms 对应正常人最大心率约 300 bpm，超过 300 的检测结果应当视为噪声。如果你处理的是运动员静息数据，200 ms 完全没有问题；如果是运动负荷测试数据，心率可能到 180 bpm，对应 R-R 间期约 333 ms，200 ms 依然安全。

3.3 这一版算法的默认参数一览

4. 检测完怎么算对：用公开数据库核对检测精度

4.1 有标注数据的评估：MIT-BIH 数据库的读取思路

算法写完不能只看两张图感觉差不多就行，要用带标注的数据算准确率。最常见的公开数据集是 MIT-BIH Arrhythmia Database，每个记录包含两路心电信号和一组由专家标注的心跳位置。标注信息按 beat 类型分类，N 代表正常心跳，V 代表室性早搏，但心跳检测这一步不用区分类型，只需要把标注里的绝对位置读出来即可。

这个数据库的原始格式不是常见的那种 CSV，而是 .dat、.hea、.atr 三个文件配合使用，需要从 PhysioNet 提供的工具或现成解析库读取。不过你现在只需要明确一点：不管格式怎么解析，最终你手里应该有两组数字，一组是算法报告的 R 波位置索引，另一组是专家标注的 R 波位置索引，接下来做的就是两组数字的匹配。

4.2 写一个评估函数：TP、FP、FN 怎么算

匹配规则我通常用容差窗口：算法检测位置与某个标注位置之间距离小于 150 ms，就认为是一次真阳性 TP。一个标注位置只能匹配一次，一个检测位置也只能匹配一次，这样算出的指标才严格。

def evaluate_detection(detections, annotations, fs, tolerance_sec=0.15): """ 计算检测结果的 TP、FP、FN 和灵敏度、阳性预测值。 detections: 算法检测到的 R 波位置索引列表 annotations: 专家标注的 R 波位置索引列表 fs: 采样率 """ tolerance = int(tolerance_sec * fs) detections = sorted(detections) annotations = sorted(annotations) tp = 0 matched_det = set() matched_ann = set() for d_idx, d in enumerate(detections): best_ann = None best_dist = float('inf') for a_idx, a in enumerate(annotations): if a_idx in matched_ann: continue dist = abs(d - a) if dist <= tolerance and dist < best_dist: best_dist = dist best_ann = a_idx if best_ann is not None: tp += 1 matched_det.add(d_idx) matched_ann.add(best_ann) fp = len(detections) - tp fn = len(annotations) - tp sensitivity = tp / len(annotations) if annotations else 0.0 ppv = tp / len(detections) if detections else 0.0 return {"TP": tp, "FP": fp, "FN": fn, "Sensitivity": sensitivity, "PPV": ppv}

逻辑说明：这里做了双重去重，一个标注只能被一个检测匹配，一个检测也只能匹配一个标注。最朴素的双重循环在数据量小时没有问题，一条记录通常只有几百次心跳，所以不需要优化。容差 150 ms 是一个在临床上被普遍接受的误差范围，它既没有严格到把一点点相位偏移都算错，也没有宽松到把误差很大的检测也放过去。

参数说明：tolerance_sec 我习惯取 0.15，但如果你做的是心率变异性分析，对 R 波位置的精确度要求更高，建议收紧到 0.05 到 0.08 秒。另外，如果你用的滤波器阶数较高或波形形态特殊，检测位置可能系统性地偏移几十毫秒，这种情况先检查滤波相位，而不是放大容差掩盖问题。

4.3 没有标注的数据怎么验证

现实项目里大部分数据没有专家标注，这时候可以换一种思路：不验证每个 R 波位置，而是验证 R-R 间期序列的合理性。正常心跳的 R-R 间期应该在 600 到 1500 ms 之间，且相邻间期变化不会超过 30% 左右。如果算法测出来的 R-R 间期出现大量小于 300 ms 的值，说明检测到了 T 波；如果间期频繁大于 2000 ms，说明漏检了心跳。

def check_rr_intervals(r_peaks, fs, min_rr=0.3, max_rr=2.0): """ 用 R-R 间期做合理性检查，返回可疑片段。 r_peaks: R 波位置索引列表 fs: 采样率 """ rr = np.diff(r_peaks) / fs # 单位换算为秒 problems = [] for i, interval in enumerate(rr): if interval < min_rr or interval > max_rr: problems.append((i, interval)) # 统计异常占比 anomaly_ratio = len(problems) / len(rr) if len(rr) else 0.0 return problems, anomaly_ratio

逻辑说明：这段代码把相邻 R 波位置的差值从采样点换算成秒，然后检查是否落在合理范围内。异常占比如果超过 5%，通常说明检测参数不适合这条数据，而不仅仅是个别噪声。这个方法不能帮你确认检测到的位置百分百正确，但能快速把明显不合格的结果筛选掉，是拿到新数据后第一个应该跑的脚本。

5. 用 Python 实现心跳检测的 5 个避坑记录：现象、原因、解决

5.1 基线漂移导致阈值整体失效

现象：滤波后的信号整体仍然像一条波浪线，R 波在波峰处被检测到，在波谷处被漏掉，检测结果忽多忽少。原因：高通截止频率设太低，假设设成 0.5 Hz，其实呼吸引起的基线漂移仍在通带内，滤波没有把它们抑制干净。解决：把 lowcut 从 0.5 Hz 提到 5 Hz。基线漂移频率集中在 0.5 Hz 以下，5 Hz 的高通虽然会损失一部分 ST 段信息，但对 R 波检测毫无影响。这个坑最容易在采样率 500 Hz 以上的数据里出现，因为低频成分在长记录中特别明显。

5.2 T 波被当成 R 波，心率直接翻倍

现象：检测到的心率大约是真实值的两倍，画出的检测点落在波形的小鼓包上。原因：T 波在滤波后仍然有一定幅度，如果信号能量曲线中 T 波的峰值超过了阈值，且不在不应期内，就会被当作一次心跳。T 波离 R 波通常有 300 ms 左右的间隔，恰好落在 200 ms 不应期之外。解决：优先检查 highcut，如果 highcut 是 20 Hz，T 波幅度明显、误检频繁，就把 highcut 降到 12 Hz；如果还有问题，再把 threshold_ratio 从 0.4 提高到 0.5。T 波比 R 波变化缓慢，斜率特征上天然有区分度，这也是为什么先做差分而不是直接用原始波形来找峰。

5.3 滤波后 R 波位置整体偏移，精度怎么都提不上去

现象：评估结果中 TP 很多但位置平均偏差超过 30 ms，拿去算心率没问题，算心率变异性就超差。原因：用了 lfilter 而不是 filtfilt。lfilter 是因果滤波，输出必然滞后于输入，滞后量跟滤波器阶数相关，阶数越高延迟越大。这个偏差在单个波形上肉眼看不出来，但和标注一对比立刻暴露。解决：换成 filtfilt 做零相位滤波，代码只有一处变化，但精度提升是质的。如果你用的是实时流式数据，filtfilt 因为需要未来数据不可用，就只能接受延迟并在时间戳上做补偿。

5.4 同一个心跳被报了两次，双峰问题

现象：算法检测出的 R 波之间有大量 80 到 100 ms 的极短间隔，画出来发现同一个 QRS 波上升沿和下降沿各触发了一次。原因：差分平方之后，QRS 波群在上升沿和下降沿都产生了明显的能量峰，如果滑动窗口时间小于 QRS 本身宽度，这两个峰就不会融合成一个。解决：把 win_sec 从 0.10 提到 0.15 秒。如果数据采样率特别高，比如 1000 Hz，还要确认 win_len 至少覆盖 100 ms 以上。另外，不应期的存在能挡住第二峰，但不应期会把漏在两个峰之间的真实心跳也挡掉，所以根本解法还是窗口宽度正确。

5.5 采样率假设写死，换数据就翻车

现象：同一套代码在一条数据上跑得很好，换一条数据后异常率飙升，检查发现采样率从 360 Hz 换成了 250 Hz，而滤波窗口和滑动积分窗口全部用固定点数。原因：代码里写死了窗口长度，没有除以实际的 fs。200 ms 不应期在 360 Hz 下是 72 个点，在 250 Hz 下是 50 个点，差 22 个点足以导致完全不同的检测行为。解决：所有时间相关参数全部用 fs 换算，不要在函数里写 magic number。这是我见过发生频率最高的问题，尤其是从网上复制代码时，别人的数据采样率跟你不一定相同。

6. 进阶验证技巧：用合成信号自检，再批量扫参数

最有效的自检工具是先构造一段标准的仿真心电信号，比如用几组高斯函数模拟 P 波、QRS 波和 T 波，再叠加噪声和基线漂移，然后用你的检测算法去测。合成信号的好处是每个 R 波位置是已知的，算法表现可以被精确量化，这比在真实数据上反复调整阈值效率高得多。要注意合成信号只是验证流程正确性，不能替代真实数据评估，因为合成波形往往比真实数据干净。

批量调参是另一个值得投入的小工具。把 threshold_ratio 从 0.2 到 0.6 按步长 0.05 扫描，把 win_sec 从 0.08 到 0.20 扫描，每次跑完整条记录，记录灵敏度和阳性预测值，最后选两个指标同时最高的组合。这个扫描在单条记录上花费的时间只有几秒到几十秒，换来的却是参数选择的确定性。我自己做每个新项目时都会把这条数据跑一遍，记录参数和结果，而不是靠感觉和运气去调，因为真实数据噪声形态千差万别，你没法凭空猜出最合适的阈值。最后再补一句：任何检测算法跑新数据集都要先画一次波形，亲眼确认检测点落在心拍上，这个过程不能省，希望帮到你。

本文还有配套的精品资源，点击获取