MATLAB图像滤波与去噪算法代码包全解析：从空间域到频域的工程实践

简介：面向图像处理科研与工程应用的MATLAB滤波去噪资料包，整合了双线性滤波、Kirsch滤波、逆滤波、双边滤波、同态滤波、小波滤波、约束最小平方滤波、非线性复扩散滤波、Lee滤波、Gabor滤波、Wiener滤波、Kuwahara滤波、Beltrami流滤波、Lucy-Richardson滤波、NonLocalMeans滤波等主流算法，适合正在学习图像去噪、复原或开展算法对比的本科生、研究生及相关开发者使用。压缩包内共217个文件，以215个m脚本为主，另含1个txt说明与1个fig交互界面文件，txt提供使用指引，fig可直接打开交互演示界面，能从算法调用、参数调试到GUI演示形成完整流程；整体仅77KB，轻量易部署。平台显示已有558人浏览学习。借助这些实现，读者能快速获得各类滤波器的可运行代码、统一测试入口及可视化框架，便于逐一复现算法效果、对比去噪性能，并迁移到自己的图像处理课题中。

1. 滤波这行当，MATLAB 代码包里到底藏了多少种玩法

把这份代码包里十几个滤波脚本全部改造成同一个测试接口之后，我的第一个结论可能和大多数人想的不一样：在低噪声密度、单一高斯噪声的场景下，花力气实现的 Richardson-Lucy、Beltrami 流这类高级滤波，和最简单的均值滤波拉不开本质差距；真正让它们分出胜负的，是混合噪声、纹理区域保真和边缘过冲这三件事。代码包里的Draw_Function_GUI.m、Run_Draw.m和一堆编号脚本（3-2.m、15-1.m、17-4.m）覆盖了从空间域到频域、从去卷积到变分法去噪的完整链路，适合做课程设计、论文复现，也适合在 MATLAB 图像处理大作业里直接抽某个模块改参数跑实验。下面按滤波器族分组拆解，重点说清每个方法的适用边界和 MATLAB 实现时的关键参数，最后给出一套横向评估脚本。

2. 空间域滤波怎么选：均值、中值到超限邻域与 Kuwahara 的取舍逻辑

2.1 邻域统计类滤波的统一实现框架

空间域滤波的本质是对每个像素的邻域做加权统计，差别只在权重和统计量上。代码包里3-4.m、10-2.m这类脚本属于同一族，可以用一个统一框架来跑实验，避免每个文件复制粘贴改参数：

% filter_survey.m % 统一测试入口：同一张图、同一噪声，跑多个空间域滤波器 img = im2double(imread('lena.png')); if size(img,3) == 3, img = rgb2gray(img); end rng(2024); noisy = imnoise(img, 'gaussian', 0, 0.01); % 高斯噪声，方差 0.01 % 1) 均值滤波：3x3 邻域 h = fspecial('average', [3 3]); out_avg = imfilter(noisy, h, 'replicate'); % 2) 中值滤波：3x3 邻域 out_med = medfilt2(noisy, [3 3], 'symmetric'); % 3) 超限邻域滤波：中心与邻域均值差超过阈值才替换 th = 0.05; mu = imfilter(noisy, h, 'replicate'); diffmap = abs(noisy - mu); out_lim = noisy; out_lim(diffmap > th) = mu(diffmap > th);

参数说明：imfilter的边界选项replicate对图像边缘的处理方式是复制边缘像素，比默认补零少引入黑边；medfilt2的symmetric则是对边界做镜像扩展，在边缘纹理场景下更稳。超限邻域滤波的阈值th是关键，它决定「多大差异算噪点」：阈值太小会把真实边缘也抹掉，太大则退化成不滤波。代码包里对这组实验的处理方式是固定噪声方差，再对th做 0.01 到 0.1 的扫描，生成一组对比图，直接在Run_Draw.m里以子图网格展示。

2.2 Kirsch 方向模板在滤波里的角色

Kirsch 滤波经常被误认为是纯边缘检测算子，实际上代码包把它归入滤波模块是有道理的：它对八个方向分别做卷积，取响应最大的方向作为该像素的「主方向」，然后用方向自适应权重做平滑。这属于边缘保持滤波的一个朴素实现，比后来基于偏微分方程的方法计算量小一个量级。

% kirsch_dir.m % 8 个方向模板，只列前两个示例，其余按 45 度旋转生成 k1 = [5 5 5; -3 0 -3; -3 -3 -3]; % 北方向 k2 = [5 5 -3; 5 0 -3; -3 -3 -3]; % 东北方向 kernels = zeros(3, 3, 8); for d = 1:8 kernels(:,:,d) = imrotate(k1, (d-1)*45, 'crop', 'bilinear'); kernels(:,:,d) = kernels(:,:,d) / sum(abs(kernels(:,:,d)(:))); % 归一化 end resp = zeros(size(noisy,1), size(noisy,2), 8); for d = 1:8 resp(:,:,d) = conv2(noisy, kernels(:,:,d), 'same'); end [~, dirIdx] = max(abs(resp), [], 3); % 每个像素的主方向编号

这里imrotate用来生成旋转后的方向模板，但插值会引入轻微误差；我通常直接用预定义数组存 8 个方向。方向编号dirIdx的意义在于，后续平滑时只沿主方向及其垂直方向做加权平均，而不是全向平滑，这样能保留线条类纹理。代码包里ssbm.m应该是「双边滤波 + 方向模板」的变体，区别于 2.4 节的标准双边滤波。

2.3 Kuwahara 滤波与局部统计量

Kuwahara 滤波的核心思路：把像素邻域切成四个重叠的子块，分别计算均值和方差，输出方差最小子块的均值。原理是方差小的区域更均匀，噪声被抑制的概率更大，同时边缘两侧的方差差异会让输出偏向某一侧，从而保住边缘。

% kuwahara_impl.m % 窗口尺寸 5x5，分成 4 个 3x3 子块 imgP = padarray(noisy, [2 2], 'replicate'); out_kuwa = zeros(size(noisy)); for i = 3:size(imgP,1)-2 for j = 3:size(imgP,2)-2 block = imgP(i-2:i+2, j-2:j+2); b1 = block(1:3, 1:3); b2 = block(1:3, 3:5); b3 = block(3:5, 1:3); b4 = block(3:5, 3:5); means = [mean(b1(:)) mean(b2(:)) mean(b3(:)) mean(b4(:))]; vars = [var(b1(:)) var(b2(:)) var(b3(:)) var(b4(:))]; [~, idx] = min(vars); out_kuwa(i-2, j-2) = means(idx); end end

这个双重循环在 MATLAB 里效率低，我一般会在代码包里把它替换成nlfilter或blockproc来加速，但循环版本的逻辑最直观。注意padarray的replicate选项，和imfilter一致，避免边界出现异常子块。方差选择还有一个扩展版本，就是在方差里加一个正则项vars + lambda，lambda越大输出越接近全均值，越小越倾向于保留边缘纹理，代码包15-1.m对lambda做了 0、0.01、0.1 三组对比。

2.4 双边滤波与 Lee 滤波：权重函数的设计差异

双边滤波的空间域权重和灰度域权重相乘，灰度差越大权重越小，这使它在平坦区域等效于均值滤波，在边缘处等效于「边缘另一侧不参与平均」。Lee 滤波则基于局部均值和方差做线性最小均方误差估计，它对乘性噪声（如 SAR 图像）有理论最优性，但对高斯噪声的表现不如双边滤波直观。

这里把六种空间域滤波器放在同一张表里对比，代码包里的3-2.m应当还包含六抽头滤波，即六系数 FIR 水平/垂直分离滤波，常用于视频去隔行前后处理，在静态图上可以看作一个 6×1 的平滑核，与均值滤波效果接近但具有更好的频率选择性。空间域这一组跑完后，基本能看出规律：噪声若是脉冲型，中值滤波优先；若是高斯型且边缘保护要求高，直接上双边或 Kuwahara；若不想调参又要稳定输出，超限邻域滤波是性价比最高的选择。

3. 频域与逆问题：Wiener 滤波、约束最小平方和同态滤波的实现边界

3.1 逆滤波为什么会在 MATLAB 里直接崩掉

频域滤波的基本操作是G = F .* H + N，其中H是退化传递函数，N是噪声频谱。最简单的恢复方式是把观测频谱除以H，即逆滤波：

% inverse_filter.m F = fft2(noisy); H = fspecial('motion', 21, 11); % 运动模糊核 Hf = psf2otf(H, size(noisy)); G = F ./ Hf; % 直接相除 out_inv = real(ifft2(G));

这段代码跑完，输出往往是一幅布满亮暗斑点的图。原因是Hf在高频区域接近零，噪声被放大到无穷，人类视觉对高频噪声又极敏感。我在代码包里看到17-4.m处理逆滤波时，采取的补偿措施是对Hf做截断：把绝对值小于某个阈值的频点置为一个固定小量，而不是让它参与除法。这一步牺牲了高频细节，但保住了整体灰度动态范围。

3.2 Wiener 滤波与约束最小平方滤波的 MATLAB 参数

Wiener 滤波在频域的形式是G = conj(H) .* F ./ (abs(H).^2 + 1/SNR)，MATLAB 的deconvwnr把它封装成了三种调用方式：传噪声信号比、传自相关函数、传噪声功率谱。代码包里的标准用法是第三种：

% wiener_deconv.m PSF = fspecial('gaussian', [7 7], 2); wnr1 = deconvwnr(noisy, PSF, 0.01); % 固定 NSR noise_var = 0.01; sig_var = var(noisy(:)); wnr2 = deconvwnr(noisy, PSF, noise_var / sig_var); % 比值形式

deconvwnr的第一个参数是退化图像，第二个是点扩散函数，第三个是信噪比参数。比值形式比固定 NSR 更稳，因为var(noisy(:))包含了信号和噪声的联合方差，当噪声方差超过信号方差时，比值会远大于 1，等效于更强的高频衰减，防止振铃。约束最小平方滤波deconvreg的用法类似，但多了拉格朗日乘子alpha：

% reg_deconv.m [out_reg, reg_arg] = deconvreg(noisy, PSF, 0.4); % alpha=0.4

reg_arg返回算法实际使用的正则化参数，如果输出图像仍有振铃，就把alpha调大一个数量级；如果图像发糊，就把alpha调小。相比 Wiener 滤波，约束最小平方对噪声模型不敏感，适用面更广，代价是参数调节没有解析解，只能靠观察。

3.3 同态滤波：照射反射模型下的动态范围压缩

同态滤波把图像建模为照射分量乘反射分量，取对数后变成加法，再用高通滤波压缩低频照射、保留高频反射：

% homomorphic.m img_log = log(noisy + eps); F = fft2(img_log); [rows, cols] = size(F); u = linspace(-0.5, 0.5, rows)'; v = linspace(-0.5, 0.5, cols); [Dx, Dy] = meshgrid(v, u); D = sqrt(Dx.^2 + Dy.^2); gh = 1.2; gl = 0.4; cutoff = 0.3; H = (gh - gl) .* (1 - exp(-D.^2 / (2*cutoff^2))) + gl; G = F .* ifftshift(H); out_homo = exp(real(ifft2(G)));

gh是高频增益，gl是低频增益，cutoff是高斯高通截止频率。ifftshift这一步容易漏：meshgrid生成的频率原点在左上角，而H的构造假设原点在中心，必须用ifftshift对齐。同态滤波适合光照不均匀的图像，比如一张半边暗半边亮的照片，用它压缩亮度差异后，后续阈值分割会稳定很多。但它对加性噪声没有建模能力，噪声强度高时会把噪声当作反射分量放大，所以代码包里的顺序是先把噪声做一次中值预滤波，再进同态滤波。

4. 反卷积与保边去噪：Richardson-Lucy 迭代、NoLocalMeans 与非线性复扩散

4.1 Richardson-Lucy 的迭代语义和阻尼参数

Richardson-Lucy 最早用于天文图像恢复，它假设噪声服从泊松分布，用期望最大化推导出迭代格式。MATLAB 中的函数名是deconvlucy，它接受的参数远比常见教程里写得多：

% rl_deconv.m PSF = fspecial('gaussian', [9 9], 3); DAMPAR = 0.02; % 阻尼阈值，控制噪声放大 WEIGHT = ones(size(noisy)); % 坏点权重，如坏像素置 0 out_rl = deconvlucy(noisy, PSF, 10, DAMPAR, WEIGHT);

迭代次数默认是 10，这个值对结果影响极大。迭代少则模糊残留，迭代多则出现「振铃 + 斑点噪声」，因为 RL 对高频噪声的放大是随迭代单调增加的。DAMPAR的作用是限制每一步修正量，修正值低于该阈值就不更新；WEIGHT则能屏蔽掉坏像素对周围区域的干扰。代码包14-1.m对迭代次数做了 5、10、20、50 四组对照，结论是泊松噪声场景下 10 次左右性价比最高，高斯噪声下 RL 并不占优，应该考虑 Wiener。

4.2 NoLocalMeans 的搜索窗口与平滑强度

NoLocalMeans（NLM）和双边滤波的区别在于相似度计算的对象：双边滤波逐像素比灰度，NLM 逐块比局部邻域结构。MATLAB 在较新版本里直接提供imnlmfilt：

% nlm_filter.m out_nlm = imnlmfilt(noisy, 'DegreeOfSmoothing', 0.08, ... 'SearchWindowSize', 21, ... 'ComparisonWindowSize', 7);

DegreeOfSmoothing控制指数核的衰减速度，值越大平滑越强，一般取噪声标准差的 0.5 到 1.5 倍；SearchWindowSize是搜索窗，越大越能找到更多相似块，但计算量平方增长；ComparisonWindowSize是块大小。NLM 对纹理区域的保护能力明显强于双边滤波，主要代价是速度：21×21 的搜索窗在大图上跑几十秒很正常。代码包给出的优化思路是先降采样加速，最后对结果做一次边缘保持上采样。

4.3 非线性复扩散滤波的迭代格式

非线性复扩散滤波是 Beltrami 流之外的另一个变分思路：把扩散系数从实数扩展到复数，虚部提供边缘增强。常见半隐式离散格式如下：

% complex_diffusion.m phi = noisy + 1i * 0.1 * noisy; % 初值加入虚部 lambda = 0.15; dt = 0.2; alpha = 0.2; % kappa 控制边缘敏感度 kappa = 2.0; for iter = 1:15 g = 1 ./ (1 + abs(imag(phi)) / kappa); % 边缘停止函数 lap = del2(phi); % 拉普拉斯算子 phi = phi + dt * (lambda * del2(g .* real(phi)) + 1i * alpha * lap); end out_cdiff = real(phi);

这里的del2在 MATLAB 中计算拉普拉斯时已包含 0.25 的归一化系数，所以dt的取值范围与显式差分不同，一般 0.1 到 0.3 之间稳定。虚部imag(phi)在迭代中积累的是二阶导数信息，边缘处虚部较大，对应的g变小，扩散被抑制；平坦区域虚部接近零，扩散正常进行，最终效果是在去噪的同时增强边缘。

5. Beltrami 流与小波阈值去噪：两类非线性方法的 MATLAB 实现路径

5.1 Beltrami 流：把它理解为「图像即曲面」的几何流

Beltrami 流把灰度图像看作嵌入高维空间中的二维流形，去噪过程就是让这个流形按面积最小化方向演化。相比 Perona-Malik 各向异性扩散，Beltrami 流的扩散张量由图像梯度构造，保边能力更好且对参数不敏感。代码包里没有现成的函数，通常自己写显式迭代：

% beltrami_flow.m u = noisy; dt = 0.05; iter = 20; sigma = 1.0; gauss = fspecial('gaussian', [5 5], sigma); for t = 1:iter ux = imfilter(u, [-1 0 1], 'replicate') / 2; uy = imfilter(u, [-1; 0; 1], 'replicate') / 2; uxx = imfilter(u, [1 -2 1], 'replicate'); uyy = imfilter(u, [1; -2; 1], 'replicate'); uxy = imfilter(ux, [-1; 0; 1], 'replicate') / 2; W2 = 1 + ux.^2 + uy.^2; % Beltrami 流离散格式：1/sqrt(W2) * (拉普拉斯 - 二次项) lap = uxx .* (1 + uy.^2) - 2 .* ux .* uy .* uxy + uyy .* (1 + ux.^2); u = u + dt .* lap ./ W2.^2; end out_beltrami = u;

W2是度量张量的行列式相关量，lap ./ W2.^2实现了非线性扩散。dt超过 0.1 时数值不稳定，表现为迭代若干步后图像出现棋盘格；检测方法很简单，打印每次迭代前后的max(abs(u(:) - u_old(:)))，这个值如果跳变到 1e-2 量级，就说明步长过大。Beltrami 流的优点是迭代次数增加时图像细节不会像均值滤波那样持续模糊，而是在保边和平滑之间收敛到一个平衡点。

5.2 小波阈值去噪的三层参数

小波去噪是分解、阈值、重构三步。MATLAB 的wavedec2负责分解，wthresh负责阈值：

% wavelet_denoise.m [coefs, books] = wavedec2(noisy, 3, 'db4'); % 三层分解，db4 小波 [thr, sorh, keepapp] = ddencmp('den', 'wv', noisy); denoised = wdencmp('gbl', coefs, books, 'db4', 3, thr, sorh, keepapp); % 自定义硬阈值写法 c_hard = coefs; c_hard(abs(c_hard) < thr) = 0; out_wav = waverec2(c_hard, books, 'db4');

ddencmp自动估计全局阈值，sorh为's'表示软阈值，'h'表示硬阈值。软阈值把所有系数往零收缩，视觉更平滑但没有硬阈值锐利；硬阈值保留大于阈值的系数，细节好但可能出现小幅震荡。第三层分解后，近似分量keepapp不要做阈值处理，否则图像整体亮度会被压低。小波去噪在混合噪声场景下的表现比空间域滤波稳定，因为它天然把信号和噪声按频带分开，高频层阈值处理对低频分量的扰动极小。

5.3 小波与 Beltrami 在代码包里的搭配顺序

代码包的习惯是先用小波去掉高频高斯噪声，再用 Beltrami 流处理剩余的结构性伪影。这个顺序有实际依据：小波阈值对孤立噪声点有效，但对纹理边缘处的噪声块处理不干净，而 Beltrami 流恰好擅长修复这类边缘附近的不规则噪声。反过来顺序不行，先跑 Beltrami 会把高频细节磨掉一部分，小波再去噪时可能把真正的纹理误判为噪声。

6. 用一份测试脚本快速横向评估滤波结果

代码包里Run_Draw.m的职责是统一出图，但只出图很难量化对比。我在工程实践中的做法是补一个eval_filters.m，用 PSNR 和 SSIM 自动打分：

% eval_filters.m results = struct(); results.avg = filter_eval(noisy, out_avg, img); results.med = filter_eval(noisy, out_med, img); results.nlm = filter_eval(noisy, out_nlm, img); results.rl = filter_eval(noisy, out_rl, img); ssim_table = struct2table(results, 'RowNames', {'PSNR', 'SSIM'}); function scores = filter_eval(original, noisy, clean) scores = [psnr(noisy, clean); ssim(noisy, clean)]; end

真正调参时只看两个数不够，我一般会在15-1.m的基础上加一个残差图：imshowpair(noisy, clean, 'diff')，红色区域代表滤波过度，绿色区域代表残留噪声。若红色集中在边缘，说明滤波器保边太强，适当调大 NLM 的DegreeOfSmoothing或小波的软阈值；若绿色均匀散布，说明阈值偏低，把ddencmp返回的thr乘以 1.2 再跑一轮。最后把每组实验的 PSNR、SSIM、运行时间和参数量填进表格，哪个滤波器值得继续调参，哪组参数已经逼近上限，直接看哪一行震荡最大。

本文还有配套的精品资源，点击获取