ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

横条纹噪声怎么去掉?四阶巴特沃斯陷波滤波器原理与MATLAB实现

横条纹噪声怎么去掉?四阶巴特沃斯陷波滤波器原理与MATLAB实现 搞图像处理的朋友十有八九都撞见过横条纹老照片扫描出来一层明暗交替的纹路工业相机采回的图像带着均匀横向干扰监控视频截个图也被扫出几道横杠。这种噪声压在图上特别碍眼空域里换中值滤波、均值滤波轮番试要么把条纹越搞越糊要么压根压不下去。我去年处理一批归档扫描件就踩了这个坑最后是转到频域用四阶巴特沃斯陷波滤波器把横条纹对应的频率成分定向压掉几分钟就搞定。这篇就把整个思路、参数原理和一套可直接复现的MATLAB代码完整分享出来给正被横条纹折磨的朋友一个能抄作业的解法。1. 横条纹噪声为什么难处理1.1 空域滤波的天然局限先说结论横条纹不是普通的随机噪声它是周期性噪声空域滤波器处理周期性噪声天生吃亏。中值滤波、均值滤波这类操作本质上是拿一个局部窗口在图像上滑动用窗口内像素的统计值替换中心像素。这个逻辑对高斯噪声、椒盐噪声这类随机分布的噪声有效因为噪声点周围的正常像素能“拉”它一把。但横条纹是全局性的、有规律的明暗变化任何局部窗口里看到的都只是条纹的一小段窗口根本没有足够的上下文判断哪些是纹理、哪些是条纹。结果就是窗口开小了压不住条纹开大了把图像的细节一并抹掉。我刚开始处理扫描件的时候也是这个思路先上3x3中值滤波发现条纹纹丝不动换成5x5、7x7条纹是淡了一点但画面上文字的边缘全糊了人脸皮肤也变成了塑料质感。折腾一下午效果图拿给同事看对方第一句话就是“这图怎么像磨皮美颜过”。这就是空域滤波处理周期性噪声的典型下场。换一个角度来看横条纹之所以难弄是因为它和图像内容在空间上是叠加在一起的你在像素域里看到的是“原始画面条纹”的混合体。只要还是在空间域里操作就永远要在“去条纹”和“保细节”之间艰难权衡而这两者在像素级上是纠缠在一起的。1.2 换到频域看横条纹的真面目傅里叶变换的思路完全不一样——它不直接在像素上动手而是把整张图拆解成不同频率、不同方向的正弦波分量。每张图像都可以看成无数个正弦波的叠加低频分量对应图像中大面积的亮度变化高频分量对应边缘、纹理这些细节。把图像变换到频域之后噪声和细节就不再是空间上纠缠的混合体而是分布在不同频率坐标上的独立成分。最关键的是周期性条纹在频域里的表现极其有辨识度。横条纹在空间上是水平方向延伸、垂直方向周期变化的也就是说它只在垂直方向上存在频率变化。对应到频谱上它不会像随机噪声那样铺满整个高频区域而是集中在垂直频率轴附近形成一对关于中心对称的明亮斑点。这个亮斑的位置直接反映了条纹的周期条纹越密亮斑离频谱中心越远条纹越疏亮斑离中心越近。看到频谱上这几个孤立的亮点思路一下就清晰了——既然横条纹的“能量”集中在这么几个点上那只要把这几个点附近的频率成分压掉图像里的条纹自然就消失了而其他频率成分几乎不受影响细节也就保住了。这就是频域去噪的核心思想不是把整幅图模糊掉而是精准打击。1.3 为什么偏偏选巴特沃斯频域去噪的方案不止一种理想陷波滤波器最简单直接在频谱上把特定区域的数值置零就行。但理想滤波器的频率响应是硬跳变的通带和阻带之间没有任何过渡这会在滤波后产生明显的振铃效应——图像边缘附近出现一圈一圈的水波纹这在视觉上甚至比原来的横条纹还难受。高斯滤波器倒是平滑了但它的过渡带比较宽选择性不够强。想要在“压得住噪声”和“保得住细节”之间找到平衡点巴特沃斯滤波器是最合适的选择。它的频率响应在通带内非常平坦过渡带又是平滑下降的而且可以通过调节阶数控制过渡带陡峭程度。阶数越高过渡带越窄越接近理想滤波器但振铃风险也随之上升。四阶是我在实际项目里反复对比后的一个经验值。四阶巴特沃斯的过渡带足够窄能把噪声频率点和周围的正常频率分得很开同时又不会像高阶那样在图像边缘附近激起明显的振铃。这个平衡点对于绝大多数的横条纹去除场景都成立这也是标题里专门强调“四阶”的原因。2. 陷波滤波器设计与参数原理2.1 从带阻到陷波不是所有频带都要切如果理解到横条纹在频域里是几个孤立的亮斑那最合理的处理方案自然不是把整个圆周切掉而是只针对这几个亮斑做定向抑制。带阻滤波器是以频谱中心为圆心、以某个频率为半径把整个圆环上的频率全部衰减。但图像中的纹理细节分布在整个频率平面里很可能在这个圆环的其它方向上有大量有用的高频信息。用带阻滤波器去处理横条纹属于“杀敌一千自损八百”。陷波滤波器Notch Filter就是专门为这种情况设计的。它的作用范围不是整个圆环而是频域中一个很小的局部区域——准确地说是以某个特定频率点为中心的一个小邻域。对于横条纹来说就是在频谱垂直方向上的两个对称亮斑附近各放一个窄带衰减区把条纹对应的频率成分压制下去其余方向的频率成分原样保留。可以说陷波滤波器是“精准手术刀”而带阻滤波器是“大范围化疗”。这也是我在项目里首选陷波方案的根本原因。横条纹的频谱特征本身就是点状的用点状的滤波器去匹配才能做到既去噪又不伤画质。2.2 四阶巴特沃斯陷波的数学表达陷波滤波器是围绕一个中心频率点设计的它的核心是计算频域中每个像素到陷波中心的距离。假设陷波中心坐标为 (u0, v0)频域中任意一点 (u, v) 到该中心的距离为Dk sqrt((u - u0)^2 (v - v0)^2)在这个距离定义的基础上四阶巴特沃斯陷波滤波器的频率响应可以写成Hk 1 - 1 / (1 (Dk / R)^(2n))其中 R 是陷波半径控制阻带的宽度n 是阶数取4。当 Dk 等于 0 时也就是频点正好落在陷波中心上Hk 0该频率被完全抑制当 Dk 等于 R 时Hk 0.5也就是半功率点此时过渡带开始起作用当 Dk 远大于 R 时Hk 趋向于 1该频率几乎不受影响。如果条纹对应的亮斑不止一对比如图像中存在多重横条纹周期那就分别对每个亮斑构造一个陷波滤波器然后把它们逐点相乘得到一个总的滤波器掩膜H_total H1 · H2 · H3 · …这种串联方式的物理含义非常直观经过第一个陷波某个噪声频率被压掉经过第二个陷波另一个噪声频率被压掉。因为每个陷波滤波器在远离自身中心时都接近于1所以它们之间的相互干扰极小串联不会对非目标频率造成额外损失。2.3 三个关键参数少一个都不行陷波滤波器看起来是几行公式但真正用起来无非是调三个参数陷波中心位置 (u0, v0)、陷波半径 R、阶数 n。这三个参数各有各的脾气。陷波中心位置决定“打哪里”。这是整个滤波过程里最关键的参数。定位不准滤波器的衰减区域就会偏离真正的噪声频率点要么条纹残存要么把不该压的频率压掉了。实际操作中这个位置不是猜出来的而是从幅度谱上直接读出来的——在MATLAB的figure窗口里用Data Cursor点一下亮斑坐标自然就出来了。陷波半径决定“打多大范围”。R取值太小阻带覆盖不了噪声频率的扩展范围条纹会残留R取值太大正常频率成分被连坐图像细节损失明显。我一般从 R 10 左右开始试观察效果再微调。阶数决定“过渡带多陡”。n1时过渡带非常平缓带外衰减慢虽然不会振铃但选择性太差n2到n4逐渐变陡n6以上过渡带已经很陡接近理想的矩形响应但振铃风险迅速上升。横条纹场景四阶是一个省心的默认值。3. MATLAB实战完整代码与逐段讲解3.1 准备工作与整体流程运行环境方面这套代码在MATLAB R2016b以上版本都可以直接跑不需要任何额外工具箱——fft2、fftshift、ifft2这些都是MATLAB的基础函数图像处理工具箱在读取图像时会用到如果只是用内置图像做测试基本环境就够用。整体流程分成六步读图与灰度化、傅里叶变换与中心化、构造频率坐标网格、生成陷波滤波器掩膜、频域滤波与逆变换、结果可视化与保存。这个流程是频域图像处理的通用框架换一种噪声类型只需要改掩膜生成那一部分其余步骤可以原封不动复用。3.2 图像读取与傅里叶变换先读入图像。如果输入是彩色图先用 rgb2gray 转成灰度图因为横条纹的频率特征在亮度通道上最明显。然后转成 double 类型fft2 对 double 类型的输入处理精度最高。傅里叶变换本身没什么可说的一个 fft2 调用就完成了。但有一个非常关键的细节fft2 计算出来的频谱零频分量在左上角直接看幅度谱是一片以左上角为最亮点的图高频分布在四周这不符合我们的操作习惯。所以第二步用 fftshift 把零频移到矩阵中心让频谱呈现“中间亮、四周暗”的直观形态所有以中心为基准的频率坐标计算都以这个中心化后的频谱为准。3.3 构造频率坐标网格这一步是整个实现里最容易懵的地方。图像的尺寸是 M 行 N 列对应的频率坐标网格也是 M 行 N 列其中每个点代表一个频率值。u 的方向对应图像的列方向水平方向v 的方向对应图像的行方向垂直方向。用 meshgrid 生成两个矩阵 U 和 VU 给出每个像素点处的水平频率坐标V 给出每个像素点处的垂直频率坐标。坐标值的范围是从 -floor(N/2) 到 N-1-floor(N/2)这个范围是跟 fftshift 后的频谱对齐的。以一幅大小为 512x512 的图像为例频谱中心点对应坐标 (0,0)最左侧的频率坐标是 -256最右侧是 255。后面算陷波中心和像素距离全部在这个坐标系里进行。3.4 生成陷波掩膜与频率响应分析滤波器的掩膜和图像一样是一个矩阵。初始化一个全1矩阵 H然后针对每一对陷波中心计算频域中每个点到该中心的距离 Dk按照四阶巴特沃斯公式算出陷波响应 Hk再把 H 与 Hk 逐点相乘。陷波中心怎么填对于横条纹来说亮斑出现在垂直频率轴上也就是 u0 0v0 ±某个值。绝对值的大小取决于条纹的周期。举个例子如果条纹的周期是16个像素在512像素高的图像里亮斑大致出现在 v0 ≈ ±32 的位置。这个数值不用精确推算直接从幅度谱上点一下就能拿到。3.5 完整可运行代码%% % MATLAB 四阶巴特沃斯陷波滤波 去除图像横条纹 % 适用版本R2016b 及以上 %% clear; clc; close all; %% 1. 读入图像并灰度化 img imread(stripes_demo.png); % 替换成你自己的图像路径 if size(img, 3) 3 gray rgb2gray(img); else gray img; end A double(gray); [M, N] size(A); %% 2. 傅里叶变换与中心化 F fft2(A); Fc fftshift(F); % 零频移到中心 S log(1 abs(Fc)); % 幅度谱便于观察 %% 3. 构造频率坐标网格 u (0:N-1) - floor(N/2); % 水平频率分量 v (0:M-1) - floor(M/2); % 垂直频率分量 [U, V] meshgrid(u, v); %% 4. 设置陷波参数核心根据你的频谱亮斑位置修改 notch_centers [ 0, 32; % 亮斑A先点频谱图确认精确坐标 0, -32; % 亮斑B关于中心对称 ]; notch_radius 12; % 陷波半径控制阻带宽度 n 4; % 四阶巴特沃斯 %% 5. 生成四阶巴特沃斯陷波滤波器掩膜 H ones(M, N); for k 1:size(notch_centers, 1) u0 notch_centers(k, 1); v0 notch_centers(k, 2); Dk sqrt((U - u0).^2 (V - v0).^2); % 四阶巴特沃斯陷波响应Dk0时衰减到0远处趋近1 Hk 1 - 1 ./ (1 (Dk / notch_radius).^(2*n)); H H .* Hk; end %% 6. 频域滤波并还原到空间域 Gc Fc .* H; G ifftshift(Gc); % 频谱中心移回左上角 g real(ifft2(G)); % 逆傅里叶变换取实部 g max(0, min(255, g)); % 防止像素值越界 %% 7. 可视化原始图、幅度谱、滤波器掩膜、去噪结果 figure(Name, 四阶巴特沃斯陷波滤波去横条纹, NumberTitle, off); subplot(2, 2, 1); imshow(uint8(A)); title(原始图像); subplot(2, 2, 2); imshow(S, []); title(频域幅度谱); subplot(2, 2, 3); imshow(H, []); title(滤波器掩膜); subplot(2, 2, 4); imshow(uint8(g)); title(去噪结果); %% 8. 保存结果 imwrite(uint8(g), denoised_result.png);把代码复制进MATLAB改一下图像路径和陷波中心参数就能直接跑通。注意 variables 区里确认 M、N 的值没有异常如果图像尺寸不是 2 的幂次完全不用慌fft2 处理任意尺寸都行坐标网格也是自动适配的。4. 参数怎么调从频谱图到完美去噪4.1 第一步用幅度谱锁定条纹频率直接跑代码最有可能出现的情况是条纹去不掉或者图像糊了根源几乎都在陷波中心定位不准。这个参数不能拍脑袋必须从幅度谱上读。运行代码到可视化那一步中间的 subplot(2,2,2) 显示的就是幅度谱。横条纹对应的亮斑一定出现在垂直方向大致在中间竖线上而且往往是两两对称的一对或几对亮点。用MATLAB figure窗口工具栏里的 Data Cursor 工具直接点一下亮斑最亮的那个点读取坐标值。注意这个坐标是像素坐标但因为我们做了 fftshift数值中心的坐标就是 (0,0)所以读出来的横坐标和纵坐标可以直接作为陷波中心的一个坐标对。在实际项目里我遇到过一次比较刁钻的情况亮斑看起来是一个点放大后发现它其实横跨了两三个像素而且垂直方向上还有轻微的弥散。这种情况下只用一个陷波滤波器效果不彻底可以在主亮斑的上下各加一个半径更小的陷波比如第三个陷波中心设为亮斑坐标周围偏移1到2个像素的点用三个陷波叠加覆盖整个噪声能量扩散区域。4.2 第二步确定陷波半径陷波半径 R 直接影响滤波作用范围。这个参数的调节规律其实很清晰R 越小对频率的选择性越精准图像细节保留越好但一旦 R 小于噪声频率扩散范围条纹就会残留一部分R 越大去条纹越彻底但周围正常频率成分被压掉的风险也越大。我的经验是从亮斑尺寸出发估算。在幅度谱中把亮斑区域放大观察亮斑直径大约覆盖多少个像素。如果亮斑看起来直径只有3到4个像素那 R 取5到8就够如果亮斑扩散到了十几个像素R 就需要12到20。拿不准的时候先取 R10 跑一遍看结果再往两个方向微调每次增减2到3个像素。另外提醒一句陷波半径不是越大越好。R 超过亮斑尺寸太多时图像里一些正常的周期性纹理比如建筑物立面的窗户阵列、织物纹理会被一并压掉画面上会出现奇怪的“洗干净”感细节像被橡皮擦擦过一样。4.3 第三步用阶数控制过渡带陡峭度阶数 n 决定的是阻带与通带之间过渡带的陡峭程度。n 越小过渡带越宽陷波作用越“温和”但压制的力度也减弱n 越大阻带越窄越深选择性越强但图像边缘处的振铃也越明显。四阶对大多数横条纹场景都是好选择。如果发现图像在去噪后边缘出现淡淡的涟漪状条纹把 n 降到2或者3通常就能解决反过来如果条纹虽然大部分被消灭了但还是隐隐约约有残留而且你确定陷波中心和半径都没问题那可以试试把 n 提到5甚至6收窄过渡带以增强压制效果。阶数和陷波半径是联动的。R 取得大、n 取得也大阻带范围就会变成一个又深又宽的“陨石坑”对图像内容的影响非常显著。我建议调整的时候先固定 n4把 R 调到合适位置然后再根据振铃情况微调 n。两个参数同时猜来猜去很容易把结果搞得越来越糟。4.4 用评价指标验证去噪效果视觉判断之外量化指标也值得跑一下。PSNR峰值信噪比和 SSIM结构相似性指标是最常用的两个图像质量评价指标如果手里有干净的原始图像做参考计算前后对比最直观% 假设 ref 是干净参考图g 是去噪结果 ref double(ref); g double(g); mse mean((ref(:) - g(:)).^2); psnr_val 10 * log10(255^2 / mse); mu_ref mean(ref(:)); mu_g mean(g(:)); sigma_ref var(ref(:)); sigma_g var(g(:)); sigma_cross mean((ref(:)-mu_ref) .* (g(:)-mu_g)); ssim_val ((2*mu_ref*mu_g 6.5025) * (2*sigma_cross 58.5225)) / ... ((mu_ref^2 mu_g^2 6.5025) * (sigma_ref sigma_g 58.5225));没有干净参考图的时候可以自己造一条“半参考”曲线取图像中包含大面积平坦区域的部分计算该区域在滤波前后的局部方差方差下降得多说明平滑力度大但如果下降得过猛说明细节损失也大。结合视觉判断比单独看一个指标靠谱得多。实际项目里我一般以目视为准指标只做辅助参考。5. 实战中的坑与对应解法5.1 去噪后出现水波纹状振铃这个坑最常见也最容易让人误判代码写错了。现象是条纹确实没了但图像边缘附近多出一圈一圈像涟漪一样的明暗变化尤其在亮度反差大的地方特别明显。根因是滤波器掩膜过渡带过陡等效于在频域里对频率成分做了“硬切”逆变换时在空间域的强边缘位置激起振铃。解法很简单把阶数 n 从4降为2或3让过渡带更平缓一些如果振铃还在同时把陷波半径 R 稍微调大1到2个像素让压制力度分散一点不要集中在某个频率点上。另一个容易被忽略的原因是陷波中心偏离真实亮斑中心。滤波器衰减区域没有精确命中噪声点导致噪声频率只被压掉一半剩下的一半在空间域里形成了频率较低的余波看起来也是条纹但间距跟原来不一样很容易被误判为振铃。5.2 条纹没去干净剩下的还很顽固如果第一次滤波后条纹还在先别急着把 R 继续加大。先回幅度谱确认亮斑的精确位置有没有选偏用 Data Cursor 读取最亮点坐标跟 notch_centers 里的值逐位核对。横条纹的亮斑通常对称分布在垂直轴上如果只填了一个非对称的坐标滤波效果会大打折扣。还有一种情况是条纹本身包含多个频率成分。幅度谱里如果能看到垂直方向上是几个离散的亮斑排成一列那就不是一个陷波能解决的问题。这时需要为每个亮斑分别设置一对陷波中心全部填入 notch_centers 数组。代码里的循环会自动为每个陷波中心生成一个滤波器并相乘不需要改主流程。另外某些图像的横条纹不是严格的水平方向可能带有微小角度幅度谱里的亮斑就会稍微偏离垂直轴。这种情况下陷波中心要跟着偏一点比如本来是 (0, 32)实际可能是 (1, 32) 或 (-1, 32)。逐像素核对坐标是这种问题最稳妥的排查方式。5.3 图像细节被连坐整体变肉这种现象通常是陷波半径 R 偏大。阻带覆盖了噪声频率周边的正常频率成分尤其是图像里原本就有一些跟条纹频率接近的纹理比如细密的布料纹理、百叶窗、栅栏之类很容易被误伤。处理思路有两条一是缩小 R每次减少2个像素观察细节恢复情况二是增加陷波中心数量、缩小每个陷波的 R用多个小陷波替代一个大陷波只精确覆盖亮斑本身。比如原来一个 R18 的大陷波可以拆成 R8 的三个小陷波分别对准亮斑中心和上下边缘的两三个像素覆盖面相同但对周围频率的影响小得多。5.4 彩色图像的横条纹怎么处理彩色图的横条纹一般也发生在亮度通道上色度通道受污染程度较轻甚至没有。如果在RGB三个通道上分别滤波一来计算量大二来三个通道之间的滤波误差容易产生偏色。推荐的做法是把RGB图像转换到YCbCr颜色空间Y通道是亮度Cb和Cr是色度。只对Y通道做陷波滤波Cb和Cr原样保留然后再转回RGB。这样处理速度更快去噪效果也自然不会出现色偏。MATLAB里用 rgb2ycbcr 和 ycbcr2rgb 两个函数就能轻松完成转换。5.5 批量处理时如何自动锁定条纹频率单张图手动调参没问题但遇到几十张同批次图像每张都打开频谱图点坐标就太折磨了。同批次图像的横条纹通常来自同一个采集设备或同一套扫描参数条纹频率几乎是一致的所以一种非常实用的做法是先取其中一张特征明显的图像手动确定陷波中心坐标然后把这组参数固定下来写一个for循环批量处理整个文件夹。如果不同图像之间的条纹频率有微小漂移可以在循环里加一步自动校正处理每张图时在幅度谱垂直方向中心线附近搜索局部极大值找到亮点位置后将其作为该张图的陷波中心。搜索范围限制在垂直轴上 u ∈ [-2, 2] 的窄带内这样既不会误检到其他方向的特征又能适应条纹频率的微小变化。这一步用 findpeaks 或者简单的 max 函数都能实现代码量不大但省事非常多。6. 写在最后经验与建议整理这套流程用到的核心思路说白了就是一句话去噪最关键的不是“猛力滤波”而是“精准定位”。空域里折腾半天解决不了的问题换到频域看一眼就明白了——条纹就是几个点把点按掉就干净了。这也提醒我处理图像问题时先分析噪声在频域里长什么样再决定用什么滤波器比拿到图像直接套滤波函数要高效得多。参数调节上我个人的习惯是固定一个变量其他变量依次搜索。先固定 n4调 RR 基本合适了再调 n 治振铃最后回到陷波中心坐标做精细校准。这套流程虽然土但不会把自己绕晕。最后再分享一个实用小技巧把最终确定的滤波器掩膜 H 保存成 .mat 文件下次遇到同类型图像直接 load 进来用连频谱分析这一步都能省了。由于它和图像的尺寸绑定记得先确认尺寸一致。这个方法我用了大半年处理同类扫描件的效率比最初手动调参翻了一倍不止。如果你手里的图像不止横条纹这一种噪声把这个陷波思路和其它滤波手段叠加使用也能组合出很多灵活的方案。
返回列表