简介面向信号与图像去噪学习者的MATLAB小波去噪代码包适合已具备MATLAB基础、希望理解多尺度分解、阈值收缩与重构流程的高校学生和工程技术人员。内容围绕小波分解、噪声估计、软硬阈值处理、系数调整及逆变换重构展开并配有PSNR、MSE等质量评价脚本便于在lena图像上复现实验、比较不同阈值策略的效果。压缩包共9个文件约708KB其中6个m脚本承担算法实现与指标计算1个pdf提供自适应小波阈值去噪与压缩的参考材料png和bmp各1个用作测试图像体量轻便、结构紧凑。已有3121人学习说明该案例在课程实验与入门研究中具有一定参考价值。读者可据此掌握小波基选择、分解层数和阈值设定对去噪结果的影响并在此基础上尝试贝叶斯收缩、自适应阈值等改进方案积累可复用的实验代码与调参思路。1. 从一段毛刺密布的振动信号说起小波去噪到底解决什么问题工位上拿到一段 1 kHz 采样的振动信号波形上全是毛刺FFT 一看主频之外铺着一层底噪直接做包络分析根本看不出特征。这种场景就是这份小波去噪 MATLAB 代码的典型用武之地。它把一维信号和二维图像的去噪流程拆成了能直接运行的脚本主干链路是wavedec多尺度分解、wthresh阈值收缩、waverec重构外层再包一层wden、wdencmp做自动阈值。开发语言是 MATLAB源码里参数集中在脚本开头小波基、层数、阈值规则都是变量改完按运行就能看到波形变化不需要从头推多分辨率分析。做信号处理课设的学生、要快速验证去噪效果的算法工程师、想拿一套基线做对比实验的研究者都能从这份代码里省下搭框架的时间。下面按「跑起来、拆链路、改参数」的顺序讲。2. 工程结构与环境把这份 MATLAB 源码跑起来的第一步拿到一份 MATLAB 源码包最容易翻车的地方不是算法本身而是路径没加对、工具箱没装全脚本一运行就报Undefined function。这一章先把目录职责、依赖环境和冒烟测试理顺后面拆代码才有落脚点。2.1 目录结构与文件清单这类小波去噪源码常见的组织方式是按「入口脚本 函数库 数据 输出」四层来分我一般会整理成下面这个样子。如果你的压缩包解压出来文件名略有出入按职责对应即可重点看denoise/下那几个函数是不是被入口脚本调用。wavelet_denoise/ ├── main_1d.m # 一维信号去噪演示入口 ├── main_2d.m # 二维图像去噪演示入口 ├── denoise/ │ ├── denoise_1d.m # 手动三段式封装分解-阈值-重构 │ ├── denoise_2d.m # 二维图像去噪封装 │ └── calc_metrics.m # SNR / PSNR / SSIM 计算 ├── data/ │ └── sample_signal.mat # 示例信号含干净信号与含噪信号 └── results/ # 运行后输出的图与指标文本上面这份清单里main_1d.m和main_2d.m是给人看的演示入口真正干活的是denoise/下的函数。sample_signal.mat里一般会存两组变量一组是理想干净信号用来算 SNR 做定量评估另一组是叠加了噪声的观测信号作为算法输入。results/目录通常由脚本自己创建saveas或print把图存进去。需要提醒的是MATLAB 不会自动递归搜索子文件夹入口脚本开头一般要加一句addpath(fullfile(pwd, denoise))否则main_1d.m调用denoise_1d时会提示找不到函数。这一点在换电脑、换路径之后尤其容易忘属于新手最先撞上的墙。2.2 运行环境与工具箱依赖这份代码的核心函数wavedec、waverec、wthresh、wden、wdencmp、wdcbm2全部来自 Wavelet Toolbox不是 MATLAB 基础安装自带的。只装了 MATLAB 主程序、没勾选工具箱的话脚本跑到第一行分解就会停。图像去噪部分还会用到imread、rgb2gray、imshow这些属于 Image Processing Toolbox 的常用函数部分基础函数在核心 MATLAB 里也有但ssim明确依赖 Image Processing Toolbox。版本方面没有硬性门槛这些函数的签名在 R2016b 之后的版本里基本稳定老版本 MATLAB 也能跑只是ssim的调用形式在不同版本间有过调整碰到参数报错就查一下本地help ssim。运行前先用几行代码确认工具箱在位比跑一半报错再回头查要省事% 检查 Wavelet Toolbox 是否可用 v ver(wavelet); if isempty(v) error(未检测到 Wavelet Toolboxwavedec / wden 等函数不可用请先安装该工具箱); end fprintf(Wavelet Toolbox 版本: %s\n, v.Version); % 图像去噪脚本额外检查 Image Processing Toolbox v2 ver(images); if isempty(v2) warning(未检测到 Image Processing Toolboxssim 等函数可能不可用); end这段代码的逻辑很直白ver(wavelet)返回一个结构体数组工具箱没装时返回空直接error中断把问题暴露在启动阶段而不是运算中途。v.Version打印出来方便排查「同一份代码在同事机器上能跑、在我这里不行」这类环境差异。第二段对 Image Processing Toolbox 只做warning因为一维去噪并不需要它没必要因为图像部分缺工具箱就挡住整条链路。2.3 用一段合成信号做冒烟测试环境确认之后先别急着上真实采集数据。合成一段已知干净波形的含噪信号跑一遍能同时验证工具箱、函数路径和基本参数是否正常这是成本最低的冒烟测试。% 冒烟测试合成含噪信号跑自动去噪比较 SNR fs 1000; % 采样率 1 kHz t (0:999) / fs; % 1 秒时长 clean sin(2*pi*50*t) 0.6*sin(2*pi*120*t); % 两个正弦叠加作为理想信号 rng(7); % 固定随机种子保证结果可复现 noisy clean 0.25*randn(size(t)); % 叠加高斯白噪声 % wden 自动去噪sqtwolog 阈值规则、软阈值、单层噪声缩放、4 层分解、db4 小波 den wden(noisy, sqtwolog, s, sln, 4, db4); snr_in 10*log10(sum(clean.^2) / sum((noisy - clean).^2)); snr_out 10*log10(sum(clean.^2) / sum((den - clean).^2)); fprintf(输入 SNR %.2f dB输出 SNR %.2f dB\n, snr_in, snr_out);这段脚本的关键在于wden的六个参数顺序是「信号、阈值规则、软硬阈值、缩放方式、分解层数、小波基」。sqtwolog对应通用阈值sqrt(2*log(N))s表示软阈值sln表示用第一层细节系数估计噪声标准差做单层缩放4是分解层数db4是小波基。SNR 的计算公式是10*log10(信号功率/噪声功率)输入 SNR 一般在几 dB 到十几 dB 之间去噪后提升 3 到 8 dB 属于正常范围如果输出 SNR 反而下降多半是阈值取得过大把信号成分也削掉了这时候先降层数、再换阈值规则排查。冒烟测试通过说明环境没问题可以进入核心链路的拆解。真实数据往往没有干净参考SNR 算不了但这一步建立的手感——什么参数下波形大概长什么样——在后面调参时很有用。3. 一维链路拆解wavedec、阈值收缩与waverec重构的三段式wden一行调用看着简单内部其实把分解、阈值、重构串成了一条流水线。想改得动参数、出了问题能定位得把这条链路拆开看。这一章按「分解、阈值、重构、自动封装」四段讲每段给出可单独运行的代码。3.1 分解层数与小波基的选择逻辑wavedec做的是多尺度分解把信号逐层拆成低频近似系数和高频细节系数每层细节系数对应一个频段。返回两个东西系数向量c和长度向量l。c的排布顺序是[cA_n, cD_n, cD_{n-1}, ..., cD_1]l记录每一段占多少元素sum(l)等于length(c)。这个结构不理解后面改系数时很容易张冠李戴。% 多尺度分解看清系数向量的排布 wname db4; % 小波基消失矩 4 maxlev wmaxlev(numel(noisy), wname); % 该小波基允许的最大分解层数 level min(4, maxlev); % 实际取 4 层但不超过上限 [c, l] wavedec(noisy, level, wname); fprintf(最大可分解层数 %d实际层数 %d\n, maxlev, level); fprintf(系数向量长度 %d长度向量和 %d\n, numel(c), sum(l)); % l(1) 是近似系数 cA_level 的长度l(2) 起依次是 cD_level ... cD_1wmaxlev根据信号长度和小波滤波器长度算出一个理论上限一般是floor(log2(N/(L-1)))量级其中L是小波滤波器长度。层数不是越多越好每分解一层细节系数长度减半到高层只剩几十个点用来估计噪声标准差的样本量不够阈值会抖同时边界效应逐层向内扩散重构后两端的失真会更明显。我一般取 3 到 5 层信号短于 256 点时取 3 层采样率高、噪声集中在高频段时可以到 5 层。小波基的选择直接影响去噪效果和边界失真。消失矩越高对平滑信号的逼近越好但滤波器越长边界处理涉及的点越多。下面这张表是几种常用基的取舍实际项目里我通常从sym4和db4起步再往两端试。小波基正交性支撑长度适用场景db4正交8通用起点平滑信号计算量小db8正交16信号较平滑、对逼近精度要求高sym4正交8近似对称边界失真比 db 系列小coif3正交18兼顾平滑性与对称性适合生物信号bior3.5双正交分解/重构不等长需要线性相位时用图像处理更常见选基的时候别只看论文里写哪个把两三种基跑同一段数据、比 SNR 和波形边缘十分钟就能有结论。db系列在信号端点处容易出现不对称的振铃sym系列对称性更好端点失真会轻一些。3.2 阈值函数软阈值、硬阈值与折中方案分解完噪声主要落在各层细节系数上阈值收缩就是把幅值小的系数往零压。wthresh提供两种基本方式软阈值和硬阈值。软阈值的处理是sign(c) * max(|c| - thr, 0)系数整体向零收缩函数连续但大系数也被砍掉一个固定量重构后幅值偏小硬阈值的处理是c .* (|c| thr)小于阈值的直接置零大于的原样保留幅值不偏但在阈值点不连续容易产生振铃。% 软阈值 vs 硬阈值同一组系数两种收缩方式 thr 0.4; % 阈值实际应由噪声估计给出 c_soft wthresh(c, s, thr); % 软阈值系数向零收缩 c_hard wthresh(c, h, thr); % 硬阈值小幅系数直接置零 fprintf(软阈值非零系数占比 %.1f%%硬阈值 %.1f%%\n, ... 100*sum(c_soft ~ 0)/numel(c_soft), 100*sum(c_hard ~ 0)/numel(c_hard));阈值thr怎么定才是去噪效果的分水岭。MATLAB 内置四种规则wden的第二个参数就是选它% 四种阈值规则对比同一段信号只看 SNR 变化 rules {sqtwolog, rigrsure, heursure, minimaxi}; for k 1:numel(rules) xd wden(noisy, rules{k}, s, sln, 4, db4); snr 10*log10(sum(clean.^2) / sum((xd - clean).^2)); fprintf(%-10s 输出 SNR %.2f dB\n, rules{k}, snr); endsqtwolog是通用阈值sqrt(2*log(N))实现简单噪声强时偏保守容易过平滑rigrsure基于 Stein 无偏风险估计逐系数算风险适合噪声水平不确定的场景但系数少的时候估计不稳heursure是前两者的折中信号信噪比高时偏向rigrsure低时偏向sqtwolog我一般拿它当默认minimaxi用极小极大准则适合信号本身较平滑、只想去掉小幅噪声的情况。真实项目里没有干净参考不能每换一次规则就算 SNR那就看波形过平滑的表现是尖峰被削圆、上升沿变缓欠去噪的表现是波形上还挂着一层高频毛刺。3.3 重构与边界效应处理重构用waverec输入是处理后的系数向量、长度向量和小波基名。这里有个容易忽略的点延拓模式是全局设置靠dwtmode控制wavedec和waverec都读这个全局状态而不是通过函数参数传。分解前设成什么模式重构时就得保持一致。% 边界延拓模式设置必须在分解之前调用 dwtmode(sym, nodisp); % 对称延拓默认值端点失真较小 % dwtmode(per, nodisp); % 周期延拓要求信号长度为 2 的幂次更稳 % dwtmode(zpd, nodisp); % 零填充端点容易突变一般不推荐 [c, l] wavedec(noisy, 4, db4); thr 0.5 * sqrt(2*log(numel(noisy))); % 通用阈值再乘个经验系数 c_thr wthresh(c, s, thr); x_den waverec(c_thr, l, db4); % 长度对齐重构长度偶尔会因延拓多出几个点按原长截断 if numel(x_den) ~ numel(noisy) x_den wkeep(x_den, numel(noisy)); fprintf(重构长度已从 %d 截断到 %d\n, numel(x_den)abs(numel(x_den)-numel(noisy)), numel(noisy)); enddwtmode(sym)是对称延拓信号两端按镜像方式补点重构后端点失真最小是绝大多数场景的默认选择。per是周期延拓频域分析里常用但它要求信号长度接近 2 的整数次幂否则内部补零会让重构长度变长这是后面避坑章节要专门讲的一条。wkeep按指定长度从中间截取是长度对不齐时的补救手段。阈值前面的系数0.5是我自己加的。通用阈值在噪声估计偏大的时候会过平滑乘一个 0.5 到 0.8 的经验系数能在保留细节和压制噪声之间找个平衡。这个系数没有理论最优值得看具体信号我在振动信号上常用 0.6在语音上常用 0.4因为语音的爆破音容易被削掉。3.4 wden 自动封装一次调用背后的参数映射把上面三段串起来就是wden干的事。它的调用形式是xd wden(x, tptr, sorh, scal, n, wname)六个参数分别对应阈值规则、软硬阈值、缩放方式、分解层数和小波基。% 自动封装等价于分解 阈值 重构但把噪声估计也代劳了 xd wden(noisy, heursure, s, mln, 4, sym4); % tptr heursure 阈值规则 % sorh s 软阈值 % scal mln 按每层细节系数分别估计噪声适合非平稳噪声 % n 4 分解层数 % wname sym4 小波基scal这个参数值得单独说。取one表示不缩放阈值直接按规则算取sln表示用第一层细节系数的噪声估计统一缩放所有层取mln表示每一层用自己的细节系数估计噪声分别缩放。信号里噪声水平沿时间轴变化时mln更稳噪声平稳时sln就够计算量还小一点。wden方便但它是黑盒参数改错了不容易定位是哪一段出的问题所以调参阶段我一般先用wden快速扫一遍锁定大概范围后再回到wavedec手动链路里逐层看系数。4. 二维图像去噪系数布局、wdcbm2自适应阈值与PSNR评估一维信号是一条线二维图像是一张面分解后的系数从「近似 各级细节」变成「近似 每级的水平、垂直、对角三个方向细节」。函数换了思路没换但系数布局和阈值策略有不少差别这一章单独拆。4.1 二维小波分解的系数布局wavedec2返回的系数向量c按[cA_n, cH_n, cV_n, cD_n, cH_{n-1}, cV_{n-1}, cD_{n-1}, ..., cH_1, cV_1, cD_1]排列cH是水平细节抓水平边缘cV是垂直细节抓垂直边缘cD是对角细节抓斜向纹理。尺寸矩阵s是(n2)×2第一行是原图尺寸之后每行是各级近似系数的尺寸。这个布局不清楚后面想只处理某一方向的细节就无从下手。% 二维分解看系数排布与尺寸矩阵 img im2double(imread(cameraman.tif)); % 示例图换成自己的灰度图 wname sym4; level 3; [c, s] wavedec2(img, level, wname); fprintf(系数向量长度 %d\n, numel(c)); fprintf(尺寸矩阵行数 %d应为 level2 %d\n, size(s,1), level2); fprintf(原图尺寸 %d x %d\n, s(1,1), s(1,2));im2double把图像转成[0,1]的 double 类型这一步不能省。原始imread读进来是uint8直接往上加高斯噪声会溢出截断去噪前图像就已经废了。s(1,:)存的是原图尺寸重构和评估都以它为准后面算 PSNR 时归一化分母也用这个。4.2 wdencmp 与 wdcbm2 的自适应阈值图像去噪有两条常用路径。一条是ddencmp给默认参数、wdencmp执行全局阈值胜在快另一条是wdcbm2按层给自适应阈值胜在细节保留。两种都值得跑一遍对比。% 路径一ddencmp wdencmp全局阈值 [thr, sorh, keepapp] ddencmp(den, wv, noisy_img); img_den1 wdencmp(gbl, noisy_img, wname, level, thr, sorh, keepapp); % 路径二wdcbm2 按层自适应阈值需要先分解拿到 c 和 s [c, s] wavedec2(noisy_img, level, wname); alpha 1.5; % 稀疏参数越小保留系数越少、去噪越狠 m 3 * prod(s(1,:)); % 每层保留系数个数的上界 thr_vec wdcbm2(c, s, alpha, m); % 返回逐层阈值向量 img_den2 wdencmp(lvd, c, s, wname, level, thr_vec, s);ddencmp(den,wv,X)三个返回值分别是阈值、软硬阈值标记、是否保留近似系数den表示去噪、wv表示小波。wdencmp(gbl,...)用单一阈值处理所有细节系数keepapp取 1 时近似系数原样保留不参与收缩。路径二里wdcbm2基于 Birgé-Massart 策略alpha控制每层保留系数的稀疏程度常取 1.5 到 3m给的是系数个数上界不是噪声标准差这一点最容易写错3*prod(s(1,:))是文档示例里的常用取值。wdencmp(lvd,...)的lvd表示层相关阈值此时thr必须是和层数对应的向量不能传标量。全局阈值实现简单但图像不同子带的噪声能量分布不一样一刀切容易在平滑区域留下块状伪影层相关阈值按子带分别设细节保留更好代价是要调alpha和m两个参数。我一般先用路径一跑个基线再换路径二看提升PSNR 涨 1 dB 以上就值得多花这点调参时间。4.3 去噪效果量化与可视化图像去噪有干净参考图时PSNR 和 SSIM 是两个直接指标。PSNR 衡量像素级误差SSIM 衡量结构相似度前者对亮度敏感后者更贴近人眼感受。% 定量评估PSNR 与 SSIM mse1 mean((img(:) - img_den1(:)).^2); psnr1 10*log10(1 / mse1); % 图像已归一化到 [0,1]分母取 1 ssim1 ssim(img_den1, img); % 需要 Image Processing Toolbox mse2 mean((img(:) - img_den2(:)).^2); psnr2 10*log10(1 / mse2); ssim2 ssim(img_den2, img); fprintf(全局阈值 : PSNR %.2f dB, SSIM %.4f\n, psnr1, ssim1); fprintf(层相关阈值: PSNR %.2f dB, SSIM %.4f\n, psnr2, ssim2); % 可视化原图、含噪图、两种去噪结果并排 figure(Position, [100 100 1000 300]); subplot(1,4,1), imshow(img), title(原图); subplot(1,4,2), imshow(noisy_img),title(含噪); subplot(1,4,3), imshow(img_den1), title(sprintf(全局 %.2fdB, psnr1)); subplot(1,4,4), imshow(img_den2), title(sprintf(层相关 %.2fdB, psnr2));PSNR 的分母取 1是因为图像已经im2double归一化到[0,1]峰值信号功率就是 1。如果忘了归一化、直接对uint8算PSNR 会虚高二十多 dB对比就失去意义。SSIM 返回[0,1]之间的值越接近 1 结构越相似去噪后一般在 0.7 到 0.95 之间具体看噪声强度和图像纹理。可视化不只是为了出图。四种图并排看能发现指标反映不出来的问题块状伪影、边缘模糊、局部过平滑这些都是 PSNR 涨了但观感变差的情况。指标用来筛参数肉眼用来做最终判断两个都看才稳。5. 避坑与排查小波去噪翻车的六个现场小波去噪的代码量不大坑却不少而且大多不是语法错误是参数和数据层面的隐性问题跑完不报错但结果不对。下面六条是我在不同项目里反复遇到的按「现象、原因、解决」写清楚碰到了直接对号入座。5.1 边界与长度对不齐现象一重构信号比原始信号长出十几个点画图时时间轴对不上。原因多半是dwtmode设成了per而信号长度不是 2 的整数次幂wavedec内部按周期延拓补了零waverec把补出来的部分也重构回去了。解决办法是把模式改回dwtmode(sym,nodisp)或者重构后加一句x_den wkeep(x_den, numel(x))按原长截断。前者从源头避免后者是补救我一般两个都做防止别人改代码时又把模式改回去。现象二图像去噪后尺寸和原图差几个像素imshow显示时边缘被裁。原因是手动改了尺寸矩阵s或者对c做了截断补零导致waverec2按错误的尺寸重构。s是分解时生成的元数据除了第一行原图尺寸要读其余不要手改。重构后尺寸不对用s(1,:)作为目标尺寸做imresize或裁剪别凭感觉补。5.2 阈值与效果类问题现象三去噪后信号过于平滑尖峰和突变全被削圆。原因是软阈值加通用阈值组合阈值偏大把幅值较小的真实成分也收缩掉了。先换rigrsure或heursure规则再把阈值乘一个 0.5 到 0.8 的系数层数从 4 降到 3三招里总有一招见效。判断标准是看上升沿过平滑时上升沿明显变缓和原始信号的陡峭程度差太多。现象四噪声压下去一些但高频细节也糊了振铃明显。硬阈值在阈值点不连续重构后容易出现振铃小波基消失矩太低对细节的逼近能力不够。把h换成s或者把db2换成sym8、coif5振铃会明显减轻。代价是计算量上升滤波器变长边界处理涉及的点更多短信号上要权衡。现象五图像去噪后出现块状伪影平滑区域一块一块的。全局阈值对不同子带用同一个阈值低频子带噪声少、被过度收缩高频子带噪声多、收缩不够重构后就出现块效应。换成wdencmp(lvd,...)配合wdcbm2做层相关阈值或者至少把keepapp设为 1 保住近似系数块状伪影会缓解很多。5.3 环境与数据类问题现象六运行报Undefined function wavedec或者ssim找不到。前者是 Wavelet Toolbox 没装后者是 Image Processing Toolbox 没装。用ver命令逐个确认缺哪个装哪个实在装不了工具箱一维去噪可以改用 Python 的 PyWavelets 做对照验证但函数名和参数顺序都不一样不能直接照搬。还有一种少见情况是工具箱装了但 License 没生效license(test,Wavelet_Toolbox)返回 0这种要找管理员处理不是代码问题。**现象七图像加噪后一片本文还有配套的精品资源点击获取
