本节书摘来自异步社区《精通Matlab数字图像处理与识别》一书中的第6章，第6.7节，作者 张铮 , 倪红霞 , 苑春苗 , 杨立红，更多章节内容可以访问云栖社区“异步社区”公众号查看
6.7 Matlab综合案例——利用频域滤波消除周期噪声
精通Matlab数字图像处理与识别
6.5~6.6节介绍了几种典型的频域滤波器，实现了频域下的低通和高通滤波，它们均可在空域下采用平滑和锐化算子实现。而本节准备给出一个特别适合在频域中完成的滤波案例，即利用频域带阻滤波器消除图像中的周期噪声。下面就来看一下这个在空域中几乎不可能完成的任务是如何在频域中实现的。
6.7.1 频域带阻滤波器
顾名思义，所谓“带阻”就是阻止频谱中某一频带范围的分量通过，其他频率成份则不受影响。常见的带阻滤波器有理想带阻滤波器和高斯带阻滤波器。
1．理想带阻滤波器
理想带阻滤波器的表达式为
式中：D 0是阻塞频带中心频率到频率原点的距离；W是阻塞频带宽度；D是（u,v）点到频率原点的距离。于是，理想带阻滤波器的频域特性曲面如图6.28所示。
2．高斯带阻滤波器
本案例中使用了高斯带阻滤波器，下面直接给出高斯带阻滤波器的表达式。
式中：D 0是阻塞频带中心频率到频率原点的距离；W是阻塞频带宽度；D是（u,v）点到频率原点的距离。于是，二维高斯带阻滤波器的频域特性曲面如图6.29所示。
3．高斯带阻滤波器的Matlab实现
根据高斯带阻滤波器的定义式（6-60），可以编写高斯带阻滤波器的生成函数如下，该函数被封装在金羽图书论坛（http://bbs.book95.com）的“金羽图书与答疑”板块与本书同名的主题帖子附件的“chapter6/code”目录下的imgaussfbrf.m文件中。
function out = imgaussfbrf(I, freq, width)
% imgaussfbrf函数 构造频域高斯带阻滤波器
% I参数 输入的灰度图像
% freq参数 阻带中心频率
% width参数 阻带宽度
[M,N] = size(I);
out = ones(M,N);
for i=1∶M
for j=1∶N
out(i,j) = 1-exp(-0.5*((((i-M/2)^2+(j-N/2)^2)-freq^2)/(sqrt(i.^2+j.^2)*width))^2);
end
end
6.7.2 带阻滤波消除周期噪声
带阻滤波器常用于处理含有周期性噪声的图像。周期性噪声可能由多种因素引入，例如图像获取系统中的电子元件等。本案例中我们人为地生成了一幅带有周期噪声的图像，而后通过观察分析其频谱特征，选择了合适的高斯带阻滤波器进行频域滤波。
1．得到周期噪声图像
通常可以使用正弦平面波来描绘周期性噪声。如下程序为Matlab示例图片pout.tif增加周期性噪声。
O = imread('pout.tif'); %读入原图像
[M,N] = size(O);
I = O;
for i=1:M
for j=1:N
I(i,j)=I(i,j)+20*sin(20*i)+20*sin(20*j); %添加周期噪声
end
end
subplot(1,2,1);
imshow(O);
title('Source');
subplot(1,2,2);
imshow(I);
title('Added Noise');
添加周期性噪声前后的区别如图6.30所示。
2．频谱分析
使用高斯带阻滤波器时，首先需要对欲处理的图像的频谱有一个了解。下面的命令得到了两幅图像的频谱。
i_f=fft2(I);
i_f=fftshift(i_f);
i_f=abs(i_f);
i_f=log(1+i_f);
o_f=fft2(O);
o_f=fftshift(o_f);
o_f=abs(o_f);
o_f=log(1+o_f);
figure(1);
imshow(o_f, [ ]); %得到图6.31（a）
title('Source');
figure(2);
imshow(i_f, [ ]); %得到图6.31（b）
title('Added Noise');
程序的运行结果如图6.31所示。
3．带阻滤波
观察图6.31（b），发现周期性图像的傅立叶频谱中出现了两对相对于坐标轴对称的亮点，它们分别对应于图像中水平和竖直方向的正弦噪声。我们构造高斯带阻滤波器的时候就需要考虑尽可能滤除具有这些亮点对应的频率的正弦噪声。注意到这4个点位于以频谱原点为中心，以50为半径的圆周上。因此，设置带阻滤波器中心频率为50，频带宽度为5，如图6.31（c）所示，滤波后的频域效果如图6.31（d）所示。
相应的程序如下。
ff = imgaussfbrf(I, 50, 5); %构造高斯带阻滤波器
figure, imshow(ff, []); %得到图6.31（c）
out = imfreqfilt(I, ff); %带阻滤波
figure, imshow(out, []); %得到图6.31（d）
subplot(1,2,1);
imshow(I);
title('Source');
subplot(1,2,2);
imshow(out);
title('Gauss Filter');
上述程序运行后得到的高斯带阻滤波器最终滤波效果如图6.32所示，我们看到周期噪声被很好地消除，这样的效果在空域中是很难实现的。
