简介Matlab心电信号峰值检测实践资源面向本科、硕士及教研人群属于Matlab基础算法与应用结合的入门级素材。资源以心电图信号处理为切入点展示如何利用Matlab完成信号读取、波形展示与峰值定位帮助读者快速建立信号处理与算法实现的基本框架适合零基础或刚接触心电分析的学习者边看边调试。压缩包共7个文件整体仅126KB包含一个.m格式主程序脚本、五张.jpg运行结果截图与一份.txt辅助说明脚本可直接运行截图呈现峰值检测输出说明则补充数据或操作要点方便对照学习。从内容预览来看程序文件以心电图峰值检测为核心覆盖从原始数据到峰值提取的完整演示流程运行结果图展示不同环节波形其中的说明文本可核对中间数据或参数配置适合理解Matlab基础语法与信号处理思路。目前已有120人学习浏览整体轻量紧凑可作为Matlab信号处理课程的小型练习或教研参考。1. 心电信号心电图峰值检测先别急着用 findpeaks拿到一份心电信号第一反应往往是打开 Matlab 敲一行findpeaks然后就等着 R 波一个接一个蹦出来。真实情况是基线漂移让波峰上下乱晃工频干扰在信号上叠加了一堆毛刺T 波在某些导联上比 R 波还高findpeaks要么把 T 波当成 QRS 峰值要么因为阈值设不准漏掉真值。做心电图峰值检测核心不是找一个“找峰函数”而是先把 QRS 波群的形态学特征搞清楚再决定滤波、阈值和不应期策略。我拆的这份 Matlab 心电信号峰值检测源码走的正是经典 Pan-Tompkins 路线适合正在做生物医学信号处理课程设计、医疗算法入门或者需要从裸信号里提取心跳位置的从业者。先去下下来再看下面的原理和代码会顺很多。2. 心电峰值检测为什么难基线漂移、工频干扰与滤波选型2.1 为什么 R 波不是单纯“局部极大值”心电信号里QRS 波群是最明显的特征R 波峰值通常也是整个心跳周期里幅值最大的点。但如果直接对原始信号做局部极大值检测你会发现三个问题第一呼吸引起的基线漂移让整个信号像坐在一艘晃动的船上幅值每几秒就缓慢起伏一次R 波的真实高度是相对基线而言的不是相对零线第二50Hz 或 60Hz 工频干扰叠加在信号上可能把每个 R 波顶部的平滑峰变成锯齿状产生多个伪局部极值第三在某些导联或某些病理信号里T 波幅值可能接近甚至超过 R 波单看高度根本没法区分。这就是为什么心电图峰值检测不能只依赖findpeaks。findpeaks的MinPeakHeight和MinPeakDistance需要你事先知道信号的大致幅值和心率范围但真实采集的信号幅值会随电极接触质量变化心率也会随运动或病理波动。阈值写死了换个受试者就翻车。Pan-Tompkins 这类经典算法的思路是用带通滤波把噪声和漂移压下去再用微分、平方、滑动积分三个操作把 QRS 的特征放大最后用自适应阈值和不应期规则去锁定 R 波位置。它不假设一个固定幅值而是让阈值跟随信号能量自动调整。我一般会先跑一遍算法看输出再回头调滤波参数因为大多数人第一步踩的坑就出在滤波上要么滤得太干净把 R 波边缘也削平了要么滤完发现基线还在。处理心电信号滤波的目标不是“光滑”而是“保留 QRS 波群的主要能量同时压掉基线漂移和肌电噪声”。提示心电信号的 QRS 主频大约在 520HzP 波和 T 波更低而基线漂移通常低于 0.5Hz工频干扰在 50/60Hz。带通范围选 515Hz 可以看到清晰的 QRS 轮廓。2.2 滤波器选型巴特沃斯还是椭圆IIR 还是 FIR做带通滤波Matlab 里最常用的就是butter搭配filtfilt其次是ellip再就是designfilt设计 FIR 等波纹滤波器。三者我都实际跑过说下感受。巴特沃斯的通带最平坦相位响应在filtfilt零相位处理下不会造成波形偏移适合做 QRS 波群检测的预处理。缺点是过渡带比较宽阶数要提到 46 阶才能把 0.5Hz 以下的基线漂移压得比较干净而阶数越高计算量越大但一段几万点的信号在 Matlab 里也就是毫秒级的事完全不用心疼性能。椭圆滤波器过渡带更窄同样的阶数能实现更陡的衰减但通带有纹波在某些信噪比场景下会让 R 波峰值略微抖动。FIR 等波纹滤波器相位线性可控性好但同样的过渡带要求阶数轻松上百对长信号来说filter的计算量会明显上升。我自己在课程设计里偏好butter(6, [5 15]/(fs/2), bandpass)再filtfilt。原因有三零相位不会让 R 波位置平移巴特沃斯在通带内增益平坦不会改变 QRS 各波形的幅值比例实现只要两行代码参数好解释写实验报告也方便。如果检测率上不去再换椭圆或 FIR 对比而不是一上来就把滤波器搞得特别复杂。还有一个细节必须注意butter的截止频率参数是归一化到奈奎斯特频率的也就是Wn 截止频率 / (fs/2)。如果你的采样率是 500Hz那 15Hz 对应的 Wn 是15/250 0.06。这个换算错了滤波器可能直接把整个信号都滤没了。fs 500; % 采样率常见 ECG 设备为 250/360/500 Hz f_low 5; % 高通截止频率压掉基线漂移 f_high 15; % 低通截止频率保留 QRS 主能量 [b, a] butter(6, [f_low f_high]/(fs/2), bandpass); ecg_filtered filtfilt(b, a, ecg_raw);这段代码里的阶数 6 是巴特沃斯带通滤波器的总阶数实际等效于 3 阶高通加 3 阶低通。filtfilt做零相位滤波会先正向滤一遍再反向滤一遍抵消相位偏移这是权力说明里最关键的点——R 波的位置不会被滤波拉偏代价是信号首尾各有一段被边界效应影响。后面避坑章节我会专门讲这段怎么处理。参数调整上如果基线漂移特别严重f_low 可以提到 68Hz如果肌电噪声明显f_high 可以降到 12Hz但降太多会把 QRS 波群本身的高频成分削掉检测率反而下降。我建议以 15Hz 为起点逐 Hz 往下试每次看漏检数不凭感觉定。3. 核心算法拆解Pan-Tompkins 的微分、平方与自适应阈值3.1 微分与平方让 R 波从背景里“跳”出来带通滤波之后QRS 波群形态被保留基线漂移基本消失但 T 波仍然存在而且某些信号里 T 波高度接近 R 波。这时候单纯看幅值不够得看“变化率”。QRS 波群的斜率远大于 T 波和 P 波R 波的下降沿和上升沿都特别陡这是它最显著的特征。Pan-Tompkins 算法的第一步是微分用一阶差分近似求导。对离散信号来说就是y[n] x[n] - x[n-1]。这一操作把 QRS 波群的高斜率成分放大成明显的正负脉冲而 T 波这种缓变的波形成分被大幅衰减。第二步是逐点平方y[n] y[n]^2。平方有两个作用一是把微分的正负脉冲全部变成正的方便后续积分累加二是进一步放大 QRS 波群相对于 T 波、P 波的幅值差距因为幅值大的点在平方后优势更明显平方是“强者愈强”。Matlab 里这两步可以这样写diff_h diff(ecg_filtered); % 一阶差分近似求导 squared diff_h .^ 2; % 逐点平方放大高频特征注意diff会让信号长度减 1后面做滑动积分或对齐索引时要记得补回这个偏移。你可以在diff前用[ecg_filtered(1); diff(ecg_filtered(:))]的方式保持长度不变我习惯在差分后squared(end1) squared(end)补一位确保索引对得上。这些预处理流程从滤波到微分、平方再到移动窗口积分是 Pan-Tompkins 的固定公式套路。不过要强调一点Matlab 里的diff是纯按样本间差分计算的。如果你的采样率不是 500Hz 而是 360Hz微分结果的幅度会整体按比例变化但由于后面要做的是相对阈值判断绝对幅度变化不影响最终检测结果这也是这套方法鲁棒的一个原因。3.2 滑动窗口积分与自适应阈值锁定峰值位置微分和平方之后信号变成了一个个窄尖峰每个 QRS 对应一个但尖峰可能有多个毛刺峰直接找极大值还是会产生误检。Pan-Tompkins 的做法是再做一次滑动窗口积分moving-window integration。积分的作用是把单个窄尖峰“抹宽”成一个平顶波多个邻近毛刺被融合成一个整体这样每个 QRS 波群在这个积分信号上只留下一个清晰的平台平台前沿的位置就有规律可循了。滑动窗口积分的公式是y[n] (1/N) * sum(x[n-N1:n])其中 N 是窗口长度。窗口长度取多少很关键经验值是取 150ms 左右的信号长度也就是 N round(0.15 * fs)。窗口太短融合不了毛刺T 波还可能产生伪平台窗口太长会把 QRS 和紧随其后的 T 波糊在一起反而造成漏检或误检。我一般从 0.15 起步看到积分波形上每个 QRS 对应一个完整平顶就满意。win_width round(0.15 * fs); % 窗口宽度约 150ms integral_signal conv(squared, ones(1, win_width)/win_width, same);用conv实现滑动求和是最快的比自己写 for 循环好得多。ones(1, win_width)/win_width就是长度为 win_width 的矩形窗same让输出和输入长度一致。这块的边界效应同样存在后面会讲。自适应阈值才是 Pan-Tompkins 的灵魂。最朴素的做法是把信号的峰值或平均值算出来乘以一个系数作为初始阈值然后每检测到一个 QRS 波群就用这个波群的峰值去更新当前阈值。常见的更新策略是threshold alpha * peak_signal; % alpha 常取 0.30.6 threshold alpha * estimated_peak (1-alpha) * old_threshold; % 递推更新第二个式子里的 alpha 是平滑系数决定新检测到的 R 波峰值对阈值的贡献权重。取值越大阈值反应越灵敏适合心率剧烈变化的场景取值太小阈值变化缓慢适合心率平稳的静息心电。实际调参时我一般把 alpha 放 0.2 到 0.5 之间试。自适应阈值同时还配了一个“不应期”规则检测到一个 R 波之后的 200ms约 0.2 秒内不再接受新的检测结果。因为心跳的生理极限决定了两次 QRS 之间必须有足够的时间间隔200ms 对应 300 次/分的极限心率正常场景下不会误杀。这个规则能有效防止同一个 QRS 波群因为积分平台的波动被重复检测。提示不要跳过自适应阈值直接写死一个阈值。用固定阈值的心电峰值检测程序换个设备或换个人就得重新调参那不能叫算法只能叫脚本。4. Matlab 全流程落地从读取信号到 QRS 位置导出4.1 读取数据与参数初始化动手前先明确输入数据格式。这个源码包里常见的测试数据有几种MIT-BIH 的.mat文件、.dat原始格式以及普通的csv或txt导出。.mat最方便直接用load就能读进工作区。.dat是 MIT-BIH 的 16 位二进制格式需要配合头文件里的采样率和增益做解析略麻烦。拿到数据后先看一眼基本信息采样率 fs、信号长度、导联数。多导联信号通常选 II 导联或 V1V5 中 QRS 最清晰的一路不要把所有导联直接当多通道数组传入检测函数。确认导联后截取一个 10 秒片段做可视化用眼睛扫一遍信号质量再决定滤波参数。load(ecg_sample.mat); % 假设里面有变量: ecg 和 fs ecg_raw ecg(:); % 统一成列向量 t (0:length(ecg_raw)-1) / fs; % 时间轴 plot(t(1:fs*10), ecg_raw(1:fs*10)); xlabel(Time (s)); ylabel(Amplitude (mV));load之后最好确认一下fs是否存在有些数据集里面采样率叫Fs或者sample_rate不提前查好后面全部按错误采样率算检测结果会全盘错位。这是最不起眼但最致命的低级错误我第一回跑 MIT 数据就栽在这里。4.2 检测主循环滤波、差分、平方、积分、阈值下面是一段可以直接跑的完整检测函数它接受原始心电信号和采样率两个参数返回 R 波峰值对应的索引位置function qrs_index detect_qrs(ecg_raw, fs) % 带通滤波 [b, a] butter(6, [5 15]/(fs/2), bandpass); ecg_f filtfilt(b, a, ecg_raw); % 微分与平方 diff_f diff(ecg_f); diff_f(end1) diff_f(end); % 补位保持长度一致 squared diff_f .^ 2; % 滑动窗口积分 win_width round(0.15 * fs); mwi conv(squared, ones(1, win_width)/win_width, same); % 自适应阈值初始值 est_peak max(mwi(ceil(fs/2):end)); % 取后半段最大值做初始估计 threshold 0.35 * est_peak; % 初始阈值系数 refractory round(0.2 * fs); % 不应期 200ms qrs_index []; last_qrs -refractory; for k win_width1 : length(mwi)-win_width if mwi(k) threshold (k - last_qrs) refractory % 在邻域内寻找精确定位点 [~, peak_pos] max(ecg_f(k-round(0.02*fs):kround(0.02*fs))); qrs_index(end1) k - round(0.02*fs) peak_pos - 1; %#okAGROW last_qrs k; % 动态更新阈值 threshold 0.35 * mwi(k) 0.65 * threshold; end end end这段代码把 Pan-Tompkins 五步走完整串起来了。逐段说明一下带通滤波用 6 阶巴特沃斯filtfilt保证零相位。微分和平方把 QRS 的高斜率特征放大。滑动积分窗口 150ms。初始阈值取积分信号后半段最大值的 35%是因为前半段可能受到滤波边界效应干扰不够可信。遍历时每检测到一个 QRS就用当前积分值按0.35 * current 0.65 * old更新阈值让阈值跟随信号幅值变化。peak_pos那段是精确定位在检测到的粗位置附近 20ms 内找原始滤波信号的最大值把 R 波位置校准到真正的峰顶而不是积分平台的中央。4.3 可视化验证与结果导出检测完成之后一定要画图验证。在原始心电信号上叠加 R 波位置标记一个个看过去比看任何统计指标都直观。导出的结果我一般存成和时间戳对齐的文本文件方便后续和其他算法对比或者喂给心率变异性分析模块。qrs_loc detect_qrs(ecg_raw, fs); qrs_time (qrs_loc - 1) / fs; % 换算成秒 figure; plot(t, ecg_raw, b); hold on; plot(qrs_time, ecg_raw(qrs_loc), rv, MarkerSize, 8, LineWidth, 2); xlabel(Time (s)); ylabel(ECG (mV)); legend(ECG, Detected R peak); grid on;画图这段的qrs_time是从索引换算成秒注意 Matlab 索引从 1 开始所以减 1。如果你处理的信号是从某个时间点开始采集的还要加上起始偏移量。检测结果导出可以用writematrix(qrs_loc, qrs_index.txt)或save(qrs_result.mat, qrs_loc)看下游需要什么格式。导出索引比导出时间更好因为时间可以由采样率事后换算索引则不会引入浮点误差。验证时我习惯把整个信号分成 30 秒一段每段单独目测数漏检和误检。超过 10 分钟的长信号不可能逐点看但可以随机抽几段或者对比已知标注数据算准确率这样比“看着整个曲线觉得不错”可靠得多。5. 峰值检测避坑指南五个实测踩坑与排查方法5.1 预处理阶段的三个高频坑坑一filtfilt 边界效应导致开头出现伪峰现象检测结果中信号开头 0.51 秒内经常出现一个假 R 波而且每次都在同一个位置。原因filtfilt在信号首尾会做边界填充边界过渡段波形严重失真幅值可能异常偏高或偏低经过微分平方积分之后形成一个伪波峰。解决检测有效范围从round(0.2*fs)之后开始前 200ms 直接丢弃不看或者用detrend先粗去基线漂移减轻边界效应。注意不要为了消边界效应把滤波改成filter那会把所有 R 波位置都加上未知的延迟更麻烦。坑二基线漂移过滤不彻底长时程信号检测率暴跌现象前 1 分钟检测挺好往后漏检越来越多积分信号的阈值被慢慢抬高。原因5Hz 高通对 0.3Hz 左右的呼吸漂移衰减有限呼吸起伏的能量被微分平方之后放大积分信号整体斜率往上抬固定比例初始阈值越来越高部分幅值略低的 QRS 就被压制了。解决f_low 从 5Hz 提到 8Hz看检测率变化或者在预处理阶段加一个medfilt1中值滤波估计基线然后减去我实际对比过medfilt1(ecg_raw, round(0.2*fs))估基线再减掉对长时程信号最管用。坑三窗口长度和微分导致的零点偏移对不齐现象检测到的 R 波位置和真实位置偏差几百毫秒肉眼看着标点落在波形旁边而不是峰顶。原因diff把信号缩短了一个点conv(..., same)的欠冲又让积分信号存在群延迟几个误差叠加导致峰值定位偏移。解决diff后手动补尾点conv用same并接受约半个窗口的群延迟最后在精确定位阶段回退到原始滤波信号ecg_f上取邻域最大值不要直接用积分信号的最大值位置。这条不解决你的 RR 间期序列算出来会带系统性偏移影响后续心率变异性分析。5.2 逻辑与调参环节的常见误操作坑四不应期设太短同一个 QRS 被重复检测现象画出图来同一个 R 波上出现两个标记间隔大约一个窗口宽度。原因滑动窗口积分信号是梯形平台平台顶部有轻微波动如果阈值略低于平台峰值可能被判断为两个独立峰值。应该期必须大于等于 QRS 波群总时长加上窗口宽度。解决把 refractory 设成round(0.2 * fs)并且在做检测时同时检查mwi(k) threshold和与上一个检测点的距离大于 refractory。日志里打印每次检测间隔凡是小于 0.2 秒的检出一律剔除观察模式就能定位问题。坑五阈值系数固定 0.35 导致强噪声场景误检现象运动心电或者电极松动片段R 波附近有大量毛刺检测结果出现一串密集假峰或漏检。原因自适应阈值确实会更新但初始阈值系数 0.35 在噪声幅值大的片段里可能低于噪声平台高度让噪声触发检测。解决引入信噪比估计——如果积分信号在某一小段内持续高于阈值但没有满足不应期规则的峰说明噪声平台可能超过有效 QRS这时候把阈值系数提上去比如 0.5。另一个实用技巧是检测前先用rms估计底噪噪声超过信号幅值 30% 的片段整体标记为“低质量”不从里面取 R 波并输出质量指示交给更上层的逻辑去处理。6. 收尾用标注数据量化检测精度别靠目测说“挺好”代码能跑、图能出之后最关键的不是调阈值而是量化评估。我把包内自带标注数据和检测结果做了一次对比才算真正知道这套算法性能到哪个量级。对比逻辑很简单把检测到的每个 QRS 位置与标注位置做匹配在 ±150ms 范围内认为命中。用真阳性TP、假阳性FP、漏检FN三个数算灵敏度和阳性预测值。注意同一标注只能匹配一个检测结果一个检测点只能匹配到一个标注这能避免一个 R 波被多个检测点重复匹配带来的虚高准确率。function [sens, ppv] eval_detection(det_loc, ref_loc, fs) toler round(0.15 * fs); % 匹配窗口 ±150ms matched_det false(size(det_loc)); matched_ref false(size(ref_loc)); for i 1:length(ref_loc) dist abs(det_loc - ref_loc(i)); [~, j] min(dist); if dist(j) toler ~matched_det(j) ~matched_ref(i) matched_det(j) true; matched_ref(i) true; end end tp sum(matched_det); fp sum(~matched_det); fn sum(~matched_ref); sens tp / (tp fn); ppv tp / (tp fp); end这段是我每次换数据、换参数都强制跑一遍的评估脚本。灵敏度和阳性预测值两个数就够了一个看漏检一个看误检。调参的时候拿这两个指标当朴素的反馈信号——把阈值系数从 0.35 调到 0.5如果灵敏度掉了 2 个点但阳性预测值涨了 5 个点那这笔交易就是划算的。用真实标注评估过一遍以后我从那以后拿到任何一份新的采集数据都强制先跑一遍eval_detection再看图再决定要不要动参数。这个顺序反不得先看图容易把注意力放在单段波形上先看指标能直接告诉你参数改动的性价比。这方法帮我避开了无数次“看起来不错换一段就翻车”的尴尬。希望帮到你。本文还有配套的精品资源点击获取