心电数据拿到手第一步不是急着上模型而是先想清楚你的信号到底有多“脏”。这是我做心电信号处理这几年最深的一个体会。这次分享一个用 Matlab 实现的心电图ECG心律失常检测流程从去噪、R 波定位到心律失常规则分类完整走一遍。这套流程适合生物医学工程、信号处理方向的初学者也适合做课程设计和大作业的同学参考——你不需要先把深度学习啃完用经典信号处理思路就能跑通一个可解释、可复现的检测方案。1. 整体设计思路为什么先用规则检测而不是直接上深度学习很多新手拿到 ECG 数据第一反应是“我是不是该用 LSTM、CNN 分类”。但真实的心电图分析场景里信号质量参差不齐标注成本也很高深度学习模型在这种数据上容易过拟合而且出问题后你很难解释清楚是哪一步出了错。经典信号处理路线的逻辑恰好相反先告诉你“信号怎么来、噪声怎么去、R 波在哪”再基于生理规则做分类。这套流程每一步都可视化、可调试遇到异常你能定位到具体的中间环节。我的整体管线分为四段数据读取滤波预处理QRS 波检测核心是 R 波定位特征提取和心律失常规则判断。每一段都有对应的 Matlab 函数组合起来就是完整的检测系统。为什么把 QRS 检测放在最关键的位置因为心律失常分类的特征无论是心率快慢、节律是否规整、有没有早搏都高度依赖准确的 R 波位置。R 波定位一旦偏差超过 20 ms后面算出来的 RR 间期、心率变异性指标全部失真。所以这一步值得花最多时间打磨。2. 环境和数据准备2.1 数据来源与读取方式我做实验用的标准数据集是 MIT-BIH 心律失常数据库里面包含 48 条双导联心电记录每条记录时长 30 分钟采样率 360 Hz幅值单位经过换算后是 mV 量级。这是心电图算法研究的事实标准你论文里写“用 MIT-BIH 验证”审稿人不会质疑数据基础。读取 MIT-BIH 数据用的 .dat 和 .hea 文件在 Matlab 里可以直接用开源的rddata.m脚本或者用 WFDB Toolbox。一段简单的读取逻辑是这样% 读取 MIT-BIH 记录以第100条为例 [signal, Fs, tm] rdsamp(100, 2); % 读取两个导联 ann rdann(100, atr); % 读取标注 % signal 是 N×2 的矩阵第一列通常为 MLII 导联幅值单位 mV Fs 360; % 采样率 360 Hz读取之后第一件事是观察信号质量。我习惯用plot把原始信号画出来看一眼基线漂移的程度、有没有明显的工频干扰再决定滤波器的参数。不要上来就滤波先肉眼观察是最快的诊断方式。2.2 信号长度的截取正片处理的时候不需要一次吃下整条 30 分钟数据那样既不直观也容易把内存撑爆。一般做法是滑动窗口处理窗口长度 10 秒或者 20 秒步长 5 秒。10 秒窗口足够容纳 10 个以上的心跳周期算心率时也有足够样本而 5 秒的步长保证相邻窗口之间有重叠防止 R 波正好落在窗口边缘被截断。winLen 10 * Fs; % 窗口长度单位采样点 stepLen 5 * Fs; % 滑动步长 for idx 1:stepLen:(length(signal) - winLen 1) seg signal(idx:idxwinLen-1, 1); % 取第一个导联 % 后续处理... end3. 预处理降噪三步把信号整干净3.1 基线漂移的去除心电信号里最常见的干扰之一就是基线漂移表现为整条曲线上下缓慢晃动频率通常在 0.5 Hz 以下。产生原因包括呼吸运动、电极移动和放大器漂移。基线漂移会让后续的幅值阈值判断出错必须去掉。我用的方法是高通滤波截止频率设在 0.5 Hz。Matlab 里设计一个 4 阶 Butterworth 高通滤波器很简单fc_high 0.5; % 截止频率 0.5 Hz [b, a] butter(4, fc_high/(Fs/2), high); ecg_high filtfilt(b, a, ecg_raw);注意必须用filtfilt做零相位滤波而不是filter。filter会引入相位延迟导致滤波后的 R 波位置和时间偏移直接影响后续检测精度。filtfilt正向反向各做一次相位延迟抵消零相位失真。3.2 工频干扰的滤除50 Hz 或 60 Hz 的工频干扰几乎是心电采集绕不开的问题特别在非屏蔽环境下信号上会叠一层细密的“毛刺”。去除工频干扰用带阻滤波器陷波器品质因数 Q 别设太高太高了滤波器可能不稳定也会把相邻频率的有用分量削掉。fc_notch 50; % 国内工频 50 Hz bw 2; % 带宽 2 Hz [b, a] iirnotch(fc_notch/(Fs/2), bw/(Fs/2)); ecg_notch filtfilt(b, a, ecg_high);如果你观察到的工频干扰并不明显这步可以跳过。经验法则是先看频谱图如果 50 Hz 处的谱峰没有明显突出就不需要陷波少一个环节少一分信号失真。3.3 肌电和高频噪声的压制肌电噪声来自肌肉收缩频率范围比较宽经常和心电信号的高频成分混在一起。为了突出 QRS 波群QRS 的频率能量主要集中在 520 Hz 左右我会再加一个低通滤波器截止频率设在 3540 Hz。fc_low 35; [b, a] butter(4, fc_low/(Fs/2), low); ecg_clean filtfilt(b, a, ecg_notch);这里有个取舍滤波越狠信号越平滑但 QRS 波的峰值幅度会被压低P 波和 T 波也可能被削弱。我试过把低通截止频率降到 25 HzR 波幅度确实掉了将近 15%导致后面的自适应阈值被迫重新标定。一般的处理建议是保留 35 Hz 以上做 QRS 检测额外做一条 15 Hz 低通的支路单独分析 P 波和 T 波两条线路互不干扰。这是我自己折腾几轮下来最稳的架构。4. QRS 波检测Pan-Tompkins 算法实现4.1 算法原理拆解Pan-Tompkins 算法是 QRS 检测的经典方案1985 年提出到现在还被广泛使用。它利用的是 QRS 波群的三个特点斜率陡、幅度高、宽度窄。算法流程是带通滤波 → 差分 → 平方 → 滑动窗口积分 → 自适应阈值。带通滤波我们在前面已经完成了实际做的时候我一般把 515 Hz 的带通单独再做一遍给检测器用而不是用预处理后的完整信号。原因是 QRS 的频率能量集中在这个区间先滤掉其他成分后续的差分和阈值判断受到的干扰会小很多。差分步是抓住 QRS 波的斜率特征一阶差分相当于求导让斜率陡峭的 QRS 波在差分信号里形成更明显的峰值。平方操作把差分结果全部转化成正数同时放大幅度差异让高幅值的 QRS 更突出。滑动窗口积分则是把一段窗口内的平方值累加起到平滑作用窗口长度我取 150 ms相当于 54 个采样点Fs360这个长度大约是 QRS 波宽度的两倍效果最理想。核心代码段长这样% 差分 diff_sig diff(ecg_band); % 一阶差分 diff_sig [0; diff_sig(:)]; % 补一位对齐长度 % 平方 squared diff_sig .^ 2; % 滑动窗口积分 winSize round(0.15 * Fs); integral filter(ones(1, winSize), 1, squared);4.2 自适应阈值和不应期处理阈值怎么定是检测器好不好用的关键。固定阈值在信号稳定时能用但遇到振幅变异大的心电片段就会出问题比如早搏后的代偿间歇期R 波幅度变化很大一个固定阈值很难两头兼顾。Pan-Tompkins 的解决方案是自适应阈值每个 R 波检测周期都根据最近一段信号的峰值更新阈值。我这里使用了双阈值两个阈值分别取信号峰值的 0.3 和 0.15高阈值定位确定的 R 波低阈值用于在高阈值没有输出时做补充检测。同时加了一个 200 ms 的不应期防止 T 波被误判成 R 波。生理上心脏在一次心搏后有一段约 200 ms 的绝对不应期这段窗口内不可能再出现新的 QRS 波所以检测器也应该尊重这个生理限制。thr_high 0.3 * max(integral(tRange)); thr_low 0.15 * max(integral(tRange)); refractory 0.2 * Fs; % 200 ms 不应期 % 峰值检测 [pksHigh, locsHigh] findpeaks(integral, MinPeakHeight, thr_high, MinPeakDistance, refractory); [pksLow, locsLow] findpeaks(integral, MinPeakHeight, thr_low, MinPeakDistance, refractory); % 合并两次检测结果以高阈值优先低阈值补充漏检是检测器最烦的问题特别是遇到一些幅度特别小的 QRS。我在自适应阈值之外加了一个回溯机制连续两个 RR 间期超过上一个 RR 间期均值的 1.5 倍时判定可能漏检了此时在当前积分数值里寻找高于 0.1 倍峰值的小峰作为漏检的 R 波补回来。这个回溯在 bradycardia心动过缓样本里尤其管用因为心率慢的时候 R 波形态差异很大固定逻辑经常漏。4.3 检测结果的可视化与评估算法跑完必须把 R 波位置标注在原始信号上肉眼看一遍确认没有多检漏检。这个环节我用一段简单的代码处理figure; plot(t, ecg_clean, b); hold on; plot(t(locsHigh), ecg_clean(locsHigh), ro, MarkerSize, 6); plot(t(locsLow), ecg_clean(locsLow), g, MarkerSize, 8); xlabel(时间 (s)); ylabel(幅值 (mV)); title(R波检测结果); legend(滤波后信号, 高阈值检测, 低阈值补充);对照数据库的标注文件可以用灵敏度Se和阳性预测率PPV两个指标评估好坏。指标低于 99% 的时候不要急着调模型先看漏检误检发生在什么波形特征上再针对性调阈值。5. 特征提取与心律失常分类规则5.1 核心特征RR 间期与瞬时心率R 波位置确定之后特征计算就顺理成章了。RR 间期就是相邻两个 R 波之间的时间间隔单位是秒。瞬时心率等于 60 除以 RR 间期单位 bpm。这些数值直接反映了心脏的节律状态。rrIntervals diff(locsR) / Fs; % RR间期单位秒 heartRate 60 ./ rrIntervals; % 瞬时心率单位bpm avgHR mean(heartRate); % 平均心率不同年龄段人群的正常静息心率范围大约在 60100 bpm。我在算法里的判定标准如下类型判定条件备注窦性心律RR 间期偏差 120 ms心率 60100正常窦性心动过缓平均心率 60 bpm常见于运动员窦性心动过速平均心率 100 bpm运动、紧张、发热等心律不齐RR 间期标准差 120 ms可进一步分析早搏PVC异常提前的 QRS 波形态宽大需进一步形态判断这个表看着简单但每一行的判据背后都有讲究。比如心律不齐的判断要算连续 30 秒以上的 RR 序列标准差只看几个心动周期很容易误判。5.2 早搏PVC的形态分析PVC 也就是室性早搏在心电图上表现为提前出现的宽大畸形 QRS 波时程一般超过 120 ms其后有代偿间歇。我在分类逻辑里专门加了一个 PVC 判别分支核心依据是 QRS 波宽度和形态相关性这两个特征。QRS 波宽度的计算以 R 波位置为中心向左向右回溯到波形起点和终点取起点和终点的时间差。形态相关性使用当前 QRS 段和正常 QRS 模板的相关系数相关系数低就说明形态偏离正常。具体实现qrsWidth zeros(length(locsR), 1); for i 1:length(locsR) left locsR(i) - round(0.04*Fs); right locsR(i) round(0.08*Fs); [~, idxStart] min(ecg_clean(left:locsR(i))); [~, idxEnd] min(ecg_clean(locsR(i):right)); qrsWidth(i) (idxEnd locsR(i) - left - 1 - idxStart) / Fs * 1000; % ms end % QRS宽度大于120ms且显著提前时判为PVC pvcIdx find(qrsWidth 120 (rrIntervals 0.8*nanmean(rrIntervals)));这个形态模板的做法有个前提就是信号里要有足够多的正常心搏来生成模板。实操中我取前 10 个检测到的 QRS 波的平均值作为初始模板之后每 5 秒更新一次。如果你处理的样本里早搏频率很高模板会被污染这时候要固定初始模板不更新或者用中位数滤波来挑选正常心搏。5.3 分类决策流程把上述规则整合起来就形成一个明确的决策树先看心率再看节律是否规整最后看形态是否异常。if avgHR 60 type 心动过缓; elseif avgHR 100 type 心动过速; else type 正常; end % 节律规整性判断 rrStd std(rrIntervals); if rrStd 0.12 type [type 合并心律不齐]; end % 早搏判断 if ~isempty(pvcIdx) type [type 合并PVC早搏]; end这套决策树写起来很直白但别小看它的实用价值。我用它在 MIT-BIH 的 100、101、103、105、109、111、115 这些样本上做过测试对正常和心动过速、心动过缓的识别准确率超过 97%PVC 的识别准确率大约 93%。相比花大量时间调深度学习模型这个规则系统胜在透明可靠。5.4 快速原型找对工具能省一半时间如果你想快速验证算法在某个数据集上的表现Matlab 自带的一些函数是很好的捷径。findpeaks的参数里MinPeakDistance对应不应期MinPeakHeight对应高度阈值配合islocalmax或isoutlier能快速做出一个粗糙版本。我的建议是先用这些内置函数跑通流程验证逻辑再回头替换成自己写的 Pan-Tompkins这样调试效率最高。6. 常见问题与排查技巧实录心电检测领域有个著名的比喻“所有算法在干净数据上都在工作真正的差距出现在噪声和异常波形面前。”我做这套流程时踩过不少坑挑典型的列出来都是能直接落地的经验。6.1 基线漂移滤除后 R 波幅度也被削了这个现象最常见。原因是我把高通滤波器的阶数设得太高比如 6 阶以上或者截止频率设到了 1 Hz。高频分量在滤掉基线漂移的同时也把 QRS 波的斜率和幅度削掉了一部分。解决方法是严格按 0.5 Hz 截止频率、4 阶以内设计另外用filtfilt避免相位偏移。6.2 低阈值补充检测框了一堆 T 波低阈值的设计本来是为了提高灵敏度但如果信号里有高耸的 T 波特别是在心率较快时T 波峰值很容易触碰到阈值。一个实用的技巧是把不应期从 200 ms 调整为 300350 ms因为正常心率下 RR 间期不可能短于 300 ms 太多。另一个技巧是给 T 波加一个形态约束T 波的积分峰值上升斜率普遍小于 QRS可以比较积分峰两侧的陡峭程度。6.3 早搏后的误检PVC 之后往往跟着一个较长的代偿间歇间隔之后的下一次心搏 R 波可能因为没有充分恢复到基线幅度被压低导致漏检。我在处理后发现最有效的方案是把阈值更新逻辑改成“滑动窗口取中位数”中位数比均值对单个异常值的抵抗力强很多。在 RR 间期序列上用medfilt1(rr, 5)做一次中值滤波能够有效把这种长间歇的误判修正过来。再看几个高频问题的速查表问题现象可能原因排查与解决R 波检测灵敏度低高通截止频率太低检查滤波参数适当提升至 0.51 Hz误检率高多检 T 波不应期太短/阈值太低加大不应期至 300 ms提高低阈值心率波动大滑动窗口长度太短窗口至少 10 秒RR 序列中值滤波早搏处漏检自适应阈值更新过快将阈值跟踪改用中位数或加长更新间隔分类结果和标注不一致窗口缩进丢失心搏用 5 秒重叠的滑动窗口保证覆盖7. 从流程到系统的扩展思路如果你把这条流程跑通了恭喜你对心电信号的理解已经比很多人扎实了。接下来还有几个值得探索的方向。一是把 MATLAB 流程封装成函数做成一个带 GUI 的心电分析小工具用App Designer搭界面导入数据、显示波形、输出诊断结果界面化和流程化之后在课程答辩里的演示效果会好很多。二是引入小波变换做 QRS 检测的替代方案。Pan-Tompkins 算法在明显噪声环境下已经很稳但遇到 QRS 波极性反转或幅度极低的场景比如某些病理性 Q 波早搏小波变换的多尺度分析只锁定特定频率的子带检测结果会更抗干扰。三是加上更多类别的心律失常分类。现在只有心动过速、心动过缓、心律不齐、PVC 早搏这四个类别如果要识别房颤、房扑、束支传导阻滞还需要在特征里加入 P 波检测判断 P 波与 QRS 波的关系。这个工作量不小但每一步都有明确的信号处理基础跟着学下来收获会很大。最后再说一点我的切身体会做心电图分析这类医学信号处理最重要的不是堆算法而是对信号本身保持敬畏。每一个 R 波标注都对应一次真实的心脏搏动如果你从没对着原始波形逐段检查过检测结果就别轻易说这个模型好。把数据可视化做好把每一步的中间产物看清楚这套流程才能真正为你所用。