Haar小波变换心电信号去噪原理与Matlab实现

简介压缩包内含一份基于Haar小波变换的心电信号去噪Matlab源码资源面向本科、硕士阶段进行信号处理与生物医学工程研究的读者可应用于课程设计、毕业设计或论文仿真场景。心电信号采集时常混有工频干扰和肌电噪声该资源演示了如何借助Haar小波变换在时频域分离噪声并保留波形关键特征对理解小波去噪原理很有帮助。包体共14个文件包含可运行的.m源程序、ECGdata.mat与ECG1.dat等实测心电数据、7张jpg及2张png运行效果图同时附说明txt文档和仿真咨询图压缩包整体约3.09MB。目前已有402人学习使用。通过该资源可完整复现心电信号读取、Haar小波分解、阈值去噪、重构及波形对比等流程配合结果图能够直观检验各环节去噪效果源码结构清晰注释明确便于学习者按需修改阈值参数、更换实验数据并进一步开展二次开发。从加载数据到输出波形对比每一步都有对应图像辅助理解可快速掌握小波基选择、分解层数与阈值规则的影响。1. 为什么心电去噪绕不开小波变换一份 10 秒的心电记录里混着 50Hz 工频、肌电和基线漂移时傅里叶滤波经常按下葫芦浮起瓢低通会把 R 波顶点磨圆高通会留下基线台阶滑动平均又会把 ST 段压低。心电本质是非平稳信号P 波、QRS 波群和 T 波在不同时刻有不同的频带固定窗函数很难同时照顾到“突变保持”和“噪声抑制”。Haar 小波变换用一组二进制伸缩平移的基函数把信号拆成多尺度细节噪声与 QRS 突变在不同尺度上分离度很好阈值收缩后重构能在不抹掉波形拐点的情况下把噪声压下去。这篇文章围绕基于 Haar 小波变换的心电信号去噪把阈值去噪原理、Matlab 最小实现、参数调节和指标验证一次讲透。适合做生物医学信号处理、毕业设计预研以及心电采集设备前期算法验证的工程师直接参考。2. Haar小波去噪的原理从分解到阈值重构2.1 小波变换为什么适合心电这种非平稳信号傅里叶变换把整段信号投影到无限长的正弦基上得到的是“全时段平均”的频谱某个时刻的瞬态突变会被稀释到整个频带里短时傅里叶变换加了窗但窗长固定低频需要长窗、高频需要短窗二者只能取一个折中。心电信号的棘手之处在于QRS 波群是毫秒级的陡峭跳变基线漂移是秒级的缓慢波动两者频率可能重叠在同一频段只是持续时间和位置完全不同。小波变换的基函数是“有限长、可伸缩、可平移”的小波簇尺度参数控制频率高低平移参数控制时间位置因此它天然具备时频局部化能力。Haar 小波是最简单的小波基本质上是方波差分计算复杂度低且不会像高阶小波那样在重构时产生过多的人工振荡。对采样率 250Hz 到 1000Hz 的常规心电采集Haar 小波有足够的时域分辨率去定位 R 波顶点去噪后的波形形态也比同阶 Daubechies 小波更容易解释。2.2 Haar小波的尺度函数与分解系数结构Haar 小波的尺度函数在 [0,1) 上恒为 1小波函数在 [0,0.5) 为 1在 [0.5,1) 为 -1。一次分解把信号分成近似系数 cA1 和细节系数 cD1近似系数是相邻两点平均再乘归一化系数细节系数是相邻两点差分。下一层继续对 cA1 分解得到 cA2 和 cD2。噪声能量主要集中在高层细节系数上而 QRS 波群在跨尺度上都有投影尤其是突变沿在 cD1、cD2 里有明显的数值峰值这就是阈值法能区分信号与噪声的基础。在 Matlab 的 Wavelet Toolbox 里wavedec返回的系数数组 C 是按“最后一层细节到第一层细节、再到最后一层近似”的顺序拼接的L 记录每段的长度。要理解去噪脚本首先要能准确切出每一层系数。[C, L] wavedec(ecg, 3, haar); % C [cD3(1...n3), cD2(1...n2), cD1(1...n1), cA3(1...nA3)] % 每层长度 len_d3 L(2); len_d2 L(3); len_d1 L(4); % 用索引切片 cD3 C(1 : L(2)); cD2 C(L(2)1 : L(3)); cD1 C(L(3)1 : L(4)); cA3 C(L(4)1 : L(5));代码里L(1)是原始信号长度L(2)到L(5)分别是前三层细节和第三层近似的长度。切片后可以单独观察某一层系数的数值分布判断噪声集中在哪一层这是调参时的第一手依据。2.3 阈值去噪的标准三步小波去噪不是简单地把高频细节置零那样会把 QRS 波群的陡峭沿一起削掉。标准做法是三步先对含噪信号做多层小波分解再对各层细节系数做阈值收缩最后用处理后的系数重构信号。阈值收缩时低频近似系数一般不动因为它承载的是信号的主体形态真正处理的是 cD1、cD2、cD3 乃至更高层的细节系数。阈值大小的核心是估计噪声标准差。常见做法是用第一层细节系数 cD1 的中位绝对偏差估计sigma median(abs(cD1)) / 0.6745 thr sigma * sqrt(2 * log(N))0.6745 来自正态分布的分位数关系这个公式对应的是“固定阈值” sqtwolog 规则。噪声越强sigma 越大阈值越高信号越长固定阈值也会缓慢增大。Matlab 的thselect函数封装了四种规则实际使用时不必手写公式。2.4 四种阈值规则的选择逻辑规则名称阈值计算方式特点心电场景适用性rigrsureStein 无偏风险估计阈值偏小保留细节多噪声较弱时效果好heursurerigrsure 与 sqtwolog 的启发式组合信噪比低时转向固定阈值心电去噪最常用sqtwolog固定阈值公式稳健但容易过平滑强噪声、基线漂移严重时minimaxi极小极大准则阈值保守保护弱信号P 波、T 波幅度小时优先表格里的“适用性”是相对而言的。心电信号形态固定QRS 幅度远大于噪声时sqtwolog 和 heursure 差别不大但 T 波、ST 段这种低幅度缓变成分对阈值很敏感rigrsure 或 minimaxi 不容易把它们压平。实际项目中我会先用 heursure 跑一版再结合第 5 章的指标判断是否需要换规则。3. 在Matlab里跑通Haar小波心电去噪的最小脚本3.1 准备数据与构造含噪心电信号去噪脚本的输入可以是 MIT-BIH 的 CSV 导出数据也可以是设备采集的 txt 文本。读入后先检查单位很多采集卡输出的是 ADC 码直接做小波分解会发现尺度相差悬殊阈值估计失效。常见做法是先转换为 mV把原始数值减去基线偏置再除以增益系数。下面的代码示例用readmatrix读入 CSV然后人工叠加 50Hz 工频和随机肌电干扰方便在已知干净信号的情况下评估效果。data readmatrix(ecg_sample.csv); % 第一列时间第二列心电 fs 360; % MIT-BIH 采样率 ecg data(:, 2) * 0.001; % 假设原始单位为 uV转成 mV t (0:length(ecg)-1) / fs; N length(ecg); % 模拟噪声50Hz工频 高斯白噪声 noise_50 0.08 * sin(2 * pi * 50 * t); noise_emg 0.15 * randn(N, 1); ecg_noisy ecg noise_50 noise_emg;采样率 fs 必须与滤波器参数匹配工频是 50Hz 还是 60Hz 取决于所在地区。0.001的单位换算系数需要根据实际硬件增益调整不换算会造成阈值整体偏小或偏大。叠加噪声的幅度要与真实采集场景接近否则去噪效果看起来很好换到真数据立刻失效。3.2 用wavedec做Haar分解分解层数先取 5 层原因是 360Hz 采样率下第 5 层细节大约对应 5.6Hz 到 11.2Hz 附近可以覆盖心电主频段的一部分同时把大部分随机噪声和工频干扰分离到前两层细节中。level 5; wname haar; [C, L] wavedec(ecg_noisy, level, wname);wavedec返回的 C 是行向量L 是长度向量。Haar 小波在这里不需要指定滤波器长度因为它是内置小波基。如果信号长度不是 2 的整数次幂wavedec会自动做边界延拓后续第 4 章会专门处理边界问题。3.3 用wthresh做软阈值处理阈值估计分两步先用第一层细节系数算噪声标准差再用wthresh对系数做收缩。软阈值函数s会把所有系数绝对值减掉阈值后归零硬阈值h则保留超过阈值的部分。心电信号我一般默认软阈值因为硬阈值在重构后容易在 R 波附近产生振铃。cD1 C(L(3)1 : L(4)); % 第1层细节 sigma median(abs(cD1)) / 0.6745; thr sigma * sqrt(2 * log(N)); C_thr C; idx L(2) : length(C); % 从第5层细节到第1层细节 C_thr(idx) wthresh(C_thr(idx), s, thr);idx从L(2)开始是为了把 5 层细节全部纳入阈值处理L(2)指向第 5 层细节的起始位置。低频近似系数cA5保留原样否则重构出来的波形会失去基线形态。阈值thr是全局标量若想逐层用不同阈值需要换成循环结构这个变体放在第 4 章讨论。3.4 用waverec重构并输出结果重构时传入收缩后的系数 C_thr 和原有的 L不需要额外指定小波基因为 L 里已经隐含了分解层数信息。重构完成后用plot对比三组信号原始干净信号、含噪信号、去噪信号。ecg_denoised waverec(C_thr, L, wname); figure; subplot(3,1,1); plot(t, ecg); title(原始心电); subplot(3,1,2); plot(t, ecg_noisy); title(含噪心电); subplot(3,1,3); plot(t, ecg_denoised); title(Haar小波去噪结果);waverec只负责重构不会检查阈值是否合理。如果去噪后波形整体缩小了一个量级通常是单位换算错误或阈值过大导致细节系数被过度收缩如果波形出现明显台阶说明某层系数切片索引写错了细节系数和近似系数混在一起被置零。逐层打印L和length(C)是排查这类问题最快的方式。4. 层数、阈值和边界怎么调参数实验与避坑4.1 分解层数选择经验公式与能量占比分解层数是最容易拍脑袋的参数。层数太少噪声和基线漂移分离不干净层数太多最后一层细节已经包含了大部分 QRS 能量收缩后会削掉信号本身。常见经验公式是floor(log2(N)) - 2它给出最大可分解层数但实际心电并不需要取满。对 360Hz、10 秒的信号N3600最大层数为 9我一般取 5 到 8 层。更稳妥的方法是用各层细节能量占比判断如果第 k 层细节能量占总能量的比例已经很小说明继续分解下去信息增量有限。下面的代码遍历 1 到 8 层计算每层细节能量占比。levels 1:8; ratio zeros(size(levels)); for k levels [Ck, Lk] wavedec(ecg_noisy, k, haar); det_energy 0; total_energy sum(Ck.^2); for j 1:k idx Lk(j)1 : Lk(j1); det_energy det_energy sum(Ck(idx).^2); end ratio(k) det_energy / total_energy; end运行后如果发现ratio(6)到ratio(8)已经平稳取那个拐点对应的层数即可。这里用平方和表示能量是因为小波系数是正交变换Parseval 定理保证系数能量等于信号能量。需要注意这个评估依赖含噪信号的噪声水平更换数据集后拐点会移动不要照搬某一篇文章里的固定层数。4.2 阈值规则与层数的组合实验不同的阈值规则和分解层数之间存在交互层数浅时rigrsure 的保守阈值可能残留较多噪声层数深时sqtwolog 的全局阈值会把 QRS 高频分量压掉。比较高效的方案是做一次小网格扫描用第 5 章的指标函数自动选优。rule_list {rigrsure, heursure, sqtwolog, minimaxi}; for level 4:7 for r 1:length(rule_list) ecg_d wden(ecg_noisy, rule_list{r}, s, mln, level, haar); snr_val calc_snr(ecg, ecg_d); % 自定义指标函数见第5章 fprintf(level%d, rule%s, SNR%.2f dB\n, ... level, rule_list{r}, snr_val); end endwden是wavedec、阈值估计和waverec的一站式封装六个参数分别是信号、阈值规则、软硬阈值、噪声尺度、层数、小波基。mln表示用小波分解后的各层噪声水平来做阈值调整比全局阈值更精细也是与“逐层阈值”最接近的封装方式。如果只看单一参数组合很难判断问题出在层数还是阈值规则上网格扫描打印的 SNR 表格能直接给出结论。4.3 边界效应与延拓方式小波分解在信号两端会引入假系数默认延拓模式是sym对称延拓对心电这种首尾不连续的波形重构后的前几十个采样点经常出现上升或下降的假象。可以用dwtmode查看当前延拓模式常用备选还有zpd零填充和ppd周期延拓。实际处理中我一般会先在原始信号两端各延拓一小段数据去噪完成后再切掉。例如以信号均值填充 0.5 秒左右的前后段让边界处的系数分布更接近信号内部避免边界失真污染 QRS 波。4.4 硬阈值振铃、软阈值过平滑与分层阈值硬阈值重构容易在 R 波两侧产生伪吉布斯振荡表现为高频毛刺软阈值则会把 T 波起点和终点“磨钝”。折中方案是分层设阈值前两层用较小的硬阈值保护 QRS 突变沿后几层用较大的软阈值过滤噪声和基线漂移。C_thr C; for j 1:level idx L(j)1 : L(j1); layer_sigma median(abs(C(idx))) / 0.6745; if j 2 layer_thr layer_sigma * sqrt(2 * log(N)) * 0.7; C_thr(idx) wthresh(C(idx), h, layer_thr); else layer_thr layer_sigma * sqrt(2 * log(N)); C_thr(idx) wthresh(C(idx), s, layer_thr); end end第 1、2 层细节系数承载了 QRS 波群的尖锐沿硬阈值配合 0.7 的折扣系数能保留更多突变信息第 3 层及以后以消噪为主软阈值更平稳。这里的 0.7 和“前两层硬阈值”是经验起点最优分层策略可以通过网格扫描确定但至少比单一全局阈值稳定得多。5. 去噪效果的量化验证与自动选参技巧5.1 SNR、RMSE、PRD三个指标的计算波形肉眼看着干净不算数需要用指标量化。信噪比 SNR、均方根误差 RMSE、百分误差 PRD 是心电去噪论文和工程验证中最常用的三个指标它们的计算都要求有干净参考信号。这个前提只存在于仿真和公开数据集中真实采集环境没有无噪真值只能退而求其次用残差平稳性判断。function [snr, rmse, prd] ecg_eval(orig, denoised) noise orig - denoised; snr 10 * log10(sum(orig.^2) / sum(noise.^2)); rmse sqrt(mean(noise.^2)); prd 100 * sqrt(sum(noise.^2) / sum(orig.^2)); endSNR 提高 3dB 以上通常能听出明显效果但对心电来说SNR 高不等于诊断信息完整ST 段抬高这类临床特征可能在 SNR 提升的同时被扭曲。RMSE 对绝对幅度敏感两个不同增益系统的 RMSE 不可直接比较PRD 是归一化指标适合不同数据段之间做横向对比。5.2 用残差自相关判断信号是否有泄漏没有真值的情况下重点检查残差含噪信号减去除噪信号是否接近白噪声。对残差做自相关如果只有零延迟处有峰值两侧快速衰减说明去噪把信号和噪声分离得比较干净如果残差在 R 波位置附近有明显周期性峰值说明一部分心电能量被当成了噪声。也可以用 FFT 看残差的频谱50Hz 处若有明显尖峰说明工频没有滤除干净低频段若有大面积突出的能量说明基线漂移只压掉了一部分。这两个检查和 SNR 结合能判断当前参数是“过度去噪”还是“去噪不足”。5.3 一个实用技巧用循环自动选最优分解层数自动选参的常规思路是把层数、阈值规则、软硬阈值三个变量做成嵌套循环用无噪参考信号计算 SNR取最优组合。这里有一个细节层数搜索范围不要从 1 开始至少从 4 开始因为低层数下噪声和信号在频带上重叠严重SNR 峰值往往在 5 到 7 层之间。best_level 1; best_snr -inf; for level 4:8 ecg_d wden(ecg_noisy, heursure, s, mln, level, haar); snr_val 10 * log10(sum(ecg.^2) / sum((ecg - ecg_d).^2)); if snr_val best_snr best_snr snr_val; best_level level; end end这段代码把“调参”变成了一个可重复执行的搜索过程。得到最优层数后再固定层数去扫描阈值规则会比同时扫描所有组合更快。需要注意的是自动选出的最优参数在另一段心电数据上未必同样最优跨患者数据验证时应在 5 到 10 条记录上分别运行同一流程取表现最稳定的参数组合投入使用。本文还有配套的精品资源点击获取