心电信号处理这条路我陆陆续续走了两年多。从最早照着论文抄代码到后来能跑通一个完整的ECG心律失常检测流程中间踩过的坑比写过的代码还多。这篇文章就把我最近一套基于Matlab的心律失常检测实现完整拆开——从预处理、R波检测到分类逻辑再到那些真正让人头疼的调参细节全部记录下来。先说清楚这套东西能干什么输入一段心电信号ECG经过滤波去噪、QRS波群检测主要是R波定位、RR间期计算最后输出心率值并区分出窦性心动过速、心动过缓、室性早搏、房颤这几类常见心律失常。对于做生物医学信号处理课设、毕业设计或者刚入门心电算法方向的朋友这份实现思路可以直接拿去复用。1. 从脏信号到可分析波形ECG预处理的三种噪声与处理方案1.1 心电信号到底有多脏三类噪声与实测现象心电图记录的是心肌细胞电活动在体表的综合电位幅度只有0.5~4mV频率集中在0.05~100Hz这个范围。可实际采到的信号里混着大量比有效信号还要抢眼的东西。我拿MIT-BIH数据库里带噪声的记录试过原始波形的基线能漂出屏幕R波峰高一会儿低一会儿你根本没法直接用固定阈值去检测。主要噪声源有三类处理思路完全不同噪声类型主要来源频率特征对检测的影响基线漂移呼吸、电极片移动、皮肤阻抗变化低于0.5Hz基线上下浮动R波幅值变化固定阈值直接失效工频干扰50Hz电网国内/ 60Hz部分国外数据50/60Hz及其谐波波形变粗毛刺明显峰值检测容易误判肌电噪声肌肉收缩肢体移动、紧张30~300Hz分散频段高频毛刺干扰QRS边界和波形形态判断1.2 预处理三步走去基线漂移、滤工频、降肌电噪声我常用的方案是三段式预处理每段解决一类噪声第一步去除基线漂移。这里有两个选择中值滤波法和高通滤波法。中值滤波的思路是取一个窗口比如0.2秒把信号中值算出来作为该点的局部基线估计再用原始信号减去这个估计。窗口宽度取200~300ms比较合适太短会把QRS成分也滤掉太长又跟不上呼吸引起的慢漂移。不过我更推荐高通滤波用截止频率0.3~0.5Hz的高通滤波器直接把低频漂移切掉。注意截止频率别设到1Hz以上否则ST段会被拉平后续做形态学分析时QRS宽度和ST段信息都会失真。第二步是工频陷波。国内数据用50Hz美国那边MIT-BIH数据库的原始记录是模拟信号数字化里面工频成分不多但如果你自己采集或用了某些公开数据集50Hz/60Hz陷波器基本必做。陷波器用iircomb或者designfilt都能实现关键是品质因数Q值别设太狠Q值过大容易把频谱掏出一个凹坑影响邻近频段。第三步是肌电噪声抑制用低通滤波截止频率100Hz左右就行。为什么不是越低越好因为QRS波群本身含有不少30~50Hz的高频成分尤其Q波的尖端转折你如果一刀切到40HzQRS波峰会变圆滑幅值也会掉后面检测出来的QRS宽度都会失真。低通只是把毛刺抹掉保留波形主体即可。1.3 Matlab滤波代码选型与参数说明以下是我在实际项目里跑通的预处理代码。注意我所有滤波都用了filtfilt做零相位滤波这一点很关键后面专门讲。fs 360; % 采样率MIT-BIH标准是360Hz % 读入原始信号假设ecg_raw是列向量 % 1. 高通滤波去基线漂移截止0.5Hz4阶Butterworth fc_hp 0.5; [b_hp, a_hp] butter(4, fc_hp/(fs/2), high); ecg_base filtfilt(b_hp, a_hp, ecg_raw); % 2. 50Hz陷波去除工频干扰 wo 50/(fs/2); bw wo/35; % Q值控制陷波带宽 [b_notch, a_notch] iirnotch(wo, bw); ecg_notch filtfilt(b_notch, a_notch, ecg_base); % 3. 低通滤波抑制高频肌电截止100Hz fc_lp 100; [b_lp, a_lp] butter(4, fc_lp/(fs/2), low); ecg_clean filtfilt(b_lp, a_lp, ecg_notch);这套预处理做完基线基本平了波形该尖的尖、该圆的圆。我习惯把滤波前后的图叠在一张figure里对比一眼就能看出信号质量有没有救回来。如果滤波后波形还有明显的波浪线多半是高通截止频率给低了或者陷波器没生效去查滤波器参数就好。2. R波检测的经典打法Pan-Tompkins算法Matlab实现详解2.1 R波检测思路把复杂波形简化成一个脉冲R波定位是整个心律失常检测的地基。心率算得准不准、RR间期算得对不对全看这一环。心电图里QRS波群是一个相对高频、高幅值的尖端P波和T波都是低幅值慢变化的波。所以R波检测最经典的思想就是把QRS波群的形态特征放大再把P波、T波和噪声压下去。Pan-Tompkins算法是1985年提出来的经典方法到现在仍然是各类QRS检测器对比的基准。它的核心流程是带通滤波 → 微分 → 平方 → 滑动窗口积分 → 自适应阈值可以说是把信号处理特征提取玩到了极致。每一步的目的都很明确带通滤波5~15HzQRS成分主要集中在这个频段带通之后T波被明显压制P波也残废了。微分突出QRS的斜率变化。R波斜率大T波和P波斜率小微分之后差距进一步拉开。平方所有值变成正值同时让大幅度值更加突出相当于给R波加权重。滑动窗口积分把每个R波变成一段山丘状波形这样可以检测R波的能量包络而不是单个尖峰抗噪声能力强很多。2.2 各环节的参数选定先说带通滤波器。Pan-Tompkins原文用的是5阶Butterworth带通通带5~15Hz我这里直接沿用。低端截止设在5Hz可以过滤掉大部分基线漂移和P波低频分量高端15Hz衰减了部分肌电高频噪声同时保留了QRS频段能量。f_low 5; f_high 15; [b_band, a_band] butter(5, [f_low f_high]/(fs/2), bandpass); ecg_band filtfilt(b_band, a_band, ecg_preprocessed);微分这一步Pan-Tompkins原文用了一个五点差分公式近似等效于4阶差分滤波器对QRS的上升沿和下降沿都比较敏感。我用Matlab的diff函数实现代码是一行但效果等同于五点差分diff_ecg diff(ecg_band); diff_ecg [diff_ecg(1); diff_ecg]; % 长度对齐然后平方squared diff_ecg .^ 2;滑动窗口积分是整个算法里最需要耐心的参数。窗口宽度按文献推荐取150ms换算成采样点win_len round(0.15 * fs)即采样率360Hz时为54个点。窗口宽度起什么作用相当于把各个样本点处的平方值累加起来形成包络。窗口太窄比如小于100ms的话一个R波可能积分出两个峰来——因为Q波的负向尖峰和R波正向尖峰分开太远窗口太宽大于200ms则两个相邻R波尤其心率快的时候会合成一个峰漏检就来了。win_len round(0.15 * fs); moving_sum conv(squared, ones(1, win_len)/win_len, same);2.3 自适应阈值的细节为什么会失灵又是怎么救回来的固定阈值在信号平稳时没问题但真实心电信号幅值会漂。有人深吸一口气心电图QRS幅值能涨30%。这时候固定阈值要么漏检要么把T波甚至是噪声当R波。Pan-Tompkins算法的精髓就在自适应阈值。我实现的是简化版自适应双阈值方案维护两个阈值一个用于检测thres_detect一个是检测到R波后的不应期内峰值记录用于动态调整。在每个R波检测到后用当前峰值更新阈值detected_peaks []; refractory_time round(0.20 * fs); % 200ms不应期防同一R波重复检测 threshold 0.4 * max(moving_sum) 0.1; % 初始阈值 for i 2 : length(moving_sum) - 1 if moving_sum(i) threshold moving_sum(i) moving_sum(i-1) moving_sum(i) moving_sum(i1) % 检查不应期距上一个检测点足够远 if isempty(detected_peaks) || (i - detected_peaks(end)) refractory_time detected_peaks(end1) i; % 用当前信号峰值的均值动态调整阈值防止信号幅值下降后漏检 if moving_sum(i) threshold * 1.2 threshold 0.4 * moving_sum(i) 0.2 * threshold; end end end end这里的关键是不应期refractory period的概念——心电生理上心肌在一次兴奋后约200ms内不会再次产生可兴奋的反应所以一个R波之后200ms内再出现峰值基本不可能是另一个正常QRS心率超过300bpm才可能出现。用200ms这层保险能有效避免因为积分波形振荡导致的同一R波重复计数。实际跑的时候你还会发现一个问题信号前1秒如果没有R波初始阈值设多少我上面用0.4 * max(moving_sum)的思路比较稳——先找整段信号的最大值做参考再乘系数。这个系数我调过很多次0.3偏灵敏会把高耸T波误检0.5偏保守某些低幅值的PVC室性早搏容易被漏掉。0.4算是个折中值后续根据信号质量再微调。2.4 从峰值点到R波位置的修正滑动窗口积分输出的峰位置是积分窗口包的重心附近不是精确的R波尖端位置。所以在积分波形上找到峰后还需要回到原始带通滤波信号或者预处理信号上去找真正的R波尖峰。这个修正步骤很多人会漏掉结果RR间期算出来总有10~20ms的系统偏差心率计算整体偏慢或偏快。我的做法以积分峰的检测点为中心前后各取20个采样点约55ms窗口在带通滤波后的信号ecg_band中找最大值点作为最终R波位置。for k 1 : length(detected_peaks) center detected_peaks(k); search_start max(1, center - 20); search_end min(length(ecg_band), center 20); [~, local_max_idx] max(ecg_band(search_start : search_end)); r_peaks(k) search_start local_max_idx - 1; end做完这一步R波位置就准了。3. 从RR间期到诊断标签心律失常分类的阈值逻辑与形态学特征3.1 先把RR间期和心率算明白R波定位好了RR间期就是一串差分值相邻R波的时间间隔采样点间隔除以采样率。假设r_peaks是N个R波位置向量rr_intervals diff(r_peaks) / fs; % 单位秒 heart_rates 60 ./ rr_intervals; % 单位bpm次/分钟 mean_hr mean(heart_rates);对多数心律失常平均心率是最基础的分类特征。正常人静息心率在60~100bpm之间。窦性心动过速就是心率持续超过100bpm除了运动、发烧这类生理性原因窦性心动过缓则是心率持续低于60bpm。但单看平均心率远远不够我还需要看RR间期的一致性和规律性。医学上心电图诊断会看节律是否规整——正常窦性是基本规整的房颤是绝对不齐的。在工程实现里这个不齐的量用SDNN和RMSSD来刻画sdnn std(rr_intervals); diff_rr diff(rr_intervals); rmssd sqrt(mean(diff_rr .^ 2));SDNN是全部RR间期的标准差正常人在50ms以下偏高说明整体变异性大。RMSSD是相邻RR间期差值的均方根它对快速的逐拍变化更敏感房颤时RMSSD会显著升高通常大于100ms。为什么房颤RMSSD会特别高因为房颤时心房无序放电心室率完全不规则有的RR间期只有0.4s心率150有的1.2s心率50相邻拍之间的差值特别大。这个特征用SDNN和RMSSD两个指标叠加很有辨识度。3.2 常见心律失常的判别规则我在这套实现里做了四类判别窦性心动过速、窦性心动过缓、室性早搏PVC、房颤粗筛。判别规则用表格整理如下类型心率特征RR间期特征其他形态特征窦性心律60~100bpmRR相对规整SDNN50ms每个QRS前有P波窦性心动过速100bpmRR仍相对规整P波存在心率快窦性心动过缓60bpmRR仍相对规整P波存在心率慢室性早搏PVC可变的提前出现一个RR间期短于平均随后代偿间歇QRS宽120ms形态异形房颤多变绝对不齐RMSSD100msP波消失基线上有f波这段分类逻辑用代码写出来大约是if mean_hr 100 label 窦性心动过速; elseif mean_hr 60 label 窦性心动过缓; elseif rmssd 100 sdnn 80 label 房颤初步判断需结合P波; else label 正常窦性心律; end这种阈值逻辑确实比较工程化医学诊断不会这么简单但用来做算法验证和课设完全够用。如果在实际数据分析场景里我建议把分类做成概率输出而不是硬标签比如基于特征向量算距离得一个置信度而不是一个if-else一锤定音。3.3 室性早搏单独拉出来说R波并不总是高大宽PVCPremature Ventricular Contraction室性早搏在心律失常检测里很特殊因为它在心电图上的表现太有辨识度了提前出现、QRS宽大畸形120ms、其前没有相关的P波、后有代偿间歇。工程上实现的关键有两个——提前量和QRS宽度。提前量当前RR间期明显短于平均RR间期的80%同时下一拍代偿间歇明显长于平均。QRS宽度在带通滤波后的信号上找到R波峰值后往前找Q波起点、往后找S波终点算两点间距。QRS正常宽度不超过120ms即每拍0.02s对应60×0.021.2个小格。PVC的QRS宽大往往超过140ms。Matlab实现QRS宽度提取qrs_width_ms zeros(size(r_peaks)); for k 1 : length(r_peaks) r_pos r_peaks(k); % 起点R峰前40ms内第一个局部极小值 search_start max(1, r_pos - round(0.04 * fs)); [~, q_idx] min(ecg_clean(search_start : r_pos)); q_start search_start q_idx - 1; % 终点R峰后60ms内第一个局部极小值 search_end min(length(ecg_clean), r_pos round(0.06 * fs)); [~, s_idx] min(ecg_clean(r_pos : search_end)); s_end r_pos s_idx - 1; qrs_width_ms(k) (s_end - q_start) / fs * 1000; end这里有个细节找局部极小值用的是预处理后的ecg_clean不是ecg_band。因为带通滤波后QRS边缘信息丢失Q波起点和S波终点会被平滑掉宽度测量不可靠。预处理信号里同时保留了低频基线和高频边缘测量宽度更接近真实值。这个差别你可以在同一段数据上对比一下非常直观。3.4 P波检测分类规则的隐藏难度上面表格里提到P波消失和P波存在实际做的时候P波检测比R波检测难得多。P波幅值只有0.1~0.2mV频率更低0.5~10Hz噪声稍大一点就淹没了。我在这套简单实现里没有做完整的P波形态分析但是一个替代方案可以告诉大家用检测到的R波位置回溯在RR间期前段前1/3处搜索局部最大值超过一定阈值就认为存在P波。% 对每个RR间期在R波之前约80~200ms窗口内检测P波 p_width round(0.12 * fs); % P波检测窗口宽度 p_detect_count 0; for k 2 : length(r_peaks) win_start r_peaks(k) - round(0.2 * fs); win_end r_peaks(k) - round(0.08 * fs); if win_start 1, continue; end segment ecg_clean(win_start : win_end); % 粗略判断是否有P波窗口内最大峰值超过平均幅值的0.15 p_detect_count p_detect_count (max(segment) 0.15 * max(ecg_clean)); end p_ratio p_detect_count / (length(r_peaks) - 1);房颤时p_ratio会显著偏低P波消失基线上是锯齿状f波这个指标可以辅助分类准确率有所提升但还不够稳健。如果想做专业一点的房颤检测建议用RR间期的熵分析或者马尔可夫模型那是另一个大课题。4. 没数据怎么练手公共数据库与仿真心电图信号4.1 MIT-BIH心律失常数据库最经典的练手数据聊实现时我默认你手里已经有心电数据了但很多读者第一步就卡在去哪弄数据。最经典的选择是MIT-BIH Arrhythmia DatabaseMIT-BIH心律失常数据库这是心电图算法验证的标准场各大论文、网上的开源项目几乎都用它。它包含48条记录每条30分钟采样率360Hz两位医师独立标注了每个拍子的类型N正常、V室性早搏、A房性早搏等标签。获取方式PhysioNet官网注册后可以下载有record.dat、annotations等多种文件格式。Matlab读起来稍微麻烦因为官方数据不是Matlab格式。我推荐两种做法第一种装PhysioNet官方的WFDB Toolbox用rdsamp直接读信号。在Matlab里一行[ecg_raw, Fs] rdsamp(100, 1); % 读取100号记录的第一导联第二种其他人打包好的.mat文件在GitHub上也能找到直接load进来用。两种方式我建议学第一种能把读取逻辑掌握在自己手里。4.2 没有真实数据时的替代方案仿真ECG信号如果只是想验证算法流程不想去下载数据库也可以在Matlab里生成模拟心电信号。有个免费的Open Source ECG Generatorecg函数可以在File Exchange上找它基于动力学模型生成包含P波、QRS、T波的合成心电采样率可设还能加噪声。仿真信号对调试流程很好用尤其是你要测试滤波器的效果——因为你可以控制噪声类型和幅度滤完一眼就能看出效果。我还是建议后期一定要换真实数据跑一遍。仿真信号太干净了R波高度整齐、P波一致算法在上面跑得漂亮不代表真实数据也能行。我见过不少同学仿真数据准确率99%换成MIT-BIH直接掉到85%甚至更低原因就是真实信号的复杂度和噪声类型仿真里根本没有。5. 调参和踩坑实录我在这套流程里翻过车的几个环节5.1 阈值失灵的真相信号质量突变最开始我跑MIT-BIH 201号记录时R波检测准确率惨不忍睹。我调了半天滤波器又试了各种阈值系数最后才发现问题不在算法而在数据——201记录前一段信号质量极差幅值骤降R波几乎和噪声融为一体。用整段数据的最大值做初始阈值时前30秒的R波全部低于阈值一个都检不出来。解决办法很土但有效把信号按5秒一段切分每段单独计算自适应阈值基数再在段与段之间做阈值平滑过渡。这样即使某一段信号质量差也只是局部检测不准不会把全局阈值带偏。另外一定要看检测结果的波形总览图。把原始波形和检测到的R波位置画在同一张图上用红色圆圈标出R位置眼睛一扫就知道哪里漏检了、哪里误检了。这个习惯我到现在都在用也是所有排错的第一步。5.2 滤波相位失真为什么必须用filtfilt滤波器有固有相位延迟普通filter函数会带来非线性相位畸变——滤波出来的波形整体平移了一段距离而且不同频率成分平移量不一样波形形状都会改变。R波检测对时间位置精度要求高RR间期差了1个采样点就是2.8ms误差所以必须用filtfilt做零相位滤波。filtfilt的原理是对信号正向滤波一次再反向滤波一次相位偏移正好抵消代价是计算量翻倍加边界效应。Matlab里实现很简单但很多人会顺手写成filter结果R波位置普遍滞后十几个采样点心率倒是没差多少RR间期差分抵消了但只要做过R波对齐分析就会露馅。5.3 采样率统一360Hz不是默认值MIT-BIH的标准采样率是360Hz但有些数据库或设备输出是250Hz、500Hz甚至128Hz。采样率不同length相关的参数全都要按比例换算哪里要用round(0.15 * fs)这种写法就不要写死一个54。有些傻瓜代码里写死了窗口长度为54一换数据就出问题。我处理其他设备数据前第一件事就是先确认fs再批量换算所有时间相关参数。此外重采样也是个坑。如果需要在不同采样率数据之间做比较用resample函数重采样到统一频率。但注意重采样会引入一定插值误差R波的幅值和位置都有细微变化对需要精确宽度的任务敏感最好在重采样之前做完形态学分析。5.4 评估指标别只报准确率做心律失常检测评估时很多人只报一个准确率99%。在类别不平衡的情况下数据里正常拍占95%准确率就是误导指标。一个永远输出正常的检测器在只有5%异常拍的数据上准确率也是95%。所以在评估时我会同时看这几个指标敏感度Sensitivity真阳性率异常样本中检出的比例漏检影响的就是它。阳性预测值Positive Predictive ValuePPV检出的阳性中真正异常的比例误检影响的就是它。F1分数敏感度和PPV的调和平均。拿R波检测来说敏感度低说明漏检多PPV低说明误检多把T波当R波。我实际跑下来质量好的记录如100、101、103号F1能到99%以上而质量差的记录如105、201、207号F1就只有90%~93%。这个分布是正常的拿到一条数据先看它的F1落在这个范围哪一端能快速判断是你的算法差还是数据本身难。6. 向前一步从离线检测到实时分析的优化空间这套实现整体是离线处理逻辑但很多实际应用——可穿戴心电监测、运动手环预警——要求实时或准实时分析。如果想把这套逻辑往实时方向推进有几个明确的优化点第一流式处理替代整段处理。现在代码里大量使用max、std这适合整段分析但实时系统只能看过去几秒的数据。改成滑动窗计算窗口长度取10秒每1秒更新一次实时性和稳定性之间能取得平衡。第二阈值初始化改成学习模式。设备刚开机时前几秒用较高的阈值并且更严格地确认R波比如要求连续两个满足条件的峰才算确认之后逐渐过渡到正常检测模式。这样可以避免开头信号质量差导致阈值被污染。第三内存和运算效率。处理两小时的长时心电记录时conv加findpeaks这种整段操作会占用不少内存。实时系统里建议用dsp.System objects结合while循环逐段处理。Matlab的DSP System Toolbox里有专门为流式信号设计的高通滤波器对象配合定时器可以做到每200ms处理一帧数据。第四个方向是深度学习替代手工特征这是另一个大坑。LSTM网络可以直接输入原始波形输出心律失常类别效果在某些数据集上优于传统阈值逻辑但可解释性差、调参成本高而且对数据量和计算资源要求不低。如果是为了做毕设或者实际产品验证我建议先把传统特征工程这套搞清楚再做神经网络对比两条腿走路。回到我对这套流程的整体评价逻辑不复杂代码量也不大但真正把它跑得稳定、结果讲得清楚需要对每个参数背后的物理和生理含义有理解。心电信号处理是个很有意思的方向R波检测做扎实了心率变异性分析、呼吸频率估算、情绪识别这些衍生应用都能顺带解锁。我把这份完整的Matlab实现流程记录在这里希望对正在入门ECG方向的你有所帮助。