《精通Matlab数字图像处理与识别》一6.4 频域滤波基础

本节书摘来自异步社区《精通Matlab数字图像处理与识别》一书中的第6章，第6.4节，作者 张铮 , 倪红霞 , 苑春苗 , 杨立红，更多章节内容可以访问云栖社区“异步社区”公众号查看

6.4 频域滤波基础

精通Matlab数字图像处理与识别

6.4.1 频域滤波与空域滤波的关系

傅立叶变换可以将图像从空域变换到频域，而傅立叶反变换则可以将图像的频谱逆变换为空域图像，即人可以直接识别的图像。这样一来，我们可以利用空域图像与频谱之间的对应关系，尝试将空域卷积滤波变换为频域滤波，而后再将频域滤波处理后的图像反变换回空间域，从而达到图像增强的目的。这样做的一个最主要的吸引力在于频域滤波的直观性特点，关于这一点稍后将进行详细的阐述。

根据著名的卷积定理：两个二维连续函数在空间域中的卷积可由其相应的两个傅立叶变换乘积的反变换而得到；反之，在频域中的卷积可由在空间域中乘积的傅立叶变换而得到，即

其中，F(u, v)和H(u, v)分别表示f(x, y)和h(x, y)的傅立叶变换，而符号<==>表示傅立叶变换对，即左侧的表达式可通过傅立叶正变换得到右侧的表达式，而右侧的表达式可通过傅立叶反变换得到左侧的表达式。

式（6-47）构成了整个频域滤波的基础，卷积的概念我们曾在第5章空间域滤波中讨论过，而式中的乘积实际上就是两个二维矩阵F(u, v)和H(u, v)对应元素之间的乘积。

6.4.2 频域滤波的基本步骤

根据式（6-47）进行频域滤波通常应遵循以下步骤。

（1）计算原始图像f(x, y)的DFT，得到F(u, v)。

（2）将频谱F(u, v)的零频点移动到频谱图的中心位置。

（3）计算滤波器函数H(u, v)与F(u, v)的乘积G(u, v)。

（4）将频谱G(u, v)的零频点移回到频谱图的左上角位置。

（5）计算第（4）步计算结果的傅立叶反变换g(x, y)。

（6）取g(x, y)的实部作为最终滤波后的结果图像。

由上面的叙述易知，滤波能否取得理想结果的关键取决于频域滤波函数H(u, v)。我们常常称之为滤波器，或滤波器传递函数，因为它在滤波中抑制或滤除了频谱中某些频率的分量，而保留其他的一些频率不受影响。本书中我们只关心其值为实数的滤波器，这样滤波过程中H的每一个实数元素分别乘以F中对于位置的复数元素，从而使F中元素的实部和虚部等比例的变化，不会改变F的相位谱，这种滤波器也因此被称为“零相移”滤波器。这样，最终反变换回空域得到的滤波结果图像g(x, y)理论上也应当为实函数。然而由于计算舍入误差等原因，可能会带有非常小的虚部，通常将虚部直接忽略。

为了更为直观地理解频域滤波与空域滤波之间的对应关系，让我们先来看一个简单的例子。6.2节中曾指出了原点处的傅立叶变换F(0, 0)实际上是图像中全部像素的灰度之和。那么如果我们要从原图像f(x, y)得到一幅像素灰度和为0的空域图像g(x, y)，就可以先将f(x, y)变换到频域F(u, v)，而后令F(0, 0) = 0（在原点移动到中心的频谱中为F(M/2, N/2)），再反变换回去。这个滤波过程相当于计算F(u, v)和如下的H(u, v)之间的乘积。

上式中的H(u, v)对应于平移过的频谱，其原点位于（M/2, N/2）。显然，这里H(u, v)的作用就是将点F(M/2, N/2)置零，而其他位置的F(u, v)保持不变。有兴趣的读者可以自己尝试这个简单的频域滤波过程，反变换之后验证g(x, y)的所有像素灰度之和是否为0。我们将在6.4.2小节详细地探讨一些具有更高实用性的频域滤波器。

6.4.3 频域滤波的Matlab实现

为方便读者在Matlab中进行频域滤波，我们编写了imfreqfilt函数，其用法同空域滤波时使用的imfilter函数类似，调用时需要提供原始图像和与原图像等大的频域滤波器作为参数，函数的输出为经过滤波处理又反变换回空域之后的图像。

注意

通常使用fftshift函数将频谱原点移至图像中心，因此需要构造对应的原点在中心的滤波器，并在滤波之后使用ifftshift函数将原点移回以进行反变换。
频域滤波算法imfreqfilt的完整实现如下，该算法被封装在金羽图书论坛（http://bbs.book95.com）的“金羽图书与答疑”板块与本书同名的主题帖子附件的“chapter6/code”目录下的imfreqfilt.m文件中。

function out = imfreqfilt(I, ff)
% imfreqfilt函数        对灰度图像进行频域滤波
% 参数I                输入的空域图像
% 参数ff                应用的与原图像等大的频域滤镜

if (ndims(I)==3) && (size(I,3)==3)   % RGB图像
    I = rgb2gray(I);
end

if (size(I) ~= size(ff))
    msg1 = sprintf('%s: 滤镜与原图像不等大，检查输入', mfilename);
    msg2 = sprintf('%s: 滤波操作已经取消', mfilename);
    eid = sprintf('Images:%s:ImageSizeNotEqual',mfilename);
    error(eid,'%s %s',msg1,msg2);
end

% 快速傅立叶变换
f = fft2(double(I));

% 移动原点
s = fftshift(f);

% 应用滤镜及反变换
out = s .* ff; %对应元素相乘实现频域滤波
out = ifftshift(out);
out = ifft2(out);

% 求模值
out = abs(out);

% 归一化以便显示
out = out/max(out(:));