
简介基于Haar小波变换的心电信号去噪Matlab源码面向信号处理、生物医学工程方向的本科生与研究生适合课设、毕设及教研实验使用。资源聚焦心电信号中工频干扰与随机噪声的抑制问题通过小波分解、阈值处理与重构完成去噪流程并配套完整可运行的代码。zip压缩包共14个文件包含m主程序、mat与dat测试数据、txt说明文档及多张运行结果图整体仅3.09MB轻量易用。已有402人学习下载便于快速上手。除主体代码外资源还提供演示图片和仿真咨询入口读者可直观对照去噪前后的波形效果并能基于示例数据修改参数、迁移到其他生物电信号处理场景可作为算法对比实验或论文复现的实用参考。1. 心电去噪的第一个选择为什么是Haar小波用 Matlab 跑 MIT-BIH 的 100.dat 心电信号第一件要决定的事不是滤波器阶数而是“你到底要留下什么”。这路信号里 QRS 波群的主能量集中在 10 Hz 附近而肌电噪声、工频干扰在时域上完全贴着 R 波叠加用普通 FIR 带通滤波会在波峰处产生振铃。我在这套源码里最看重的就是它选了 Haar 小波——这是唯一同时具备正交、紧支撑和对称性的小波基分解过程只有加减和移位没有任何浮点乘法软硬件成本都低。对本科生的课程设计、研究生的预处理环节甚至嵌入式实时心电监测它都是一个能讲清楚原理、又能直接出结果的选择。我会从噪声构成讲起把 process.m 的每一段逻辑拆开最后给出验证指标和批量处理写法。2. ECG噪声成分与Haar小波分解原理2.1 三类主要噪声及其频段分布心电信号的有效频带是 0.05~100 Hz其中 QRS 波群能量集中在 5~20 HzP 波和 T 波则在 5 Hz 以下。实际采集时叠加的噪声按频段可以分成三类基线漂移来自呼吸、肢体运动和电极极化频率通常低于 0.5 Hz工频干扰在电网频率 50 Hz部分地区 60 Hz处表现为正弦纹波肌电干扰来自肌肉收缩频谱从几赫兹一直扩展到 2000 Hz与心电信号频带重叠最大处理起来最麻烦。噪声类型频段典型来源对诊断形态的影响基线漂移0.05~0.5 Hz呼吸、电极活动ST 段整体偏移工频干扰50/60 Hz ± 1 Hz电网耦合、接地不良波形叠加正弦纹波肌电干扰5~2000 Hz肌肉收缩R 波变钝、出现毛刺如果用巴特沃斯带通滤波器去压制工频干扰在 50 Hz 附近的陷波会把 QRS 波的高频拐角削掉R 波幅度直接损失 10% 以上。小波变换的优势在于它把信号分解成不同尺度下的逼近和细节可以在不触碰低频逼近系数的前提下只对特定尺度的细节系数做处理这正是心电信号去噪里选择小波的核心理由。2.2 Haar小波的尺度函数与小波函数Haar 小波是 1909 年提出的正交小波尺度函数在 [0,1) 区间取值为 1小波函数在 [0,0.5) 为 1、在 [0.5,1) 为 -1。对离散序列做 Haar 分解就是不断把相邻两个采样点做平均和差分设信号 x [x1, x2, x3, x4]一层分解后的低频系数是 (x1x2)/√2 和 (x3x4)/√2高频系数是 (x1-x2)/√2 和 (x3-x4)/√2。这里面的 √2 是归一化因子用来保证变换前后能量守恒去掉它虽然不影响滤波趋势的理解但重构时幅值会整体偏大。Mallat 算法就是把这种两两配对的多分辨分析递归进行每一层的低频系数再作为下一层的输入继续分解第 k 层的细节系数捕获的是 2^k 个采样点尺度上的变化。对 fs360 Hz 的 100.dat 信号来说一层分解对应约 0.0028 秒的时间分辨率足够分辨 QRS 波群的陡峭上升沿。2.3 多分辨分析噪声落在哪一层对 fs360 Hz 的信号做 5 层 Haar 分解各层频带大致如下第 1 层细节 cD1 覆盖 90~180 Hz第 2 层细节 cD2 覆盖 45~90 Hz第 3 层 cD3 覆盖 22.5~45 Hz第 4 层 cD4 覆盖 11.25~22.5 Hz第 5 层逼近 cA5 覆盖 0~5.625 Hz。这样一个频率映射关系可以直接对照 2.1 节的噪声表去判断工频干扰主要落在 cD2 附近肌电干扰散落在 cD1 和 cD2基线漂移则全部沉在 cA5 里。去噪策略因此变得非常清晰保留 cA5 或在 cA5 上做轻微的高通处理避免破坏 ST 段和 T 波形态对 cD2 之前的高频细节层做阈值收缩把幅度小的噪声系数压掉cD3 和 cD4 与 QRS 波群的主频带重合这一层的阈值要放得宽一些免得把 R 波削平。这套“分尺度看噪声”的思路是理解 process.m 后续代码的钥匙所谓小波去噪本质是信噪比在尺度域分离后对噪声主导的系数做幅度处罚。3. process.m源码拆解212格式读取、阈值处理与重构3.1 数据源与MIT-BIH的212格式包里的数据有两种形态ECGdata.mat 是已经准备好的 Mat 变量直接用 load 读入100.dat 和 ECG1.dat 是 MIT-BIH 数据库的原始记录。MIT-BIH 的 .dat 文件是 212 格式两个 12 位采样值打包在 3 个字节里第一位采样取第 1、2 字节的低 12 位第二位采样取第 2 字节的高 4 位加上第 3 字节的全部 8 位。如果直接 fread 成 int16 再截断得到的数据顺序是乱的R 波位置会完全对不上。我一般在 process.m 里先做数据源判断文件存在就优先加载 .mat否则走 212 解包路径if exist(ECGdata.mat, file) load(ECGdata.mat, ecg); % 变量名以实际文件为准 fs 360; % MIT-BIH 标准采样率 else fid fopen(100.dat, rb); raw fread(fid, Inf, uint8); fclose(fid); nPair floor(length(raw) / 3); ecg zeros(nPair * 2, 1); for i 1:nPair b raw((i-1)*3 1 : i*3); s1 b(1) bitshift(bitand(b(2), 15), 8); % 低12位 s2 bitshift(b(2), -4) bitshift(b(3), 4);% 高12位 if s1 2048, s1 s1 - 4096; end % 符号扩展 if s2 2048, s2 s2 - 4096; end ecg(2*i - 1) s1; ecg(2*i) s2; end ecg ecg / 200; % 幅值标定单位mV end这段代码里的 bitshift(bitand(b(2), 15), 8) 是把第二字节的低四位提出来左移八位后和第一字节拼出第一采样bitshift(b(2), -4) 则是把第二字节的高四位挪到低位再与第三字节拼出第二采样。阈值 2048 对应 12 位补码的符号位超过即视为负数除以 200 是因为 MIT-BIH 数据规定的增益是 200 ADC 单位每 mV。这套写法虽然循环逐点处理略慢但胜在直观跑 10 分钟的数据也没有明显卡顿之后所有去噪处理都基于这个 ecg 向量。3.2 核心去噪wavedec wthresh waverecprocess.m 的去噪主线是三个 Matlab Wavelet Toolbox 函数。wavedec 完成多级分解wthresh 做阈值收缩waverec 重建信号。源码中的核心段大致如下wavelet haar; level 5; [C, L] wavedec(ecg, level, wavelet); % 多级小波分解 % Donoho通用阈值基于绝对中位差估计噪声水平 sigma median(abs(C)) / 0.6745; thr sigma * sqrt(2 * log(length(ecg))); % 细节系数全部软阈值收缩逼近系数不动 Ctmp C; for k 1:level idx sum(L(1:k)) 1 : sum(L(1:k1)); Ctmp(idx) wthresh(Ctmp(idx), s, thr); end % 重构 ecg_den waverec(Ctmp, L, wavelet);wavedec 的返回参数 C 是自顶向下的系数拼接向量低频逼近系数在最前面然后是各层细节系数L 向量记录了每一段的长度所以取第 k 层细节的索引范围要看 L(1:k) 的累加和。噪声标准差先对全体系数取中位数再除以 0.6745这是小波去噪里最常用的鲁棒估计因为噪声系数占多数中位数不受大系数影响。阈值公式里的 sqrt(2*log(N)) 来自 Donoho 极小极大风险推导N 取整个信号长度即可。wthresh 的 soft 模式对所有超过阈值的系数做整体收缩数学形式是 sign(x)·(|x|-thr)结果信号连续但幅值略被压缩hard 模式是超过阈值留原值、低于阈值归零波形更锐利但阈值点不连续。心电去噪若要保留 R 波幅值soft 之后的幅值损失约等于 thrhard 的伪 Gibbs 现象则在重构后表现为局部毛刺二者各有代价具体怎么选在第 4 章展开。3.3 波形绘图与残差观察我把绘图部分放在同一个脚本里三行 subplot 分别显示原始、去噪、残差t (0:length(ecg)-1) / fs; figure(Color, w); subplot(3,1,1); plot(t, ecg); ylim([-1.5 2]); title(原始ECG信号); subplot(3,1,2); plot(t, ecg_den); ylim([-1.5 2]); title(Haar小波去噪结果); subplot(3,1,3); plot(t, ecg - ecg_den); title(去噪残差); xlabel(时间/s);残差图的价值在于如果残差里仍能看出周期性的 QRS 波形轮廓说明阈值设得偏大把心拍主体一起削掉了如果残差是被压扁的白噪声形状且幅度都在阈值附近说明参数基本合理。另外建议把 ylim 固定在原始数据相同的幅值范围这样对比时不会因为坐标轴自适应而误判去噪后的幅值大小。包里 7 张运行结果图对应的是不同数据段和不同 level 参数下的输出换参数重跑时可以直接拿这些图做参照。4. 阈值策略、分解层数与边界延拓调优4.1 软硬阈值与折中方案阈值方式表达式优点心电场景的风险硬阈值y x·I(xthr)软阈值y sign(x)·max(x-thr,0)折中阈值y sign(x)·max(x-a·thr,0)多数 ECG 去噪论文用软阈值因为后续要计算 SNR 和 PRD连续信号在误差指标上更稳定。但如果你的目标是提取 R 波做心率计算我更推荐折中阈值a 取 0.5~0.8这样比软阈值少损失一半幅度又保留了硬阈值对噪声的抑制能力。参数 a 用网格搜索在 [0.3, 0.9] 之间每 0.1 扫一次取 PRD 最小的组合即可这个做法在源码注释里也提到了。4.2 分解层数如何对应采样率分解层数决定最低频细节层覆盖多宽的频带一般原则是让工频干扰落在一个细节层而不是逼近层。对 fs360 的信号4 层分解时 cD2 中心约 45~90 Hz工频 50 Hz 正好落在里面所以 4 层就够工频抑制5 层分解把最低频带降得更细适合基线漂移明显的记录。超过 6 层之后有一个容易忽视的问题Haar 滤波器在深度分解时边界效应加剧逼近系数里会掺入边界反射重构后的 QRS 波起点和终点会出现小幅 overshoot。包里的 100.dat 信号质量尚可处理时我建议从 level5 起步如果切换数据发现 ST 段扭曲优先减少一层而不是继续增加层数。这个现象也解释了为什么原始代码里没有超过 6 层的预置。4.3 延拓模式与分段处理wavedec 默认使用 dwtmode(sym) 对称延拓信号两端按镜像方式补充对具有对称形态的 QRS 波最友好。如果修改成 dwtmode(zpd) 零填充去噪后的首尾 10~20 个点会明显下垂改成 dwtmode(per) 周期延拓则要求信号长度能整除 2^level否则报错或自动重新填充。提示dwtmode 是全局设置修改一次会对后续所有 wavedec 调用生效跑完记得改回 sym。在 360 Hz 下一个典型心拍不足 200 个采样点如果按整段数据处理边界问题只影响一个心拍左右。但当数据长到数分钟甚至数小时内存占用和重构时间会明显增长我一般先用 buffer 把信号切成每段 4096 点4096 可被 2^5 整除满足 per 模式约束逐段去噪后拼接。拼接处前后各留 64 点重叠用线性交叉淡入淡出消除接缝segLen 4096; hop segLen - 64; nSeg floor((length(ecg) - segLen) / hop) 1; ecg_out zeros(size(ecg)); for k 1:nSeg idx (k-1)*hop 1 : (k-1)*hop segLen; seg ecg(idx); seg_den seg_denoise(seg, haar, level, thr); % 独立去噪 if k 1 prevEnd (k-1)*hop - 63 : (k-1)*hop; w linspace(0, 1, 64); % 权重递增 ecg_out(prevEnd) ecg_out(prevEnd) .* (1-w) seg_den(1:64) .* w; end ecg_out(idx) seg_den; end重叠区的权重 w 用 linspace(0,1,64) 生成前一帧权重从 1 衰减到 0后一帧从 0 升到 1线性交叉只在 64 点窗口内做对 QRS 的时域宽度来说这个接缝非常短。注意这种分段方式要求去噪函数完全独立不要在分段间共享滤波器状态。整个 process.m 的脚本结构就是数据读取——参数设定——分段循环去噪——绘图验证四段式。5. 去噪验证与批量处理SNR、PRD指标和并行循环5.1 量化指标怎么算去噪效果不能只看图要算两个数字以残差为噪声的估计noise ecg - ecg_den; SNR 10 * log10(sum(ecg.^2) / sum(noise.^2)); PRD 100 * sqrt(sum(noise.^2) / sum(ecg.^2)); fprintf(SNR %.2f dB, PRD %.2f%%\n, SNR, PRD);SNR 是信噪比的能量比换算成分贝PRD 是百分均方根差两个指标同源只是呈现方式不同。PRD 小于 5% 时波形形态基本满足临床读图要求5%~10% 之间的信号适合做心率变异性分析但要观察残差里是否还有周期性心拍痕迹。若 PRD 偏高且残差出现周期峰降低阈值或换用折中阈值并减小 a若 PRD 偏低但 P 波细节被抹平说明阈值过小噪声和信号被一起压掉了。5.2 批量处理多个数据文件包里同时有 100.dat 和 ECG1.dat 两段数据用循环把两段数据一起跑完是最直接的验证方式files {100.dat, ECG1.dat}; for f 1:length(files) [ecg, fs] load_ecg(files{f}); % 复用 3.1 节的读取函数 [ecg_den, SNR, PRD] ecg_denoise(ecg, haar, 5, soft, 1.0); fprintf(%s: SNR%.2f dB, PRD%.2f%%\n, files{f}, SNR, PRD); end把读取和去噪分别封装成 load_ecg 与 ecg_denoise 两个函数后整段脚本的可读性会明显提升也方便把阈值类型、分解层数作为参数传进来做多组对比。如果记录超过 30 条可以考虑用 parfor 并行循环小波工具箱的函数对并行循环是线程安全的唯一要注意的是每个 worker 都要能访问 Matlab 路径下的小波工具箱部署到集群时记得用 addAttachedFiles 附上自定义函数。有一个非常实用的判断技巧跑完之后把不同 level 下的 SNR 列成一行正常情况下 SNR 随层数先升后降峰值所在层数就是这段信号的最佳分解深度。如果 SNR 单调上升说明信号里低频噪声占主导还需要加一层如果单调下降很可能是阈值取大连 QRS 波都被当成噪声处理掉了。包里那 7 张运行结果图基本就是不同参数组的输出对比重新调参后按这个标准检查去噪前后的 R 波幅值差控制在 10% 以内残差呈现均匀白噪声形态这组参数就可以定下来了。本文还有配套的精品资源点击获取