简介这份资源是面向生物医学工程、电子信息及相关专业学生的心电信号分析课程设计参考资料以PDF形式呈现课题二“基于MATLAB平台的心电信号分析系统设计及仿真”的完整方案适合正在做数字信号处理课程设计或信号与系统综合实验的读者参考。资源包共1个PDF文件体积约1019KB内容围绕MATLAB程序设计展开。文档完整梳理了从MIT-BIH数据库读取数字心电信号、还原实际波形、对非均匀采样数据做线性插值到依据心电信号频域特点设计低通与高通滤波器、滤除基线漂移与工频干扰再到对滤波前后信号做频谱分析与比较的整套流程并给出Simulink动态建模与仿真的设计思路附录部分还整理了心电信号读取原理、滤波器指标选取与课程设计报告撰写要求。目前已有357人学习可为选题开题、程序编写与结果分析提供一份可对照的工程实例。1. 心电信号分析系统为什么值得用 MATLAB 从零搭一遍做信号处理方向课题的人往往先找现成工具箱跑一遍图出来挺漂亮可一旦被追问“这个 R 波凭什么标在这个位置”“把滤波器阶数从 3 改成 5 结果会怎样”就答不上来了。基于 MATLAB 平台把心电信号分析系统从数据读取、噪声去除、QRS 检测一路搭到分类识别和仿真界面暴露出来的问题会非常具体采样率不同时间参数要重算滤波器相位没处理 T 波就会整体平移自适应阈值不收敛得去翻积分窗而不是改检测逻辑。这套系统设计的输入是一条原始心电记录输出是干净波形、R 波位置序列、心率与心率变异性指标外加一个点一下就能跑的仿真界面。它覆盖了生理信号处理里最典型的那几类噪声和评价方法适合生物医学工程、电子信息方向的课程设计或毕业设计也适合想练一练 MATLAB 工程化写法的开发者。2. 心电信号读入、噪声建模与 MATLAB 工程骨架设计心电信号分析系统的第一道坎不是算法而是数据从哪来、长什么样。公开的心电数据里MIT-BIH 心律失常数据库是课程设计用得最多的一份它每条记录由三个文件组成.hea描述采样率、导联和增益.dat存二进制采样值.atr存专家标注的心拍位置和类型。搞清楚这三个文件的对应关系后面所有评价指标才有参照物。2.1 MIT-BIH 记录在 MATLAB 里的读取与通道选择常见做法是装 WFDB Toolbox for MATLAB用它提供的rdsamp和rdann两个函数把信号和标注读进来。这两行代码看着简单但通道选错、时长截错是新手最容易踩的坑。% 读取 MIT-BIH 100 号记录的前 10 秒 % 该记录采样率 360 Hz10 秒对应 3600 个采样点 [signal, Fs, tm] rdsamp(mitdb/100, 1:3600); % 读取同一记录的专家标注anntype 里 N 表示正常心拍V 表示室早 [ann, anntype] rdann(mitdb/100, atr);rdsamp返回三个量signal是 N×C 的矩阵C 是导联数MIT-BIH 多数记录是两导联第一导联通常是 MLIIFs是采样率100 号记录是 360 Hztm是时间列单位秒。rdann返回的ann是以采样点为单位的 R 波位置索引anntype是与之等长的字符数组。这里有个细节要注意rdann读出来的是整条记录的标注如果只取了前 3600 点做分析必须自己按索引范围裁剪否则参考位置和检测结果对不上。2.1.1 采样率对后续所有时间参数的影响采样率决定了每一个以秒为单位的时间参数要换算成多少个点。360 Hz 下 150 ms 积分窗是 54 个点1000 Hz 下就是 150 个点。很多人在自己录的数据上跑得好好的代码换到 MIT-BIH 上就漏检一堆原因就是积分窗、不应期这些参数还按原来 1000 Hz 写的。噪声/干扰主要频段常用处理关键参数基线漂移0.05–1 Hz高通滤波或中值滤波截止 0.5 Hz中值窗 0.6 s工频干扰50 Hz ±1 HzIIR 陷波Q 值 30–40肌电噪声20–100 Hz低通滤波截止 40 Hz电极接触噪声突变阶跃中值滤波或形态滤波窗长 0.2 s2.2 工频干扰、基线漂移与肌电噪声的滤波链处理顺序上我一般按“先去基线漂移再去工频最后压高频肌电”来排。原因是基线漂移幅度能达到毫伏级先不处理会让后面所有阈值判断全部失真。用 Butterworth 带通一次性把 0.5 Hz 以下和 40 Hz 以上切掉再用 IIR 陷波打掉 50 Hz这是最省事也最稳的组合。% 0.5-40 Hz 带通切掉基线漂移和高频肌电 [b_bp, a_bp] butter(3, [0.5 40]/(Fs/2), bandpass); ecg_bp filtfilt(b_bp, a_bp, ecg_raw); % 50 Hz 工频陷波Q 值决定陷波宽度 wo 50/(Fs/2); bw wo/35; [b_nc, a_nc] iirnotch(wo, bw); ecg_clean filtfilt(b_nc, a_nc, ecg_bp);butter的第一个参数是阶数3 阶在绝大多数场景够用阶数往上加到 5、6 会有更陡的过渡带但相位非线性也更严重。filtfilt做的是前向加后向两次滤波等效零相位代价是有效阶数翻倍、边界会有短暂瞬态。如果这里图省事用filter输出的 QRS 波在时间轴上会整体滞后十几到几十毫秒后面跟专家标注比对时就会系统性错位。iirnotch的第二个参数是归一化带宽写成wo/35等价于 Q 值 35。Q 越大陷波越窄对 50 Hz 附近有用信号损伤越小但对工频频率漂移的容忍度也越低。国内电网频率相对稳定Q 取 30–40 是合理区间。2.3 工程目录结构与函数接口约定系统设计到后半程代码会膨胀到十几个文件这时候目录怎么分比算法本身更影响调试效率。我习惯按“数据—处理—检测—特征—界面”五层拆目录data/放原始记录和缓存src/preprocess/放滤波函数src/detect/放 QRS 相关src/features/放特征提取和分类src/gui/放界面回调results/落盘中间结果和评价指标。每个函数只干一件事输入输出用结构体串起来方便单独拿一段信号复现问题。function out ecg_preprocess(raw, Fs, cfg) %ECG_PREPROCESS 心电信号预处理带通 陷波 % raw 原始信号 N×1Fs 采样率cfg 可覆盖默认参数 % 返回 out.clean 去噪信号、out.fs 采样率、out.cfg 实际配置 arguments raw (:,1) double Fs (1,1) double {mustBePositive} cfg struct struct(bp,[0.5 40],notch,50,order,3,Q,35) end [b, a] butter(cfg.order, cfg.bp/(Fs/2), bandpass); x filtfilt(b, a, raw); wo cfg.notch/(Fs/2); [bn, an] iirnotch(wo, wo/cfg.Q); out.clean filtfilt(bn, an, x); out.fs Fs; out.cfg cfg; end用arguments块声明参数好处是配置项集中在一处界面传参时不用记位置还能顺便做类型和正数校验。这里把cfg设成带默认值的结构体界面上想改截止频率只需要传struct(bp,[1 35])覆盖一项其余保持默认。提示把每次运行实际生效的cfg一起存进结果文件过两周回头复现某个图的时候能省掉大量猜参数的时间。3. QRS 波检测、心率计算与阈值参数整定QRS 检测是整个心电信号分析系统的枢纽。R 波位置一旦错心率、RR 间期、HRV 指标、心拍分类的切窗位置全都会跟着错。Pan-Tompkins 之所以成为事实标准是因为它用一套纯时域的因果处理链能在嵌入式环境里实时跑同时自适应阈值让它对幅度变化有不错的鲁棒性。3.1 Pan-Tompkins 五步法在 MATLAB 中的逐步实现五步依次是5–15 Hz 带通、微分、平方、滑动窗积分、自适应双阈值。带通这一步是为了让 QRS 的 10 Hz 左右主频段凸出来同时把 P 波和 T 波压下去。微分把波形的斜率信息放大平方让所有值非负并进一步拉开 QRS 与其他波的差距积分把能量摊平成一个包络。% Pan-Tompkins 前四级处理 [b, a] butter(3, [5 15]/(Fs/2), bandpass); x filtfilt(b, a, ecg_clean); % 第一步带通 d [0; diff(x)]; % 第二步一阶差分近似微分 sq d.^2; % 第三步平方非负且放大能量差 win max(1, round(0.150 * Fs)); % 第四步150 ms 积分窗 kernel ones(1, win) / win; integ filter(kernel, 1, sq); % 用 filter 保留因果性不引入前视这里刻意用filter而不是filtfilt。积分这一步在真实系统里是边采边算的用零相位滤波等于用了未来数据仿真指标会虚高移植到实时场景就崩。150 ms 这个窗长要覆盖整个 QRS 波群宽度太短包络会碎成多个峰太长会让相邻心拍粘连取 120–200 ms 之间都能用。3.2 自适应阈值与不应期两个最容易调错的参数阈值策略维护两个滑动均值信号峰均值SPKI和噪声峰均值NPKI。当前阈值由两者插值得到检出真峰就更新SPKI其余采样点更新NPKI。更新系数取 0.125 是原始文献里的经验值它决定了阈值跟随速度。SPKI max(integ(1:2*Fs)); % 用前 2 秒初始化信号峰均值 NPKI mean(integ(1:2*Fs)) / 2; THR NPKI 0.25 * (SPKI - NPKI); % 初始阈值 refract round(0.2 * Fs); % 200 ms 不应期 peaks []; for n 2:length(integ)-1 isLocalMax integ(n) THR integ(n) integ(n-1) integ(n) integ(n1); if isLocalMax (isempty(peaks) || (n - peaks(end)) refract) peaks(end1) n; % 记录一次检出 SPKI 0.125*integ(n) 0.875*SPKI; else NPKI 0.125*integ(n) 0.875*NPKI; end THR NPKI 0.25 * (SPKI - NPKI); % 每个点都重算阈值 end0.25这个系数控制阈值贴近噪声还是贴近信号。调小比如 0.15更灵敏室早这类低幅心拍不容易漏但 T 波高点容易被误判成 R 波调大比如 0.4更保守误检少但漏检多。不应期是硬约束物理上两次心搏不可能在 200 ms 内发生。它的作用是兜底挡住 T 波误检——T 波距离前一个 R 波通常在 250–400 ms200 ms 刚好卡在中间。心率特别快比如运动状态超过 150 bpm时可以把不应期收到 180 ms但再往下就会开始吃掉真实心拍。参数典型值调大后的后果调小后的后果积分窗长150 ms相邻心拍包络粘连单个 QRS 碎裂成多峰阈值系数0.25低幅心拍漏检T 波误判为 R 波不应期200 ms快心率下丢失心拍T 波被当成 R 波SPKI 更新系数0.125阈值跟随慢幅度突变处失准阈值抖动检出位置跳动3.3 从 R 波序列算心率与 HRV 时频域指标拿到 R 波索引后心率是 RR 间期的倒数。这里有一个常见错误直接对 RR 序列做 FFT。RR 间期本身是非均匀采样的时间序列样本点之间的时间间隔并不相等直接 FFT 出来的频谱在频率轴上会整体偏移。rr diff(peaks) / Fs; % RR 间期单位秒 hr 60 ./ rr; % 瞬时心率 bpm hrv.sdnn std(rr) * 1000; % 总体标准差单位 ms hrv.rmssd sqrt(mean(diff(rr).^2)) * 1000; % 相邻差值均方根 hrv.pnn50 sum(abs(diff(rr)) 0.05) / (numel(rr)-1) * 100; % 频域先把 RR 序列按时间插值成均匀序列再估计功率谱 t_rr cumsum(rr); t_rr t_rr - t_rr(1); t_u (0:0.25:t_rr(end)); rr_u interp1(t_rr, rr(2:end), t_u, spline); [pxx, f] pwelch(rr_u - mean(rr_u), hamming(64), 32, 256, 4); % 4 Hz 均匀采样 hrv.lf bandpower(pxx, f, [0.04 0.15], psd); hrv.hf bandpower(pxx, f, [0.15 0.40], psd); hrv.lfhf hrv.lf / hrv.hf;interp1用样条插值把不等间隔的 RR 序列重采样到 4 Hz这是 HRV 频域分析的标准前置。pwelch的第三个参数是重叠点数取窗长一半第四个参数是 FFT 点数决定频率分辨率。bandpower直接按频段积分LF 反映交感和副交感共同活动HF 反映副交感活动两者比值是临床上常用的自主神经平衡指标。4. 特征提取、BP 神经网络分类与 MATLAB 仿真界面搭建检测阶段解决的是“心拍在哪”分类阶段要回答“这个心拍属于哪一类”。MIT-BIH 的标注里有正常、室早、房早、左束支阻滞、右束支阻滞等类型做五分类是对课程设计比较合适的规模。4.1 小波变换提取 QRS 形态特征只用幅值和时间间隔做特征很难区分形态相近的室早和束支阻滞。小波变换的好处是能在不同尺度上分别看细节低频尺度反映 QRS 的整体轮廓高频尺度反映切迹和毛刺。wname db4; level 4; seg ecg_clean(idx-100 : idx150); % 以 R 波为中心取 250 点窗口 [c, l] wavedec(seg, level, wname); % 4 层小波分解 d4 wrcoef(d, c, l, wname, 4); % 重构第 4 层细节系数对应低频形态 feat [ ... max(seg), min(seg), ... % 峰值与谷值反映 QRS 幅度 rms(d4), ... % 第 4 层细节能量 trapz(abs(seg)), ... % 波形面积对宽窄敏感 kurtosis(seg)]; % 峰度反映波形尖锐程度db4的小波函数形状窄而尖和 QRS 波群形态接近是小波检测常用的选择。窗口取 R 波前 100 点后 150 点是因为 Q 波在 R 波之前这个窗口能完整包住 QRS 加一部分 ST 段。特征向量维度控制在 10–20 维比较合适维度过高在小样本上容易过拟合。4.2 用 BP 神经网络做心拍五分类MATLAB 的深度学习工具箱里patternnet专做模式分类它把输出层设成 softmax、损失函数用交叉熵比自己拼网络省事。输入是特征矩阵按列排也就是每列一个样本。% X: 特征维度 × 样本数Y: one-hot 标签 类别数 × 样本数 net patternnet([12 8]); % 两个隐藏层节点数 12 和 8 net.trainFcn trainlm; % Levenberg-Marquardt小样本收敛快 net.trainParam.epochs 300; % 最大迭代轮数 net.trainParam.goal 1e-4; % 目标均方误差 net.trainParam.max_fail 12; % 验证集连续不降就早停 net.divideParam.trainRatio 0.70; net.divideParam.valRatio 0.15; net.divideParam.testRatio 0.15; [net, tr] train(net, X, Y); y_pred net(X_test); % 输出为各类概率 [~, cls] max(y_pred, [], 1); % 取概率最大的类别 plotperform(tr); % 可视化 bp 神经网络拟合曲线trainlm在几百到几千样本量下通常几十轮就收敛max_fail是早停保护验证集误差连续 12 轮不下降就停防止在训练集上过拟合。plotperform画出的三条曲线就是常说的拟合曲线训练曲线一直降而验证曲线抬头说明网络开始记样本了得减少隐藏层节点或者增加样本。训练参数取值影响隐藏层结构[12 8]层数过多在千级样本上必过拟合学习算法trainlm内存占用高样本过万时改用 trainscgmax_fail12调小收敛快但可能停在次优解数据划分比例7:1.5:1.5测试集太少评估方差大如果想让参数搜索省点力气可以用 MATLAB 优化工具箱里的bayesopt把隐藏层节点数和max_fail一起当超参搜目标设成验证集准确率通常几十次迭代就能明显好过手调。4.3 App Designer 仿真界面与信号发生器仿真输入界面部分用 App Designer 搭左边放坐标区显示原始波形和去噪波形右边放 R 波列表和 HRV 指标顶部放导入按钮和几个旋钮。一个常被忽略的细节是调试界面逻辑时不该每次都去读数据库文件最好内置一个合成心电信号发生器参数可控出问题一眼能看出是哪一环。function sig synth_ecg(Fs, dur, hr_bpm) %SYNTH_ECG 合成带基线漂移和白噪声的心电信号用于界面仿真 t (0:1/Fs:dur-1/Fs); rr 60 / hr_bpm; sig 0.15 * sin(2*pi*0.3*t); % 呼吸引起的基线漂移 for k 1:round(dur/rr) c round((k-1)*rr*Fs); idx c (1:round(0.08*Fs)); idx idx(idx numel(t)); off 0.02*Fs; sig(idx) sig(idx) 1.2 * exp(-((1:numel(idx))-off).^2 / (2*(0.008*Fs)^2)); end sig sig 0.02 * randn(size(t)); % 白噪声 end这个合成函数把三个关键量都做成了入参采样率、时长、心率。0.008*Fs控制 R 波宽度0.02*Fs控制 R 波在窗口内的位置。改心率就能测快慢心率下的检测表现加大randn系数就能测低信噪比下的鲁棒性。界面上的“仿真”按钮调它比反复读真实记录做对比测试快得多。注意App Designer 的回调里不要直接跑耗时几秒的循环界面会卡成假死。把处理逻辑封成独立函数用drawnow在关键节点刷新或者拆成定时器分批执行。5. 检测性能评估与仿真排错把 QRS 检出率调到 99% 以上算法写完了不等于系统成立得拿指标说话。QRS 检测的标准评价有三个量敏感度 Se、阳性预测值 PPV 和检测错误率 DER。Se 衡量真值里被检出的比例PPV 衡量检出结果里有多少是真的DER 是两者的综合误差。function [Se, PPV, DER] eval_detect(ref, det, Fs) %EVAL_DETECT 以 150 ms 匹配窗评估 QRS 检测性能 tol round(0.150 * Fs); TP 0; for i 1:numel(ref) % 每个真值找最近检出 if any(abs(det - ref(i)) tol) TP TP 1; end end FP numel(det) - TP; FN numel(ref) - TP; Se TP / numel(ref); PPV TP / numel(det); DER (FP FN) / numel(ref); endtol取 150 ms 是通用约定因为不同标注者对同一个 R 波的位置判断本身就有几十毫秒的差异。这个函数只能在整段记录上跑不能只喂前 10 秒——短段上的 PPV 方差非常大前面漏一个后面多一个数字就难看。仿真跑出来结果不对时我一般按固定顺序排查而不是乱改参数。现象可能原因验证方式检出位置整体滞后用了 filter 而非 filtfilt对比滤波前后同一 R 波索引幅度小的心拍全部漏检阈值系数偏大或 SPKI 初始化过高打印 THR 随时间的曲线一个心拍报出两三个位置积分窗太短或漏了不应期判断检查相邻检出点间隔是否小于 200 ms结果发散、数值爆炸滤波器阶数过高或未去均值看滤波输出是否超出原信号量级换一份记录指标骤降采样率变了但时间参数没重算打印 Fs 和 win、refract 的取值最后一步是批量化验证。把 MIT-BIH 里十到二十条记录一次性跑完统计总体 Se 和 PPV比在单条记录上抠参数有意义得多。写一个循环脚本把每条记录的结果和实参一起存成 mat 文件跑完直接汇总如果总体 Se 低于 99%重点看是不是某几条低信噪比记录拖了后腿而不是去动已经在多数记录上表现正常的阈值。这种“先看分布再调参数”的顺序能把排查时间从几小时压到十几分钟。本文还有配套的精品资源点击获取
简介这份资源是面向生物医学工程、电子信息及相关专业学生的心电信号分析课程设计参考资料以PDF形式呈现课题二“基于MATLAB平台的心电信号分析系统设计及仿真”的完整方案适合正在做数字信号处理课程设计或信号与系统综合实验的读者参考。资源包共1个PDF文件体积约1019KB内容围绕MATLAB程序设计展开。文档完整梳理了从MIT-BIH数据库读取数字心电信号、还原实际波形、对非均匀采样数据做线性插值到依据心电信号频域特点设计低通与高通滤波器、滤除基线漂移与工频干扰再到对滤波前后信号做频谱分析与比较的整套流程并给出Simulink动态建模与仿真的设计思路附录部分还整理了心电信号读取原理、滤波器指标选取与课程设计报告撰写要求。目前已有357人学习可为选题开题、程序编写与结果分析提供一份可对照的工程实例。1. 心电信号分析系统为什么值得用 MATLAB 从零搭一遍做信号处理方向课题的人往往先找现成工具箱跑一遍图出来挺漂亮可一旦被追问“这个 R 波凭什么标在这个位置”“把滤波器阶数从 3 改成 5 结果会怎样”就答不上来了。基于 MATLAB 平台把心电信号分析系统从数据读取、噪声去除、QRS 检测一路搭到分类识别和仿真界面暴露出来的问题会非常具体采样率不同时间参数要重算滤波器相位没处理 T 波就会整体平移自适应阈值不收敛得去翻积分窗而不是改检测逻辑。这套系统设计的输入是一条原始心电记录输出是干净波形、R 波位置序列、心率与心率变异性指标外加一个点一下就能跑的仿真界面。它覆盖了生理信号处理里最典型的那几类噪声和评价方法适合生物医学工程、电子信息方向的课程设计或毕业设计也适合想练一练 MATLAB 工程化写法的开发者。2. 心电信号读入、噪声建模与 MATLAB 工程骨架设计心电信号分析系统的第一道坎不是算法而是数据从哪来、长什么样。公开的心电数据里MIT-BIH 心律失常数据库是课程设计用得最多的一份它每条记录由三个文件组成.hea描述采样率、导联和增益.dat存二进制采样值.atr存专家标注的心拍位置和类型。搞清楚这三个文件的对应关系后面所有评价指标才有参照物。2.1 MIT-BIH 记录在 MATLAB 里的读取与通道选择常见做法是装 WFDB Toolbox for MATLAB用它提供的rdsamp和rdann两个函数把信号和标注读进来。这两行代码看着简单但通道选错、时长截错是新手最容易踩的坑。% 读取 MIT-BIH 100 号记录的前 10 秒 % 该记录采样率 360 Hz10 秒对应 3600 个采样点 [signal, Fs, tm] rdsamp(mitdb/100, 1:3600); % 读取同一记录的专家标注anntype 里 N 表示正常心拍V 表示室早 [ann, anntype] rdann(mitdb/100, atr);rdsamp返回三个量signal是 N×C 的矩阵C 是导联数MIT-BIH 多数记录是两导联第一导联通常是 MLIIFs是采样率100 号记录是 360 Hztm是时间列单位秒。rdann返回的ann是以采样点为单位的 R 波位置索引anntype是与之等长的字符数组。这里有个细节要注意rdann读出来的是整条记录的标注如果只取了前 3600 点做分析必须自己按索引范围裁剪否则参考位置和检测结果对不上。2.1.1 采样率对后续所有时间参数的影响采样率决定了每一个以秒为单位的时间参数要换算成多少个点。360 Hz 下 150 ms 积分窗是 54 个点1000 Hz 下就是 150 个点。很多人在自己录的数据上跑得好好的代码换到 MIT-BIH 上就漏检一堆原因就是积分窗、不应期这些参数还按原来 1000 Hz 写的。噪声/干扰主要频段常用处理关键参数基线漂移0.05–1 Hz高通滤波或中值滤波截止 0.5 Hz中值窗 0.6 s工频干扰50 Hz ±1 HzIIR 陷波Q 值 30–40肌电噪声20–100 Hz低通滤波截止 40 Hz电极接触噪声突变阶跃中值滤波或形态滤波窗长 0.2 s2.2 工频干扰、基线漂移与肌电噪声的滤波链处理顺序上我一般按“先去基线漂移再去工频最后压高频肌电”来排。原因是基线漂移幅度能达到毫伏级先不处理会让后面所有阈值判断全部失真。用 Butterworth 带通一次性把 0.5 Hz 以下和 40 Hz 以上切掉再用 IIR 陷波打掉 50 Hz这是最省事也最稳的组合。% 0.5-40 Hz 带通切掉基线漂移和高频肌电 [b_bp, a_bp] butter(3, [0.5 40]/(Fs/2), bandpass); ecg_bp filtfilt(b_bp, a_bp, ecg_raw); % 50 Hz 工频陷波Q 值决定陷波宽度 wo 50/(Fs/2); bw wo/35; [b_nc, a_nc] iirnotch(wo, bw); ecg_clean filtfilt(b_nc, a_nc, ecg_bp);butter的第一个参数是阶数3 阶在绝大多数场景够用阶数往上加到 5、6 会有更陡的过渡带但相位非线性也更严重。filtfilt做的是前向加后向两次滤波等效零相位代价是有效阶数翻倍、边界会有短暂瞬态。如果这里图省事用filter输出的 QRS 波在时间轴上会整体滞后十几到几十毫秒后面跟专家标注比对时就会系统性错位。iirnotch的第二个参数是归一化带宽写成wo/35等价于 Q 值 35。Q 越大陷波越窄对 50 Hz 附近有用信号损伤越小但对工频频率漂移的容忍度也越低。国内电网频率相对稳定Q 取 30–40 是合理区间。2.3 工程目录结构与函数接口约定系统设计到后半程代码会膨胀到十几个文件这时候目录怎么分比算法本身更影响调试效率。我习惯按“数据—处理—检测—特征—界面”五层拆目录data/放原始记录和缓存src/preprocess/放滤波函数src/detect/放 QRS 相关src/features/放特征提取和分类src/gui/放界面回调results/落盘中间结果和评价指标。每个函数只干一件事输入输出用结构体串起来方便单独拿一段信号复现问题。function out ecg_preprocess(raw, Fs, cfg) %ECG_PREPROCESS 心电信号预处理带通 陷波 % raw 原始信号 N×1Fs 采样率cfg 可覆盖默认参数 % 返回 out.clean 去噪信号、out.fs 采样率、out.cfg 实际配置 arguments raw (:,1) double Fs (1,1) double {mustBePositive} cfg struct struct(bp,[0.5 40],notch,50,order,3,Q,35) end [b, a] butter(cfg.order, cfg.bp/(Fs/2), bandpass); x filtfilt(b, a, raw); wo cfg.notch/(Fs/2); [bn, an] iirnotch(wo, wo/cfg.Q); out.clean filtfilt(bn, an, x); out.fs Fs; out.cfg cfg; end用arguments块声明参数好处是配置项集中在一处界面传参时不用记位置还能顺便做类型和正数校验。这里把cfg设成带默认值的结构体界面上想改截止频率只需要传struct(bp,[1 35])覆盖一项其余保持默认。提示把每次运行实际生效的cfg一起存进结果文件过两周回头复现某个图的时候能省掉大量猜参数的时间。3. QRS 波检测、心率计算与阈值参数整定QRS 检测是整个心电信号分析系统的枢纽。R 波位置一旦错心率、RR 间期、HRV 指标、心拍分类的切窗位置全都会跟着错。Pan-Tompkins 之所以成为事实标准是因为它用一套纯时域的因果处理链能在嵌入式环境里实时跑同时自适应阈值让它对幅度变化有不错的鲁棒性。3.1 Pan-Tompkins 五步法在 MATLAB 中的逐步实现五步依次是5–15 Hz 带通、微分、平方、滑动窗积分、自适应双阈值。带通这一步是为了让 QRS 的 10 Hz 左右主频段凸出来同时把 P 波和 T 波压下去。微分把波形的斜率信息放大平方让所有值非负并进一步拉开 QRS 与其他波的差距积分把能量摊平成一个包络。% Pan-Tompkins 前四级处理 [b, a] butter(3, [5 15]/(Fs/2), bandpass); x filtfilt(b, a, ecg_clean); % 第一步带通 d [0; diff(x)]; % 第二步一阶差分近似微分 sq d.^2; % 第三步平方非负且放大能量差 win max(1, round(0.150 * Fs)); % 第四步150 ms 积分窗 kernel ones(1, win) / win; integ filter(kernel, 1, sq); % 用 filter 保留因果性不引入前视这里刻意用filter而不是filtfilt。积分这一步在真实系统里是边采边算的用零相位滤波等于用了未来数据仿真指标会虚高移植到实时场景就崩。150 ms 这个窗长要覆盖整个 QRS 波群宽度太短包络会碎成多个峰太长会让相邻心拍粘连取 120–200 ms 之间都能用。3.2 自适应阈值与不应期两个最容易调错的参数阈值策略维护两个滑动均值信号峰均值SPKI和噪声峰均值NPKI。当前阈值由两者插值得到检出真峰就更新SPKI其余采样点更新NPKI。更新系数取 0.125 是原始文献里的经验值它决定了阈值跟随速度。SPKI max(integ(1:2*Fs)); % 用前 2 秒初始化信号峰均值 NPKI mean(integ(1:2*Fs)) / 2; THR NPKI 0.25 * (SPKI - NPKI); % 初始阈值 refract round(0.2 * Fs); % 200 ms 不应期 peaks []; for n 2:length(integ)-1 isLocalMax integ(n) THR integ(n) integ(n-1) integ(n) integ(n1); if isLocalMax (isempty(peaks) || (n - peaks(end)) refract) peaks(end1) n; % 记录一次检出 SPKI 0.125*integ(n) 0.875*SPKI; else NPKI 0.125*integ(n) 0.875*NPKI; end THR NPKI 0.25 * (SPKI - NPKI); % 每个点都重算阈值 end0.25这个系数控制阈值贴近噪声还是贴近信号。调小比如 0.15更灵敏室早这类低幅心拍不容易漏但 T 波高点容易被误判成 R 波调大比如 0.4更保守误检少但漏检多。不应期是硬约束物理上两次心搏不可能在 200 ms 内发生。它的作用是兜底挡住 T 波误检——T 波距离前一个 R 波通常在 250–400 ms200 ms 刚好卡在中间。心率特别快比如运动状态超过 150 bpm时可以把不应期收到 180 ms但再往下就会开始吃掉真实心拍。参数典型值调大后的后果调小后的后果积分窗长150 ms相邻心拍包络粘连单个 QRS 碎裂成多峰阈值系数0.25低幅心拍漏检T 波误判为 R 波不应期200 ms快心率下丢失心拍T 波被当成 R 波SPKI 更新系数0.125阈值跟随慢幅度突变处失准阈值抖动检出位置跳动3.3 从 R 波序列算心率与 HRV 时频域指标拿到 R 波索引后心率是 RR 间期的倒数。这里有一个常见错误直接对 RR 序列做 FFT。RR 间期本身是非均匀采样的时间序列样本点之间的时间间隔并不相等直接 FFT 出来的频谱在频率轴上会整体偏移。rr diff(peaks) / Fs; % RR 间期单位秒 hr 60 ./ rr; % 瞬时心率 bpm hrv.sdnn std(rr) * 1000; % 总体标准差单位 ms hrv.rmssd sqrt(mean(diff(rr).^2)) * 1000; % 相邻差值均方根 hrv.pnn50 sum(abs(diff(rr)) 0.05) / (numel(rr)-1) * 100; % 频域先把 RR 序列按时间插值成均匀序列再估计功率谱 t_rr cumsum(rr); t_rr t_rr - t_rr(1); t_u (0:0.25:t_rr(end)); rr_u interp1(t_rr, rr(2:end), t_u, spline); [pxx, f] pwelch(rr_u - mean(rr_u), hamming(64), 32, 256, 4); % 4 Hz 均匀采样 hrv.lf bandpower(pxx, f, [0.04 0.15], psd); hrv.hf bandpower(pxx, f, [0.15 0.40], psd); hrv.lfhf hrv.lf / hrv.hf;interp1用样条插值把不等间隔的 RR 序列重采样到 4 Hz这是 HRV 频域分析的标准前置。pwelch的第三个参数是重叠点数取窗长一半第四个参数是 FFT 点数决定频率分辨率。bandpower直接按频段积分LF 反映交感和副交感共同活动HF 反映副交感活动两者比值是临床上常用的自主神经平衡指标。4. 特征提取、BP 神经网络分类与 MATLAB 仿真界面搭建检测阶段解决的是“心拍在哪”分类阶段要回答“这个心拍属于哪一类”。MIT-BIH 的标注里有正常、室早、房早、左束支阻滞、右束支阻滞等类型做五分类是对课程设计比较合适的规模。4.1 小波变换提取 QRS 形态特征只用幅值和时间间隔做特征很难区分形态相近的室早和束支阻滞。小波变换的好处是能在不同尺度上分别看细节低频尺度反映 QRS 的整体轮廓高频尺度反映切迹和毛刺。wname db4; level 4; seg ecg_clean(idx-100 : idx150); % 以 R 波为中心取 250 点窗口 [c, l] wavedec(seg, level, wname); % 4 层小波分解 d4 wrcoef(d, c, l, wname, 4); % 重构第 4 层细节系数对应低频形态 feat [ ... max(seg), min(seg), ... % 峰值与谷值反映 QRS 幅度 rms(d4), ... % 第 4 层细节能量 trapz(abs(seg)), ... % 波形面积对宽窄敏感 kurtosis(seg)]; % 峰度反映波形尖锐程度db4的小波函数形状窄而尖和 QRS 波群形态接近是小波检测常用的选择。窗口取 R 波前 100 点后 150 点是因为 Q 波在 R 波之前这个窗口能完整包住 QRS 加一部分 ST 段。特征向量维度控制在 10–20 维比较合适维度过高在小样本上容易过拟合。4.2 用 BP 神经网络做心拍五分类MATLAB 的深度学习工具箱里patternnet专做模式分类它把输出层设成 softmax、损失函数用交叉熵比自己拼网络省事。输入是特征矩阵按列排也就是每列一个样本。% X: 特征维度 × 样本数Y: one-hot 标签 类别数 × 样本数 net patternnet([12 8]); % 两个隐藏层节点数 12 和 8 net.trainFcn trainlm; % Levenberg-Marquardt小样本收敛快 net.trainParam.epochs 300; % 最大迭代轮数 net.trainParam.goal 1e-4; % 目标均方误差 net.trainParam.max_fail 12; % 验证集连续不降就早停 net.divideParam.trainRatio 0.70; net.divideParam.valRatio 0.15; net.divideParam.testRatio 0.15; [net, tr] train(net, X, Y); y_pred net(X_test); % 输出为各类概率 [~, cls] max(y_pred, [], 1); % 取概率最大的类别 plotperform(tr); % 可视化 bp 神经网络拟合曲线trainlm在几百到几千样本量下通常几十轮就收敛max_fail是早停保护验证集误差连续 12 轮不下降就停防止在训练集上过拟合。plotperform画出的三条曲线就是常说的拟合曲线训练曲线一直降而验证曲线抬头说明网络开始记样本了得减少隐藏层节点或者增加样本。训练参数取值影响隐藏层结构[12 8]层数过多在千级样本上必过拟合学习算法trainlm内存占用高样本过万时改用 trainscgmax_fail12调小收敛快但可能停在次优解数据划分比例7:1.5:1.5测试集太少评估方差大如果想让参数搜索省点力气可以用 MATLAB 优化工具箱里的bayesopt把隐藏层节点数和max_fail一起当超参搜目标设成验证集准确率通常几十次迭代就能明显好过手调。4.3 App Designer 仿真界面与信号发生器仿真输入界面部分用 App Designer 搭左边放坐标区显示原始波形和去噪波形右边放 R 波列表和 HRV 指标顶部放导入按钮和几个旋钮。一个常被忽略的细节是调试界面逻辑时不该每次都去读数据库文件最好内置一个合成心电信号发生器参数可控出问题一眼能看出是哪一环。function sig synth_ecg(Fs, dur, hr_bpm) %SYNTH_ECG 合成带基线漂移和白噪声的心电信号用于界面仿真 t (0:1/Fs:dur-1/Fs); rr 60 / hr_bpm; sig 0.15 * sin(2*pi*0.3*t); % 呼吸引起的基线漂移 for k 1:round(dur/rr) c round((k-1)*rr*Fs); idx c (1:round(0.08*Fs)); idx idx(idx numel(t)); off 0.02*Fs; sig(idx) sig(idx) 1.2 * exp(-((1:numel(idx))-off).^2 / (2*(0.008*Fs)^2)); end sig sig 0.02 * randn(size(t)); % 白噪声 end这个合成函数把三个关键量都做成了入参采样率、时长、心率。0.008*Fs控制 R 波宽度0.02*Fs控制 R 波在窗口内的位置。改心率就能测快慢心率下的检测表现加大randn系数就能测低信噪比下的鲁棒性。界面上的“仿真”按钮调它比反复读真实记录做对比测试快得多。注意App Designer 的回调里不要直接跑耗时几秒的循环界面会卡成假死。把处理逻辑封成独立函数用drawnow在关键节点刷新或者拆成定时器分批执行。5. 检测性能评估与仿真排错把 QRS 检出率调到 99% 以上算法写完了不等于系统成立得拿指标说话。QRS 检测的标准评价有三个量敏感度 Se、阳性预测值 PPV 和检测错误率 DER。Se 衡量真值里被检出的比例PPV 衡量检出结果里有多少是真的DER 是两者的综合误差。function [Se, PPV, DER] eval_detect(ref, det, Fs) %EVAL_DETECT 以 150 ms 匹配窗评估 QRS 检测性能 tol round(0.150 * Fs); TP 0; for i 1:numel(ref) % 每个真值找最近检出 if any(abs(det - ref(i)) tol) TP TP 1; end end FP numel(det) - TP; FN numel(ref) - TP; Se TP / numel(ref); PPV TP / numel(det); DER (FP FN) / numel(ref); endtol取 150 ms 是通用约定因为不同标注者对同一个 R 波的位置判断本身就有几十毫秒的差异。这个函数只能在整段记录上跑不能只喂前 10 秒——短段上的 PPV 方差非常大前面漏一个后面多一个数字就难看。仿真跑出来结果不对时我一般按固定顺序排查而不是乱改参数。现象可能原因验证方式检出位置整体滞后用了 filter 而非 filtfilt对比滤波前后同一 R 波索引幅度小的心拍全部漏检阈值系数偏大或 SPKI 初始化过高打印 THR 随时间的曲线一个心拍报出两三个位置积分窗太短或漏了不应期判断检查相邻检出点间隔是否小于 200 ms结果发散、数值爆炸滤波器阶数过高或未去均值看滤波输出是否超出原信号量级换一份记录指标骤降采样率变了但时间参数没重算打印 Fs 和 win、refract 的取值最后一步是批量化验证。把 MIT-BIH 里十到二十条记录一次性跑完统计总体 Se 和 PPV比在单条记录上抠参数有意义得多。写一个循环脚本把每条记录的结果和实参一起存成 mat 文件跑完直接汇总如果总体 Se 低于 99%重点看是不是某几条低信噪比记录拖了后腿而不是去动已经在多数记录上表现正常的阈值。这种“先看分布再调参数”的顺序能把排查时间从几小时压到十几分钟。本文还有配套的精品资源点击获取
