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

图像处理与识别是当今计算机科学中一个热门研究方向，其应用广泛，发展前景乐观。Matlab是MathWorks公司开发的一款工程数学计算软件。由于其强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计，以及便捷的与其他程序和语言接口的功能，Matlab 已成为当今国际上科学界最具影响力、最有活力的软件。本书正是紧密结合大量的 Matlab 代码和案例展开对数字图像处理和识别技术的介绍。

让计算机理解所“看”到的东西是一件非常神秘和令人兴奋的事情。但要掌握好数字图像处理与识别技术却并非易事，它的理论性较强、门槛较高，在各个高校中，这门课程多是作为计算机专业研究生的选修课程。要理解该领域的知识，读者需要具有一定的数学基础，除此之外，还涉及信号处理、统计分析、模式识别和机器学习等专业领域知识，因此，令很多人望而却步。

其实“难以理解”的关键在于缺乏必要的先序知识，造成了读者在相关知识上难以跨越的鸿沟。我们在撰写本书过程中，对于可能造成读者理解困难的地方，均给出了必要的先序知识，尽量定性地去描述；对于那些并不一目了然的结论均给出了思路和解释；对于某些非常专业、已经超过本书讨论范围的相关知识也都在全书最后给出了参考文献，供有兴趣的读者进一步学习和研究。

本书的宗旨是，向读者介绍知识的同时，培养读者的思维方法，使读者知其然还要知其所以然，并在解决实际问题中能有自己的想法。

全书共分15章，主要内容介绍如下。

第1～2章介绍了数字图像处理的基础知识和Matlab数字图像编程基础，使读者第一步能够建立起对于数字图像本质的正确认识，了解和掌握必要的术语，并且熟悉本书自始至终需要使用的工具Matlab，重点介绍其数字图像处理工具箱。

第3～4章分别介绍了图像的灰度变换和几何变换。通过灰度变换可以有效改善图像的外观，并在一定程度上实现图像的灰度归一化；几何变换则主要应用在图像的几何归一化和图像校准当中。总体而言，这些内容大多作为图像的前期预处理工作的一部分，是图像处理中相对固定和程式化的内容。

第5～6章分别从空间域和频率域两个角度去考量图像增强的各个主要方面。图像增强作为数字图像处理中相对简单却最具艺术性的领域之一，可理解为根据特定的需要突出一幅图像中的某些信息，同时削弱或去除某些不需要的信息的处理方法。其主要目的是，使处理后的图像对某种特定的应用来说，比原始图像更适用。

第7章小波变换继第6章之后继续在频率域中研究图像。傅立叶变换一直是频率域图像处理的基石，它能用正弦函数之和表示任何分析函数，而小波变换则基于一些有限宽度的基小波，这些小波不仅在频率上是变化的，而且具有有限的持续时间。例如，对于一张乐谱，小波变换不仅能提供要演奏的音符，而且说明了何时演奏等细节信息，但傅立叶变换只提供了音符，局部信息在变换中丢失。

第8章图像复原与图像增强相似，其目的也是改善图像质量。但是图像复原是试图利用退化过程的先验知识使已被退化的图像恢复本来面目，而图像增强是用某种试探的方式改善图像质量，以适应人眼的视觉与心理。引起图像退化的因素包括由光学系统、运动等造成的图像模糊，以及源自电路和光学因素的噪声等。图像复原是基于图像退化的数学模型，复原的方法也建立在比较严格的数学推导上。

第9章是本书中相对独立的一章，以介绍色彩模型之间的相互转换，以及彩色图像处理方面的基本概念和基本方法为主。随着基于互联网的图像处理应用在不断增长，彩色图像处理已经成为一个重要领域。

第10～12章（形态学处理、边缘检测与图像分割、特征提取）是从单纯图像处理向图像识别（机器视觉）的过渡，这一阶段的特点是，输入是图像，输出则是在识别意义上我们感兴趣的图像元素。形态学处理是提取图像元素的有力技术，它在表现和描述形状方面非常有用；分割过程则将一幅图像划分为组成部分或目标对象；研究特征提取则是要将前面提取出来的图像元素或目标对象表示为适合计算机后续处理的数值形式，最终形成能够直接供分类器使用的特征。

第13章在前面知识的基础之上，引出了机器视觉的前导性内容，给出了解决识别问题的一般思路。

第14～15章（人工神经网络、支持向量机）介绍了两种十分强大的分类技术，并介绍了人脸识别技术及经典案例应用。

● 具备必要的数学基础的相关专业的本科生、研究生。

● 工作在图像处理和识别领域一线的工程技术人员。

● 对于数字图像处理和机器视觉感兴趣的并且具备必要预备知识的所有读者。

● 读者应具备一定的数学基础，如高等数学知识、少量的线性代数基本概念加上对于概率理论主要思想的理解（识别部分）。

本书所有Matlab实例的源代码均可在bbs.book95.com网站的“金羽图书与答疑板块”板块与本书同名的主题帖子附件中提供下载。虽然本书中的所有例子都已经在Windows XP、Windows 2003和Windows 7等操作系统下的Matlab2006到MatlabR2011a的各个版本中测试通过，但由于交稿时间要求和笔者水平的局限，也有存在 Bug 的可能，即便正确也很可能存在更加优化的算法或更加合理的程序结构，如发现任何上述问题，请您不吝告知本书的作者（aaron@book95.com），以便我们做出改进。

首先要感谢我的授业恩师——南开大学的白刚教授和天津大学的赵政教授，是他们引导我进入了图像处理与机器视觉的研究领域。同时，他们在我写作过程中的指点和教诲确保了本书内容的严谨。

本书的第7章、第8章由王艳艳和赵国宇编写，在此向他们表示感谢；同时感谢我的好友王艳平提供并调试了许多实例代码；感谢我的同事——闵卫东、陈香凝、任淑霞、姚清爽、王佳欣、孙连坤、孙学梅、张振和王作为等为本书所做出的工作；感谢我昔日的师弟和师妹——王杉、闫丽霞、刘旭、赵国宇和李宏鹏等参与了部分章节的编写和修改；感谢罗小科先生为本书制作了很多插图；感谢我的兄长张钊为本书提供了部分照片；还要感谢徐超、王欣和郭朋博士等为本书的编写提出了很多的宝贵意见和建议。

最后感谢我的妻子马宏、儿子张垚淼以及我的家人，没有你们的鼓励和支持就不会有我的这部作品。

空间域和频率域为我们提供了不同的视角。在空间域中，函数的自变量（x，y）被视为二维空间中的一点，数字图像f（x，y）即为一个定义在二维空间中的矩形区域上的离散函数；换一个角度，如果将f（x，y）视为幅值变化的二维信号，则可以通过某些变换手段（如傅立叶变换、离散余弦变换、沃尔什变换和小波变换等）在频域下对它进行分析。

第5章详细介绍了空间域图像增强的有关知识，紧接着本章就从频域的角度去看待和分析图像增强问题，相信这一定会使您对图像增强的理解更加深刻。

● 傅立叶变换的数学基础

● 快速傅立叶变换

● 频率域图像增强

● 高通滤波器和低通滤波器

● 美女与猫——交换两幅图像的相位谱

● 利用频域滤波消除周期噪声

在很多情况下，频率域滤波和空间域滤波可以视为对于同一个图像增强问题的殊途同归的两种解决方式。而在另外一些情况下，有些增强问题更适合在频域中完成（6.7节），有些则更适合在空域中完成。我们常常根据需要选择是工作在空间域还是频率域，并在必要时在空间域和频率域之间相互转换。

傅立叶变换提供了一种变换到频率域的手段，由于用傅立叶变换表示的函数特征可以完全通过傅立叶反变换进行重建，不丢失任何信息，因此它可以使我们工作在频率域，而在转换回空间域时不丢失任何信息。

要理解傅立叶变换，掌握频率域滤波的思想，必要的数学知识是不能跳过的。为便于理解，我们将尽可能定性地去描述。其实傅立叶变换所必需的数学知识对于一个理工科大学二年级以上的学生来说是很有限的，高等数学中傅立叶级数的知识加上线性代数中基和向量空间的概念就足够了。下面就从一维情况下的傅立叶级数开始进行介绍。

法国数学家傅立叶发现任何周期函数只要满足一定条件（狄利赫里条件）都可以用正弦函数和余弦函数构成的无穷级数，即以不同频率的正弦和余弦函数的加权和来表示，后世称为傅立叶级数。

对于有限定义域的非周期函数，可以对其进行周期拓延，从而使其在整个扩展定义域上为周期函数，从而也可以展开为傅立叶级数。

周期为T的函数 f (t)的三角形式傅立叶级数展开为

，而u=1/T，是函数f(x)的频率。

a_k和 b_k称为傅立叶系数。稍后在学习傅立叶级数的复数形式时还将介绍傅立叶系数的另一种形式。事实上，傅立叶系数正是我们在6.2.2小节傅立叶变换中所关心的对象。

于是，周期函数f(t)就与下面的傅立叶序列产生了一一对应，即

图6.1所示形象地显示出了这种频率分解，左侧的周期函数f（x）可以由右侧函数的加权和来表示，即由不同频率的正弦和余弦函数以不同的系数组合在一起。

原函数f（x）（左），其傅立叶展开为一系列不同频率的正弦、余弦函数的加权和（右）

从数学上已经证明了，傅立叶级数的前N项和是原函数f(t)在给定能量下的最佳逼近。

图6.2为我们展示了对于一个方波信号函数采用不同的N值的逼近情况。随着N的增大，逼近效果越来越好。但同时也注意到，在f (x)的不可导点上，如果只取式（6-1）右边的无穷级数中的有限项之和作为

，那么

在这些点上会有起伏，对于图6.2（a）的方波信号尤为明显，这就是著名的吉布斯现象。

除上面介绍的三角形式外，傅立叶级数还有其他两种常用的表现形式，即余弦形式和复指数形式。借助欧拉公式，上述3种形式可以很方便地进行等价转化，本质上它们都是一样的。

复指数傅立叶级数即我们经常说的傅立叶级数的复数形式，因其具有简洁的形式（只需一个统一的表达式计算傅立叶系数），在进行信号和系统分析时通常更易于使用；而余弦傅立叶级数可使周期信号的幅度谱和相位谱意义更加直观，函数的余弦傅立叶级数展开可以解释为f(x)可以由不同频率和相位的余弦波以不同系数组合在一起来表示，而在三角形式中相位是隐藏在系数A_n和b_n中的。下面主要介绍复指数傅立叶级数，在后面的傅立叶变换中要用到的正是这种形式。关于余弦傅立叶级数的有关知识，感兴趣的读者请参考附录Ⅲ。

傅立叶级数的复指数形式为

由式（6-4）和（6-5）可见，复指数傅立叶级数形式比较简洁，级数和系数都可以采用一个统一的公式计算。有关如何由式（6-1）推导出傅立叶级数复指数形式（6-4）的过程，由于这里我们感兴趣的并非傅立叶级数本身，就不在正文中给出了，详细的内容可参考附录Ⅱ，只要您相信不同的展开形式之间本质上是等价的，并对复指数形式的傅立叶级数展开建立了一个基本的形式上的认识就足以继续阅读和理解后面的内容了。

对于定义域为整个时间轴（−∞< t <∞）的非周期函数 f（t），此时已无法通过周期拓延将其扩展为周期函数，这种情况下就要用到傅立叶变换。

由F(u)我们还可以通过傅立叶反变换获得f（t）。

式（6-6）和式（6-7）即为我们通常所说的傅立叶变换对，6.1 节中提到的函数可以从它的反变换进行重建正是基于上面的傅立叶变换对。

由于傅立叶变换与傅立叶级数涉及两类不同的函数，在很多数字图像处理的书中通常对它们分别进行处理，并没有阐明它们之间存在的密切联系，这给很多初学者带来了困扰，实际上我们不妨认为周期函数的周期可以趋向无穷大，这样可以将傅立叶变换看成是傅立叶级数的推广。

仔细地观察式（6-6）和式（6-7），对比复指数形式的傅立叶级数展开公式式（6-4），注意到在这里傅立叶变换的结果F(u)实际上相当于傅立叶级数展开中的傅立叶系数，而反变换公式式（6-7）则体现出不同频率复指数函数的加权和的形式，相当于复指数形式的傅立叶级数展开公式，只不过这里的频率u变为了连续的，所以加权和采用了积分的形式。这是因为随着作为式（6-5）的积分上下限的T向整个实数定义域扩展，即T→∞,频率u则趋近于du（因为u=1/T），导致原来离散变化的u的连续化。

一维函数 f(x)（其中 x=0, 1, 2 ,… ,M-1）的傅立叶变换的离散形式为

相应的反变换为

由于一维情况下很多性质更为直观，我们更青睐于分析一维离散傅立叶变换，而由此得出的这些结论都可顺利推广至二维。一些有用的性质如下。

● 仔细观察式（6-8）和式（6-9），注意到在频域下变换F(u)也是离散的，且其定义域仍为0～M−1，这是因为F(u)的周期性，即

● 考虑式（6-9）中的系数1/M，在这里该系数被放在反变换之前，实际上它也可以位于式（6-8）的正变换公式中。更一般的情况是只要能够保证正变换与反变换之前的系数乘积为1/M即可。例如，两个公式的系数可以均为

。

● 为了求得每一个F(u)（u=0, 1, 2,…,M-1），都需要全部M个点的 f(x)参与加权求和计算。对于M个u，则总共需要大约M ² 次计算。对于比较大M（在二维情况下对应着比较大的图像），计算代价还是相当可观的，我们会在下一节快速傅立叶变换中来研究如何提高计算效率的问题。

有了之前的基础，下面我们将傅立叶变换及其反变换推广至二维。对于二维连续函数，傅立叶变换为

类似的，其反变换为

在数字图像处理中，我们关心的自然是二维离散函数的傅立叶变换，下面直接给出二维离散傅立叶变换（Discrete Fourier Transform, DFT）公式。

相对于空间域（图像域）的变量x、y，这里的u、v是变换域或者说是频率域变量。同一维中的情况相同，由于频谱的周期性式（6-13）只需对u值（u=0, 1, 2, …,M−1）及 v值（v=0, 1, 2, …N−1）进行计算。同样，系数1/MN的位置并不重要，有时也放在正变换之前，有时则在正变换和反变换前均乘以系数

。

根据式（6-13），频域原点位置的傅立叶变换为

显然，这是f(x, y)各个像素的灰度之和。而如果将系数1/MN放在正变换之前，则F(0, 0)对应于原图像f(x, y)的平均灰度。F(0, 0)有时被称作频谱的直流分量（DC）。

我们之前曾指出了一维函数可以表示为正弦（余弦）函数的加权和形式；类似的，二维函数f(x, y)可以分解为不同频率的二维正弦（余弦）平面波的按比例叠加。图6.3（a）中给出了一幅简单的图像，可将它视为以其灰度值作为幅值的二维函数，如图 6.3（b）所示，根据式（6-13），它可以分解为如图6.3（c）所示的不同频率和方向的正弦（余弦）平面波的按比例叠加（只给出了一部分）。比如图6.3（c）中第一行中间的平面波为sin(Y)，而第二行右面的平面波则为sin(X+2Y)，而第三行最后的一个为sin(2X+2Y)。

下面，我们再来定义傅立叶变换的幅度谱、相位谱以及功率谱。

显然，幅度谱关于原点具有对称性，即|F (−u,−v)|=|F (u, v)|。

通过幅度谱和相位谱，我们可以还原F (u, v)。

● 功率谱（谱密度）

式中：Re(u, v)和Im(u, v)分别为F(u,v)的实部和虚部。

幅度谱又叫频率谱，是图像增强中关心的主要对象，频域下每一点（u,v）的幅度|F(u, v)|可用来表示该频率的正弦（余弦）平面波在叠加中所占的比例，如图6.4所示。幅度谱直接反映频率信息，是频域滤波中的一个主要依据。

图6.4所示幅度谱中的A、B、C、D四点的幅值分别为四周的 4 个正弦平面波在的加权求和中的权值（混合比例）。注意这4个正弦平面波的方向和频率。

相位谱表面上看并不那么直观，但它隐含着实部与虚部之间的某种比例关系，因此与图像结构息息相关。

由于对于和空域等大的频域空间下的每一点（u, v），均可计算一个对应的|F(u, v)|和φ(u, v)，因此可以像显示一幅图像那样显示幅度谱和相位谱。图6.5（b）、（c）分别给出了图6.5（a）中图像的幅度谱和相位谱，获得它们的方法请参考6.3节中傅立叶变换实现的相关内容，关于幅度谱和相位谱的一个非常有趣的例子请参考例6.2。

▲图6.5 circuit.tif幅度谱和相位谱。幅度谱和相位谱都将（0,0）点移到了中心

无论是傅立叶变换、离散余弦变换还是小波变换，其本质都是基的变换。下面首先让我们一起回顾一下线性代数中基和向量空间的相关知识。

在三维欧氏向量空间中，某向量

可以由3个复数{v₁,v₂,v₃}来定义，常常记作

=(v₁,v₂,v₃)，这3个复数与3个正交单位向量{

}相联系。实际上，有序集{v₁,v₂,v₃}表示向量

的3个标量分量，也就是系数；而3个正交单位向量{

}即为该三维欧氏空间的基向量。我们称该空间为这3个基向量所张成的空间，任何该空间中的向量

均可由这3个基向量的线性组合（加权和）表示为

也可以用矩阵的形式来表示该向量。

在上面的叙述中涉及了向量的正交，这是向量代数中一个非常重要的概念。为了说明正交的概念，让我们首先回顾一下向量点积（数量积），两个向量的点积定义为

式中：

，表示向量

的模；θ为向量

和

之间的夹角；上标T表示转置。

此时，如果

=0，则称这两个向量

和

互相正交。由式（6-22）可知，两非零向量正交则cosθ=0 ，说明其夹角为90°（垂直）。

接下来，定义一个向量在另一个向量方向上的投影或分量为

式中：

为向量

单位化后的单位向量，模为1，方向与

相同。式（6-23）说明如果需要得到某向量在给定方向上的分量，只需计算该向量与给定方向单位向量的点积。

图6.6能够帮助我们理解上述内容，图6.6（a）中为一个三维空间中的向量

以及3个单位正交基向量

；图6.6（b）中给出了向量

在

方向的投影

；在图6.6（c）中，根据矢量加法的平行四边形法则，向量

被分解为3个正交基向量

线性组合，显然可以表示为

=(v₁,v₂,v₃)的形式。

将三维向量空间中基与投影的概念推广至N维向量空间。任何一个该空间中的N×1向量均可由N个基向量

,...,

的线性组合来表示，记作

式中：分量v_i为向量

在

方向的投影。

式（6-24）称为对

的重构，式（6-25）称为对

的分解。

而N个单位基向量之间满足两两正交关系，即

尽管上面的向量分解与重构的问题比较基础，但它与傅立叶变换与反变换之间的关系却十分紧密。事实上，它们在形式上有着惊人的相似，唯一不同的是这里的向量空间变成了函数空间，向量

变成了函数f(x)，而基向量

,…,

也相应地变成了基函数。对比式（6-24）～（6-25）和式（6-8）～（6-9）的形式不难看出，式（6-25）的分解过程即相当于傅立叶变换，而式（6-26）的重构过程则恰恰相当于傅立叶反变换。也就是说，相应函数空间中的任意函数均可以由该函数空间中的一组基函数的加权和来表示。观察式（6-8）容易发现，这里的基函数的形式为e⁻ⁱ ^2πux ，我们用下面的等式来表示函数的正交性。

至此，读者应该已经理解了傅立叶变换的实质——基的转换。对于给定函数f(x)，关键是选择合适的基，使得f(x)在这组基下表现出我们需要的特性。当某一组基不满足要求时，就需要通过变换将函数转换到另一组基下表示，方可得到我们需要的函数表示。常用的变换有傅立叶变换（以正弦和余弦函数为基函数）、小波变换（以各种小波函数为基函数）、离散余弦变换以及 Walsh 变换等。实际上，我们在第12章中将指出，特征降维中常用的主成份分析法（K-L变换）本质上也是一种基的转换。

6.2节介绍了离散傅立叶变换（DFT）的原理，但并没有涉及其实现问题，这主要是因为DFT的直接实现效率较低。在工程实践中，我们迫切地需要一种能够快速计算离散傅立叶变换的高效算法，快速傅立叶变换（FFT）便应运而生。本节将给出快速傅立叶变换算法的原理及其实现细节。

之所以提出快速傅立叶变换（FFT）方法，是因为在计算离散域上的傅立叶变换时，对于N点序列，它的DFT变换与反变换对定义为

于是不难发现，计算每个u值对应的F(u)需要N次复数乘法和N-1次复数加法。因此，为了计算长度为N的序列的快速傅里叶变换，共需要执行n² 次复数乘法和N(N-1)次复数加法。而实现1次复数相加至少需要执行2次实数加法，执行1次复数相乘则可能需要至多4次实数乘法和2次实数加法。如果使用这样的算法直接处理图像数据，则运算量会大得惊人，更无法实现实时处理。

然而，离散傅立叶变换的计算实质并没有那么复杂。在离散傅里叶变换的运算中有大量重复运算。上面的变量W_N是一个复变量，但是可以看出它具有一定的周期性，实际上它只有N个独立的值。而这N个值也不是完全相互独立的，它们又具有一定的对称关系。关于变量W_N的周期性和对称性，我们可以做如下总结。

式（6-29）是W矩阵中元素的某些特殊值，而式（6-30）则说明了W矩阵元素的周期性和对称性。利用W的周期性，DFT运算中的某些项就可以合并；而利用W的对称性，则可以仅计算半个W序列。而根据这两点，我们就可以将一个长度为N的序列分解成两个长度为N/2的序列并分别计算DFT，这样就可以节省大量的运算量。我们将在讲述常见的FFT算法后分析节省的运算量。

这正是快速傅立叶变换（Fast Fourier Transform，FFT）的基本思路——通过将较长的序列转换成相对短得多的序列来大大减少运算量。

目前流行的大多数成熟的 FFT 算法的基本思路大致可以分为两大类，一类是按时间抽取的快速傅立叶算法（Decimation In Time, DIT-FFT），另一类是按频率抽取的快速傅立叶算法（Decimation In Freqency, DIF-FFT）。这两种算法思路的基本区别如下。

按时间抽取的FFT算法是基于将输入序列f(x)分解（抽取）成较短的序列，然后从这些序列的DFT 中求得输入序列的 F(u)的方法。由于抽取后的较短序列仍然可分，所以最终仅仅需要计算一个很短的序列的 DFT。在这种算法中，我们主要关注的是当序列的长度是 2 的整数次幂时，如何能够高效地进行抽取和运算的方法。

而按频率抽取的FFT算法是基于将输出序列F(u)分解（抽取）成较短的序列，并且从f(x)计算这些分解后的序列的 DFT。同样，这些序列可以继续分解下去，继续得到更短的序列，从而可以更简便地进行运算。这种算法同样是主要针对2的整数次幂长度的序列的。

从本章前面对DFT的介绍和本节开头的分析可知，随着序列长度的减小，FFT运算的复杂度将以指数规律降低。

本节主要讨论序列长度是2的整数次幂时的DFT运算，这称为基−2FFT。除了基−2FFT，还有基4−FFT和基−8FFT，甚至还有基−6FFT。那些算法的效率比基−2FFT更高，但应用的范围更狭窄。事实上，很多商业化的信号分析库都是使用混合基FFT的。那样的程序代码更加复杂，但效率却高得多，而且应用范围更广。本书从学习和研究的角度，仅介绍最常见的按时间抽取的基−2FFT算法。

对于基−2的FFT，可以设序列长度为N=2^L 。由于N是偶数，我们可将这个序列按照项数的奇偶分成2组。分组的规律如下式所示。

则f(x)的傅立叶变换F(u)可以表示为f(x)的奇数项和偶数项分别组成的序列的如下变换形式。

因为

，所以上式可以继续化简为

容易发现，上式的第一项为f(2r)的N/2点DFT，而第二项的求和部分为f(2r+1)的N/2点DFT （序列f(2r)和序列f(2r+1)的周期均为N/2）。也即

这里，我们用f_偶(u)和f_奇(u)分别表示f(2r)和f(2r+1)的N/2点DFT。

而且，根据DFT序列的周期性特点，还可得到如下式子成立。

并且，由于

，我们还可以得出

将式（6-36）和式（6-37）代入式（6-38），并根据式（6-35），得

这是一个递推公式，它就是 FFT 蝶形运算的理论依据。该公式表明，一个偶数长度序列的傅立叶变换可以通过它的奇数项和偶数项的傅立叶变换得到，从而可以将输入序列分成两部分分别计算，并按公式相加/相减。而在这个运算过程中，实际上只需要计算

，u=0,1,2,3,…N/2−1。

因此，一个8点按时间抽取的FFT算法的第一步骤如图6.7所示。

图6.7所示是根据式（6-39）绘制的，这一算法也可以用图6.8抽象地表示出来。

由于我们讨论的是基-2 的 FFT 算法，N/2 一般应是偶数，因此得到的序列还可以继续分解，分解过程可以一直持续到每个序列只需要2点的DFT。这样只需要如下的运算即可计算这一DFT值，这一运算是FFT的基本运算，称为蝶形运算。

这一基础单元是对初始输入序列进行傅立叶变换操作的第一步，即2点时的FFT。把这个基本的DFT运算和上面的抽象化蝶形运算比较，可以发现它们的基本结构是完全一致的。在蝶形算法中，我们可以只计算一次

，而后分别与f₁相加和相减，从而每一次蝶形算法只需1次复数乘法和2次复数加法（从复杂度分析的角度，相减当然也可看作是一次加法）。并且，注意到

，因此可以进一步简化计算。尤其第一级蝶形运算更是可以完全简化为单纯的复数加减法。

一个8点FFT的完整计算过程如图6.10所示，请思考这个过程与DFT过程的区别，以及这个过程所需的算法复杂度和存储空间问题。稍后我们将讨论这个问题。

用基-2的时间抽取FFT算法比直接计算DFT的效率高得多。在计算长度为N=2^L 序列的FFT时，在不对复数乘法进行额外优化的情况下，所需运算量分析如下。

对于每一个蝶形运算，我们需要进行1次复数乘法和2次复数加法。而FFT运算的每一级都含有N/2 = 2^L-1 个蝶形运算单元。因此，完成L级FFT运算共需要的复数乘法次数m_cm和复数加法m_ca数目分别为

而本节开头提到，实现同样长度序列的DFT运算则需要n² 次复数乘法和N(N-1)次复数加法，这远远多出FFT算法的所需。近似比较FFT和DFT运算的算法复杂度可知

因此，在N或L取值增大时，FFT运算的优势更加明显。例如，当N=2¹⁰时，C(N)=102.4，即FFT算法的速度是DFT的102.4倍。

此外，从占用的存储空间看，按时间抽取的FFT算法也远比DFT算法节约。一对复数进行完蝶形运算后，就没有必要再次保留输入的复数对。因此，输出对可以和输入对放在相同的存储单元中。所以，只需要和输入序列大小相等的存储单元即可。也就是一种“原位运算”。

但是，经过观察上面的8点FFT运算全过程，我们发现，如果要使用这种“原位”运算，输入序列就必须按照倒序存储。由于f(x)是逐次抽取的，所以必须对原输入码列倒转位序，得到的次序相当于是原序列编号的二进制码位倒置。也即将原序列编号按照二进制表示，并且将二进制的所有位次序颠倒，就得到了在实际的输入序列中应该使用的排序位置。下面同样以8点FFT为例说明码位倒置的方法，如表6.1所示。

按照表6.1中的顺序排列输入数据，就可以方便地进行原位运算，以节约内存空间。

离散反傅立叶变换的形式与离散傅立叶变换很相似,首先比较它们的公式形式。

离散反傅立叶变换IDFT的公式为

离散傅立叶变换DFT的公式为

观察式（6-44）和式（6-45）发现，只要把DFT算子中的

换成

，并在前面乘以1/N，即可得到IDFT的算子。于是考虑使用复数共轭方式建立两者之间的联系，推导离散反傅立叶变换的公式如下。

因此，我们只需先将F(u)取共轭，就可以直接使用FFT算法计算IFFT了。

N维快速傅立叶变换用于对高维信号矩阵执行傅立叶频谱分析操作。其中二维的快速傅立叶变换常常用于数字图像处理。

N 维快速傅立叶变换是由一维 FFT 组合而成的，其运算实质就是在给定二维或多维数组的每个维度上依次执行一维FFT，并且使用“原位运算”的方法。在开始之前，算法将输入直接复制到输出上，所以之后在每个维度上执行 FFT 的原位操作都不会改变原本的输入数组，同时也使这个算法输出的数组和输入的数组拥有同样的大小和维度。也就是说，如果对一幅灰度图像执行二维快速傅立叶变换操作，得到的结果也将是一个二维数组。

Matlab中提供了fft2和ifft2函数分别计算二维傅立叶变换和反变换，它们都经过了优化，运算速度非常快；另一个与傅立叶变换密切相关的函数是fftshift，常需要利用它来将傅立叶频谱图中的零频点移动到频谱图的中心位置。

下面分别介绍这3个函数。

该函数用于执行二维快速傅立叶操作，因此可以直接用于数字图像处理。调用语法为

● X为输入图像。

● m和n分别用于将X的第一和第二维规整到指定的长度。当m和n均为2的整数次幂时，算法的执行速度要比m和n均为素数时更快。

● Y是计算得到的傅立叶频谱，是一个复数矩阵。

计算abs(Y)可得到幅度谱，计算angle(Y)可得到相位谱。

在fft2函数输出的频谱分析数据中，是按照原始计算所得的顺序来排列频谱的，而没有以零频为中心来排列，因此造成了零频在输频谱矩阵的角上，显示幅度谱图像时表现为4个亮度较高的角（零频处的幅值较高），如图6.11（a）所示。

fftshift函数利用了频谱的周期性特点，将输出图像的一半平移到另一端，从而使零频被移动到图像的中间。其调用语法为

● X为要平移的频谱。

● dim指出了在多维数组的哪个维度上执行平移操作。

● Y是经过平移的频谱。

利用fftshift函数对图6.11（a）中的图像平移后的效果如图6.11（b）所示。

下面给出对于二维图像矩阵，fftshift函数的平移过程，如图6.12所示。

可见，输出矩阵被分为了4个部分，其中1、3两部分对换，2、4两部分对换。这样，原来在角上的零频率点（原点）位置就移动到了图像的中央位置。而dim参数则可以指定在多维数组的哪个维度上执行对换操作。例如，对于矩阵而言，dim取1和2的情形分别如图6.13所示。

该函数用于对图像（矩阵）执行逆傅立叶变换。输出矩阵的大小与输入矩阵相同。调用形式为

● X为要计算反变换的频谱。

● m、n的意义与fft2中相同。

● Y是反变换后得到的原始图像。

在执行IFFT2函数之前，如果曾经使用FFTSHIFT函数对频域图像进行过原点平移，则还需要使用IFFTSHIFT将原点平移回原位置。

[例6.1] 幅度谱的意义

下面的程序展示了如何利用fft2进行二维快速傅立叶变换。为了更好地显示频谱图像，需要利用3.3节中学习过的对数变换来增强频谱。

I1 = imread('cell.tif');　%读入原图像

fcoef = fft2(double(I1));　%做fft变换

spectrum = fftshift(fcoef);　%将零点移到中心

temp =log(1+abs(spectrum));　%对幅值做对数变换以压缩动态范围

I2 = imread('circuit.tif');　%读入原图像

fcoef = fft2(double(I2));　%做fft变换

spectrum = fftshift(fcoef);　%将零点移到中心

temp =log(1+abs(spectrum));　%对幅值做对数变换以压缩动态范围

上述程序的运行结果如下。

由图6.14所示可以看出，图6.14（b）中的cell.tif图像较为平滑，而在其傅立叶频谱中，低频部分对应的幅值较大；而对图6.14（d）中细节复杂的的图像circuit.tif，灰度的变化趋势更加剧烈，相应的频谱中高频分量较强。

事实上，由于图6.14（b）图中基本只存在水平和垂直的线条，导致了在输出的频谱中亮线集中存在于水平和垂直方向（并且经过原点）。具体地说，原图像中的水平边缘对应于频谱中的竖直亮线，而竖直边缘则对应着频谱中的水平响应。我们不妨这样理解，水平方向的边缘可以看作在竖直方向上的灰度值的矩形脉冲，而这样的矩形脉冲可以分解为无数个竖直方向正弦平面波的叠加，从而对应频域图像中的垂直亮线；而对于竖直方向的边缘，情况是类似的。

通过例6.1可以发现一些频谱与其空间域图像之间的联系。实际上，低频（频谱图像中靠近中心的区域）对应着图像的慢变化分量；高频（频谱图像中远离中心的区域）对应着一幅图像中较快变化的灰度级，常常对应着图像细节，如物体的边缘和噪声等。就拿图6.14（c）的电路图像来说，电路板的灰度较为一致的背景区域就对应着频谱的低频部分，而横竖电路线条的灰度变换则是相对高频的成份，且灰度变换越剧烈，就对应着越高的频域分量。

我们在6.2.3小节曾给出了幅度谱和相位谱的定义，并对其作用进行了简单的介绍。为了进一步加深读者对幅度谱和相位谱的认识，这里给出一个关于它们的有趣的例子。

[例6.2]美女与猫——交换两幅图像的相位谱

图6.15（a）、（b）中分别是一张美女的照片和一张猫的照片，这里我们准备交换这两幅图像的相位谱，即用美女的幅度谱加上猫的相位谱，而用猫的幅度谱加上美女的相位谱，然后根据式（6-18），通过幅度谱和相位谱来还原傅立叶变换F(u,v)，再经傅立叶反变换得到交叉相位谱之后的图像。根据6.2.2中关于幅度谱和相位谱各自作用的讨论，您能想到这样做将会产生怎样的结果吗?

A = imread('../beauty.jpg');

B = imread('../cat.jpg');

% 求傅立叶变换

Af = fft2(double(A));

Bf = fft2(double(B));

% 分别求幅度谱和相位谱

% 交换相位谱并重建复数矩阵

AfR = AfA .* cos(BfB) + AfA .* sin(BfB) .* i;

BfR = BfA .* cos(AfB) + BfA .* sin(AfB) .* i;

% 傅立叶反变换

AR = abs(ifft2(AfR));

BR = abs(ifft2(BfR));

title('美女原图像');

title('猫的原图像');

title('美女的幅度谱和猫的相位谱组合');

title('猫的幅度谱和美女的相位谱组合');

程序运行结果如图6.15所示。

通过这个示例，我们可以发现，交换相位谱之后，反变换之后得到的图像内容与其相位谱对应的图像一致，这就印证了我们之前关于相位谱决定图像结构的论断。而图像中整体灰度分布的特性，如明暗、灰度变化趋势等则在比较大的程度上取决于对应的幅度谱，因为幅度谱反映了图像整体上各个方向的频率分量的相对强度。

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

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

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

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

根据式（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小节详细地探讨一些具有更高实用性的频域滤波器。

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

通常使用fftshift函数将频谱原点移至图像中心，因此需要构造对应的原点在中心的滤波器，并在滤波之后使用ifftshift函数将原点移回以进行反变换。

频域滤波算法 imfreqfilt 的完整实现如下，该算法被封装在金羽图书论坛（ http://bbs.book 95.com）的“金羽图书与答疑”板块与本书同名的主题帖子附件的“chapter6/code”目录下的imfreqfilt.m文件中。

function out = imfreqfilt(I, ff)

% imfreqfilt函数　对灰度图像进行频域滤波

% 参数I　输入的空域图像

% 参数ff　应用的与原图像等大的频域滤镜

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

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);

% 快速傅立叶变换

f = fft2(double(I));

% 应用滤镜及反变换

out = s .* ff; %对应元素相乘实现频域滤波

out = ifftshift(out);

% 归一化以便显示

out = out/max(out(:));

在频谱中，低频主要对应图像在平滑区域的总体灰度级分布，而高频对应图像的细节部分，如边缘和噪声。因此，图像平滑可以通过衰减图像频谱中的高频部分来实现，这就建立了空间域图像平滑与频域低通滤波之间对应关系。

最容易想到的衰减高频成份的方法就是在一个称为“截止频率”的位置“截断”所有的高频成份，将图像频谱中所有高于这一截止频率的频谱成份设为0，低于截止频率的成份设为保持不变。能够达到这种效果的滤波器如图6.16所示，我们称之为理想低通滤波器。如果图像的宽度为M，高度为N，那么理想低通频域滤波器可形式化地描述为

其中D₀表示理想低通滤波器的截止频率，滤波器的频率域原点在频谱图像的中心处，在以截止频率为半径的圆形区域之内的滤镜元素值全部为1，而该圆之外的滤镜元素值全部为0。理想低通滤波器的频率特性在截止频率处十分陡峭，无法用硬件实现，这也是我们称之为理想的原因，但其软件编程的模拟实现较为简单。

理想低通滤波器可在一定程度上去除图像噪声，但由此带来的图像边缘和细节的模糊效应也较为明显，其滤波之后的处理效果比较类似于5.3.1中的平均平滑。实际上，理想低通滤波器是一个与频谱图像同样尺寸的二维矩阵，通过将矩阵中对应较高频率的部分设为0，较低频率的部分（靠近中心）设为1，可在与频谱图像相乘后有效去除频谱的高频成份（由于是矩阵对应元素相乘，频谱高频成份与滤波器中的0相乘）。其中0与1的交界处即对应滤波器的截止频率。

利用我们编写的imidealflpf函数可以得到截止频率为freq的理想低通滤波器。

imidealflpf 函数的完整实现如下，该函数被封装在金羽图书论坛（http://bbs.book95.com）的“金羽图书与答疑”板块与本书同名的主题帖子附件的“chapter6/code”目录下的imidealflpf.m文件中。

function out = imidealflpf(I, freq)

% imidealflpf函数　构造理想的频域低通滤波器

% I参数　输入的灰度图像

% freq参数　低通滤波器的截止频率

% 返回值：out – 指定的理想低通滤波器

if (sqrt(((i-M/2)^2+(j-N/2)^2))>freq)

下面我们仍以第 5 章中图像平滑中曾使用过的婴儿老照片 baby_noise.bmp 为例，使用频域理想低通滤波器进行处理，相应的Matlab代码如例6.3所示。

[例6.3]理想低通滤波实现。

I = imread('../baby_noise.bmp'); %读入原图像

ff = imidealflpf(I, 20);

out = imfreqfilt(I, ff);

% 计算FFT并显示

temp = fft2(double(I));

temp = fftshift(temp);

temp = log(1 + abs(temp));

title('Ideal LPF, freq=20');

% 计算FFT并显示

temp = fftshift(temp);

temp = log(1 + abs(temp));

title(' Ideal LPF, freq=20');

ff = imidealflpf(I, 40);

out = imfreqfilt(I, ff);

title('Ideal LPF, freq=40');

% 计算FFT并显示

temp = fftshift(temp);

temp = log(1 + abs(temp));

title(' Ideal LPF, freq=40');

ff = imidealflpf(I, 60);

out = imfreqfilt(I, ff);

title('Ideal LPF, freq=60');

% 计算FFT并显示

temp = fftshift(temp);

temp = log(1 + abs(temp));

title(' Ideal LPF, freq=60');

上述程序的运行效果如图6.17和图6.18所示。

从图6.17和图6.18可见，当截止频率非常低时，只有非常靠近原点的低频成份能够通过，图像模糊严重；截止频率越高，通过的频率成份就越多，图像模糊的程度越小，所获得的图像也就越接近原图像。但可以看出，理想低通滤波器并不能很好地兼顾噪声滤除与细节保留两个方面，这和空域中采用平均模板时的情形比较类似。下面，我们将介绍频域的高斯低通滤波器，并比较它与理想低通滤波器的处理效果。

高斯低通滤波器的频率域二维形式由下式给出。

高斯函数具有相对简单的形式，而且它的傅立叶变换和傅立叶反变换都是实高斯函数。为了简单起见，下面仅给出一个一维高斯函数的傅立叶变换和傅立叶反变换作为例子，式（6-52）告诉我们一维高斯函数的傅立叶反变换和正变换仍为高斯函数，该式的证明留给有兴趣的读者自己完成（提示：可以利用高斯分布的概率密度函数在定义域上积分为1的性质）。

式中：σ是高斯曲线的标准差。

频域和与之对应的空域一维高斯函数的图形如图6.19所示。

从图6.19可以发现，当σ增大时，H(u)的图像倾向于变宽，而h(x)的图像倾向于变窄和变高。这也体现了频率域和空间域的对应关系。频率域滤波器越窄，滤除的高频成份越多，图像就越平滑（模糊）；而在空间域，对应的滤波器就越宽，相应的卷积模板越平坦，平滑（模糊）效果就越明显。

我们在图6.20中分别给出σ取值为3和1时的频域二维高斯滤波器的三维曲面表示，可以看出频域下的二维高斯滤波器同样具有上述一维情况时的特点。

根据上面二维高斯低通滤波器的定义，可以编写高斯低通滤波器的生成函数如下，该函数被封装在金羽图书论坛（http://bbs.book95.com）的“金羽图书与答疑”板块与本书同名的主题帖子附件的“chapter6/code”目录下的imgaussflpf.m文件中。

function out = imgaussflpf(I, sigma)

% imgaussflpf函数　构造频域高斯低通滤波器

% I参数　输入的灰度图像

% sigma参数　高斯函数的Sigma参数

out(i,j) = exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);

下面的例6.4给出了针对图像baby_noise.bmp，sigma不同取值时高斯低通滤波的Matlab程序。

[例6.4]高斯低通滤波实现。

I = imread('../baby_noise.bmp'); %读入原图像

ff = imgaussflpf (I, 20);

out = imfreqfilt(I, ff);

% 计算FFT并显示

temp = fft2(double(I));

temp = fftshift(temp);

temp = log(1 + abs(temp));

title('Gauss LPF, sigma=20');

% 计算FFT并显示

temp = fftshift(temp);

temp = log(1 + abs(temp));

title(' Gauss LPF, sigma=20');

ff = imgaussflpf (I, 40);

out = imfreqfilt(I, ff);

title('Gauss LPF, sigma =40');

% 计算FFT并显示

temp = fftshift(temp);

temp = log(1 + abs(temp));

title(' Gauss LPF, sigma =40');

ff = imgaussflpf (I, 60);

out = imfreqfilt(I, ff);

title('Gauss LPF, sigma =60');

% 计算FFT并显示

temp = fftshift(temp);

temp = log(1 + abs(temp));

title(' Gauss LPF, sigma =60');

上述程序的运行后得到的滤波效果如图6.21所示。

图6.21中各幅图像所对应的频域图像如图6.22所示。显然，高斯滤波器的截止频率处不是陡峭的。

高斯低通滤波器在Sigma参数取40的时候可以较好地处理被高斯噪声污染的图像，而相比于理想低通滤波器而言，处理效果上的改进是显而易见的。高斯低通滤波器在有效抑制噪声的同时，图像的模糊程度更低，对边缘带来的混叠程度更小，从而使高斯低通滤波器在通常情况下获得了比理想低通滤波器更为广泛的应用。

图像锐化可以通过衰减图像频谱中的低频成份来实现，这就建立了空间域图像锐化与频域高通滤波之间对应关系。

temp = log(1 + abs(temp));

从6.5.2小节高斯低通滤波器中的H(u)图像，可以发现滤波器中心频率处的值即为其最大值1，如果需要做相反的滤波操作，滤除低频成份而留下高频成份，则可以考虑简单地使用如下表达式来获得一个高斯高通滤波器。

因此，高斯高通滤波器的频域特性曲线如图6.23所示（仍旧以一维情况为例）。

根据上面二维高斯高通滤波器的定义，可以编写高斯高通滤波器的生成函数如下，该函数被封装在金羽图书论坛（http://bbs.book95.com）的“金羽图书与答疑”板块与本书同名的主题帖子附件的“chapter6/code”目录下的imgaussfhpf.m文件中。

function out = imgaussfhpf(I, sigma)

% imgaussfhpf函数　构造频域高斯高通滤波器

% I参数　输入的灰度图像

% sigma参数　高斯函数的Sigma参数

out(i,j) = 1 - exp(-((i-M/2)^2+(j-N/2)^2)/2/sigma^2);

下面的例6.5给出了针对Matlab示例图像coins.png，sigma取不同值时高斯高通滤波的Matlab程序。

[例6.5]高斯高通滤波实现。

I = imread('coins.png');

ff = imgaussfhpf (I, 20);

out = imfreqfilt(I, ff);

% 计算FFT并显示

temp = fft2(double(I));

temp = fftshift(temp);

temp = log(1 + abs(temp));

title('Gauss HPF, sigma=20');

% 计算FFT并显示

temp = fftshift(temp);

temp = log(1 + abs(temp));

title('Gauss HPF, sigma=20');

ff = imgaussfhpf (I, 40);

out = imfreqfilt(I, ff);

title('Gauss HPF, sigma =40');

% 计算FFT并显示

temp = fftshift(temp);

title('Gauss HPF, sigma =40');

ff = imgaussfhpf(I, 60);

out = imfreqfilt(I, ff);

title('Gauss HPF, sigma =60');

% 计算FFT并显示

temp = fftshift(temp);

temp = log(1 + abs(temp));

title('Gauss HPF, sigma =60');

上述程序运行后得到的滤波效果如图6.24所示。

滤波前后对应的频域图像如图6.25所示。

高斯高通滤波器可以较好地提取图像中的边缘信息，Sigma参数取值越小，边缘提取越不精确，会包含越多的非边缘信息；Sigma参数取值越大，边缘提取越精确，但可能包含不完整的边缘信息。

频域拉普拉斯算子的推导可以从一维开始，由傅立叶变换的性质可知

因此拉普拉斯算子的傅立叶变换计算如下。

因此有下式成立。

也即频域的拉普拉斯滤波器为

根据频域图像频率原点的平移规律，将上式改写为

根据式（6-58），可以编写拉普拉斯频域滤波器的生成函数如下，该函数被封装在金羽图书论坛（http://bbs.book95.com）的“金羽图书与答疑”板块与本书同名的主题帖子附件的“chapter6/code”目录下的imlapf.m文件中。

function out = imlapf(I)

% imlapf函数　构造频域拉普拉斯滤波器

% I参数　输入的灰度图像

out(i,j) = -((i-M/2)^2+(j-N/2)^2);

下面的例6.6给出了对Matlab示例图像coins.png进行频域拉普拉斯滤波的Matlab程序。

[例6.6] 拉普拉斯滤波实现。

I = imread('coins.png');

out = imfreqfilt(I, ff);

temp = fft2(double(I));

temp = fftshift(temp);

temp = log(1 + abs(temp));

title('Laplace Filter');

temp = fftshift(temp);

temp = log(1 + abs(temp));

title('Laplace Filter');

上述程序运行后得到的滤波效果如图6.26所示。

得到的滤波前后的频域图像如图6.27所示。

6.5～6.6 节介绍了几种典型的频域滤波器，实现了频域下的低通和高通滤波，它们均可在空域下采用平滑和锐化算子实现。而本节准备给出一个特别适合在频域中完成的滤波案例，即利用频域带阻滤波器消除图像中的周期噪声。下面就来看一下这个在空域中几乎不可能完成的任务是如何在频域中实现的。

顾名思义，所谓“带阻”就是阻止频谱中某一频带范围的分量通过，其他频率成份则不受影响。常见的带阻滤波器有理想带阻滤波器和高斯带阻滤波器。

理想带阻滤波器的表达式为

式中：D₀是阻塞频带中心频率到频率原点的距离；W是阻塞频带宽度；D是（u,v）点到频率原点的距离。于是，理想带阻滤波器的频域特性曲面如图6.28所示。

本案例中使用了高斯带阻滤波器，下面直接给出高斯带阻滤波器的表达式。

式中：D₀是阻塞频带中心频率到频率原点的距离；W是阻塞频带宽度；D是（u,v）点到频率原点的距离。于是，二维高斯带阻滤波器的频域特性曲面如图6.29所示。

根据高斯带阻滤波器的定义式（6-60），可以编写高斯带阻滤波器的生成函数如下，该函数被封装在金羽图书论坛（http://bbs.book95.com）的“金羽图书与答疑”板块与本书同名的主题帖子附件的“chapter6/code”目录下的imgaussfbrf.m文件中。

function out = imgaussfbrf(I, freq, width)

% imgaussfbrf函数　构造频域高斯带阻滤波器

% I参数　输入的灰度图像

% freq参数　阻带中心频率

% width参数　阻带宽度

out(i,j) = 1-exp(-0.5*((((i-M/2)^2+(j-N/2)^2)-freq^2)/(sqrt(i.^2+j.^2)*width))^2);

带阻滤波器常用于处理含有周期性噪声的图像。周期性噪声可能由多种因素引入，例如图像获取系统中的电子元件等。本案例中我们人为地生成了一幅带有周期噪声的图像，而后通过观察分析其频谱特征，选择了合适的高斯带阻滤波器进行频域滤波。

通常可以使用正弦平面波来描绘周期性噪声。如下程序为Matlab示例图片pout.tif增加周期性噪声。

O = imread('pout.tif'); %读入原图像

I(i,j)=I(i,j)+20*sin(20*i)+20*sin(20*j); %添加周期噪声

title('Added Noise');

添加周期性噪声前后的区别如图6.30所示。

使用高斯带阻滤波器时，首先需要对欲处理的图像的频谱有一个了解。下面的命令得到了两幅图像的频谱。

imshow(o_f, [ ]); %得到图6.31（a）

imshow(i_f, [ ]); %得到图6.31（b）

title('Added Noise');

程序的运行结果如图6.31所示。

观察图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）

title('Gauss Filter');

上述程序运行后得到的高斯带阻滤波器最终滤波效果如图 6.32 所示，我们看到周期噪声被很好地消除，这样的效果在空域中是很难实现的。

在6.4.1小节我们曾探讨了频域滤波与空域滤波之间的关系。这里则要更进一步，来研究频域滤波器与空域滤波器之间的内在联系。

频域滤波较空域而言更为直观，频域下滤波器表达了一系列空域处理（平滑、锐化等）的本质，即对高于/低于某一特定频率的灰度变化信息予以滤除，而对其他的灰度变化信息基本保持不变。这种直观性增加了频域滤波器设计的合理性，使得我们更容易设计出针对特定问题的频域滤波器，就如在6.7节中我们利用了带阻滤波器实现了对图像中周期噪声的滤除，而想直接在空域中设计出一个能够完成如此滤波任务的滤波器（卷积模板）是相当困难的。

为了得到合适的空域滤波器，我们很自然地想到可以首先设计频域滤波器 H(u, v)，而后根据卷积定理式（6-61)，将H(u, v)反变换至空域后就得到了空域中滤波使用的卷积模板h(x, y)，从而解决了空域滤波器的设计难题。

然而，直接反变换得到的空域卷积模板h(x, y)同 H(u, v)等大，从而与图像f(x, y)具有相同的尺寸。而模板操作十分耗时，要计算这样大的模板与图像的卷积将是非常低效的。在第3章中我们使用的都是很小的模板（如 3×3, 5×5, 7×7等），因为这样的模板在空域中才具有滤波效率上的优势。一般来说，如果空域模板中的非零元素数目小于132（大约13×13见方），则直接在空域中计算卷积较为划算，否则直接利用H(u, v)在频域下滤波更为合适。

在实际中我们发现，利用以全尺寸的空域滤波器h(x, y)为指导设计出的形状与之类似的小空域卷积模板，同样可以取得类似于频域滤波器 H(u, v)的滤波效果。这就为从频域出发，最终设计出具有实用价值的空域模板提供了一种完美的解决方案。

式（6-52）给出的高斯频域低通滤波器 H(u)及与其构成傅立叶变换对儿的空域高斯模板 h(x)正好印证了上述结论。从图6.19上来看，H(u)越窄，h(x)就越宽。而频域低通滤波器H(u)越窄，说明能够通过的频率越低，被截断的高频成份也就越多，从而使滤波处理后原函数f(x)变得平滑；而空域下以越宽的模板h(x)与函数f(x)卷积则同样会产生平滑的效果。再进一步以h(x)的形状为指导，就可以得到曾在高斯平滑中使用的高斯模板式（5-5）。

傅立叶级数的收敛性：满足狄利赫里条件的周期函数表示成的傅里叶级数都收敛。狄利赫里条件如下。

（1）在任何周期内，x(t)须绝对可积。

（2）在任一有限区间中，x(t)只能取有限个最大值或最小值。

（3）在任何有限区间上，x(t)只能有有限个第一类间断点。

式（6-1）可转化为

则式(6-62)可表示为

可将上式合并，得到傅立叶级数的复数形式。

为求得系数c_n,将式（6-2）代入式（6-63），得

同样的，可将结果合并写为

即傅立叶系数的复数形式。

式中: u表示f(x)的频率，即

由复指数形式的傅立叶级数

∵c_n=|c_n|e^iθn ，其中

即c_n的虚部比c_n的虚部，为n次频率的谐波的相位。

又∵|c_-n|=|c_n|，而θ_-n=−θ_n⇒c_n=|c_n|e^-iθn

根据欧拉公式展开2个复指数项易得傅立叶级数展开的复指数形式。

精通Matlab数字图像处理与识别/张铮等编著.--北京：人民邮电出版社，2013.4

ISBN　978-7-115-30463-6

Ⅰ.①精…　Ⅱ.①张…　Ⅲ.①Matlab软件—应用—数字图象处理　Ⅳ.①TN911.73

中国版本图书馆CIP数据核字（2012）第307788号

本书将理论知识、科学研究和工程实践有机结合起来，内容涉及数字图像处理和识别技术的方方面面，包括图像的点运算、几何变换、空域和频域滤波、小波变换、图像复原、形态学处理、图像分割以及图像特征提取的相关内容；同时，对于机器视觉进行了前导性的探究，重点介绍了两种目前在工程技术领域非常流行的分类技术——人工神经网络（ANN）和支持向量机（SVM），并讲解了人脸识别的热点技术。

全书结构紧凑，内容深入浅出，讲解图文并茂，适合计算机、通信和自动化等相关专业的本科生、研究生，以及工作在图像处理和识别领域一线的广大工程技术人员参考使用。

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

图书目录:

详情

图书摘要

相关图书

相关文章

相关课程