第5章频率域图像处理

频率域图像是把空间域图像像素的灰度值表示成随位置变化的空间频率，并以频谱(也称为频谱图))的形式表示图像信息分布特征的一种表示方式。频率域图像处理是指在图像的频率域中对图像进行某种处理的方法，这种方法以傅里叶变换为基础，即先通过傅里叶变换把图像从空间域变换到频率域，然后用频率域方法对图像进行处理，处理完后再利用傅里叶反变换把图像变换回空间域。
本章首先介绍面向空间域图像的二维离散傅里叶变换及其频谱特性，然后介绍频率域图像处理的基本实现思路、基于频率域的图像噪声消除(频率域低通滤波)和基于频率域的图像增强(频率域高通滤波)方法，最后介绍面向某些特殊应用的带阻滤波和带通滤波。
5．1二维离散傅里叶变换
1822年，法国工程师傅里叶(Fourier)指出： 一个“任意”的周期函数f(x)都可以分解为无穷多个不同频率的正弦和/或余弦的和，即傅里叶级数。求解傅里叶级数的过程就是傅里叶变换。傅里叶级数和傅里叶变换又统称为傅里叶分析或谐波分析。
离散傅里叶变换(discrete Fourier transform，DFT)描述的是离散信号的一维时域或二维空间域表示与频域表示的关系，是信号处理和图像处理中的一种最有效的数学工具之一，在频谱分析、数字滤波器设计、功率谱分析、传递函数建模、图像处理等方面具有广泛的应用。
5．1．1二维离散傅里叶变换的定义和傅里叶频谱
由于一幅数字图像可以描述成一个二维函数，所以下面仅介绍应用于图像处理的二维离散傅里叶变换。
1． 二维离散傅里叶变换的定义
设f(x,y)是在空间域上等间隔采样得到的M×N的二维离散信号，x和y是离散实变量，u和v为离散频率变量，则二维离散傅里叶变换对一般定义为

F(u,v)=1MN∑M－1x=0∑N－1y=0f(x,y)exp－j2πxuM+yvN
(u=0,1,…,M-1； v=0,1,…,N-1)(5．1)
f(x,y)=1MN∑M－1u=0∑N－1v=0F(u,v)expj2πuxM+vyN
(x=0,1,…,M-1； y=0,1,…,N-1)(5．2)

其中，式(5．1)为傅里叶正变换，式(5．2)为傅里叶反(逆)变换。
在图像处理中，有时为了讨论上的方便，取M=N，这样二维离散傅里叶变换对就定义为

F(u,v)=1N∑N－1x=0∑N－1y=0f(x,y)exp－j2π(xu+yv)N
(u,v=0,1,…,N-1)(5．3)
f(x,y)=1N∑N－1u=0∑N－1v=0F(u,v)expj2π(ux+vy)N
(x,y=0,1,…,N-1)(5．4)

其中，exp［－j2π(xu+yv)/N］是正变换核，exp［j2π(ux+vy)/N］是反变换核。




二维离散傅里叶变换的频谱和相位角定义为

|F(u,v)|=R2(u,v)+I2(u,v)(5．5)
(u,v)=arctan［I(u,v)/R(u,v)］(5．6)

2． 图像的傅里叶频谱特性及频谱图
傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，从物理效果看，傅里叶变换是将图像从空间域转换到频率域。一幅空间域的图像f(x,y)变换到频率域F(u,v)后，其傅里叶频谱|F(u,v)|可以以频谱图的形式予以显示。比如，图5．1(a)的图像的频谱如图5．1(b)所示； 图5．2(a)的图像的频谱如图5．2(b)所示。


图5．1图像的傅里叶频谱示例1




图5．2图像的傅里叶频谱示例2


在图像的傅里叶频谱中，原空间域图像上的灰度突变部位、图像结构复杂的区域、图像细节及干扰噪声等信息集中在高频区，对应于图5．1(b)和图5．2(b)的傅里叶频谱的中间部位； 原空间域图像上灰度变化平缓部位(即空间域图像的平坦区域)的信息集中在低频区，对应于图5．1(b)和图5．2(b)的傅里叶频谱的4个角部分。即在傅里叶频谱图上，比较亮的4个角反映的是原图像的低频特性，亮度越大说明低频能量越大； 比较暗的中间部位反映的是原图像的高频特性。
按照图像空间域和频率域的对应关系，空间域中的强相关性，即图像中一般都存在有大量的平坦区域，使得图像平坦区域中的相邻或相近像素一般趋向于取相同或相近的灰度值，它们在频率域中对应于低频部分，即对应于傅里叶频谱的4个角部分。由于低频部分能量较集中，因而在频谱图上的视觉效果较亮。
5．1．2二维离散傅里叶变换的重要性质
二维离散傅里叶变换的性质包括线性性、可分离性、平均值性质、周期性、共轭对称性、空间位置和空间频率的平移性、旋转性、尺度变换性、卷积性质等，下面仅介绍几种比较重要且在本书的有关内容中涉及的性质。
1．  变换系数矩阵
由二维离散傅里叶反变换式(5．4)可知，由于u和v均有0，1，…，N-1的N个可能的取值，所以f(x,y)由N2个频率分量组成，所以每个频率分量都与一个特定的(u,v)值相对应； 且对于某个特定的(u,v)值来说，当(x,y)取遍所有可能的值(x=0，1，…，N-1； y=0，1，…，N-1)时，就可得到对应于该特定的(u,v)值的一个变换系数矩阵： 

expj2π0u+0vN
expj2π0u+1vN
…
expj2π0u+(N-1)vN

expj2π1u+0vN
expj2π1u+1vN
…
expj2π1u+(N-1)vN

︙︙︙

expj2π(N-1)u+0vN
expj2π(N-1)u+1vN
…
expj2π(N-1)u+(N-1)vN
(5．7)


显然，对应于不同(u,v)值的变换系数矩阵共有N2个，且它们与f(x,y)无关。
2． 可分离性
式(5．3)和式(5．4)的二维离散傅里叶变换对可写成如下的分离形式： 

F(u,v)=1N∑N－1x=0∑N－1y=0f(x,y)exp－j2πyvNexp－j2πxuN
(u,v=0,1,…,N-1)(5．8)
f(x,y)=1N∑N－1u=0∑N－1v=0F(u,v)expj2πvyNexpj2πuxN
(x,y=0,1,…,N-1)(5．9)


上述的分离表示形式说明，可以连续运用两次一维离散傅里叶变换来实现一个二维离散傅里叶变换。以式(5．8)为例，可先沿y轴方向进行一维的(列)变换而求得

F(x,v)=1N∑N－1y=0f(x,y)exp－j2πvyN
(v=0,1,…,N-1)(5．10)

然后再对F(x,v)沿x方向进行一维的(行)变换而得到最后结果

F(u,v)=1N∑N－1x=0F(x,v)exp－j2πuxN
(u,v=0,1,…,N-1)(5．11)

3． 平均值
一幅图像的灰度平均值可表示为

f-=1N2∑N－1x=0∑N－1y=0f(x,y)(5．12)

如果将u=v=0代入式(5．3)，可得

F(0,0)=1N∑N－1x=0∑N－1y=0f(x,y)(5．13)

所以，一幅图像的灰度平均值可由离散傅里叶变换在原点处的值求得，即

f-=1NF(0,0)(5．14)


对于M×N的图像f(x,y)和二维离散傅里叶变换对的一般定义式(5．1)和式(5．2)，图像的灰度平均值公式为

f-=1MNF(0,0)(5．15)

4． 周期性
对于M×N的图像f(x,y)和二维离散傅里叶变换对的一般定义式(5．1)和式(5．2)，F(u,v)的周期性定义为

F(u,v)=F(u+mM,v+nN)(m,n=0,±1，±2,…)(5．16)

5． 共轭对称性
设f(x,y)为实函数，则其傅里叶变换F(u,v)具有共轭对称性，表示为

F(u,v)=F(－u,－v)(5．17)
|F(u,v)|= |F(－u,－v)|(5．18)

6． 平移性
对于M×N的图像f(x,y)和二维离散傅里叶变换对的一般定义式(5．1)和式(5．2)，设用符号表示函数与其傅里叶变换的对应性，则傅里叶变换的平移性可表示为

f(x,y)expj2πu0xM+v0yNF(u－u0,v－v0)(5．19)
F(u,v)exp－j2πx0uM+y0vNf(x－x0,y－y0)(5．20)


式(5．19)说明，给函数f(x,y)乘以一个指数项，就相当于把其变换后的傅里叶频谱在频率域进行平移。式(5．20)说明，给傅里叶频谱F(u,v)乘以一个指数项，就相当于把其反变换后得到的函数在空间域进行平移。
5．1．3图像的傅里叶频谱特性分析
1． 图像傅里叶频谱关于(M/2,N/2)的对称性
设f(x,y)是一幅大小为M×N的图像，根据离散傅里叶变换的周期性公式(5．16)，有

F(u,v)＝F(u+M,v+N)(5．21)


再根据式(5．21)和离散傅里叶变换的共轭对称性公式(5．18)就可得

F(u,v)=F(M－u,N－v)(5．22)


下面以图5．3为参照，进一步说明图像的傅里叶变换频谱关于(M/2,N/2)对称性。
根据式(5．22)，对于u=0有
当u=0、v=0时，F(0,0)=|F(M,N)|； 
当u=0、v=1时，F(0,1)=|F(M,N－1)|； 
当u=0、v=2时，F(0,2)=|F(M,N－2)|； 
……
当u=0、v=N/2时，|F（0，N/2）|=|F（M，N/2）|。
同理，对于v=0有 
当u=0、v=0时，F(0,0)=F(M,N)； 
当u=1、v=0时，F(1,0)=F(M－1,N)； 
当u=2、v=0时，F(2,0)=F(M－2,N)； 
……
当u=M/2、v=0时，|F（M/2，0）|=|F（M/2，N）|。
以上分析表明，如图5．3所示的傅里叶频谱的A区和D区关于坐标(M/2,N/2)对称。


图5．3频谱图对称性图示

更一般地，
当u∈［0，M/2］和v∈［0，N/2］时，有(M-u)∈［M/2，M］和(N-u)∈［N/2，N］； 
当u∈［M/2，M］和v∈［0，N/2］时，有(M-u)∈［0，M/2］和(N-u)∈［N/2，N］。
即傅里叶频谱的A区和D区关于坐标(M/2,N/2)对称； 傅里叶频谱的B区和C区关于坐标(M/2,N/2)对称。图5．1(b)和图5．2(b)即是原点坐标位于(0，0)的傅里叶频谱关于(M/2,N/2)对称的两个例子。
2．  原点坐标平移到(M/2,N/2)后的傅里叶频谱
由于图像的低频分量比较集中，在频谱图的4个角所占的区域都较小，所以不利于图像的频谱特性分析。根据傅里叶频谱的周期性和平移性，如图5．4(c)所示，当把傅里叶频谱图的原点从(0，0)平移至(M/2,N/2)时，以(M/2,N/2)为原点截取大小为M×N的区间，就可得到一个低频分量位于中心的图像频谱图，如图5．4(d)所示。


图5．4将坐标原点平移到(M/2,N/2)后的傅里叶频谱特征分析示意图


按照前面关于“图像傅里叶频谱关于(M/2,N/2)的对称性”的分析结论，即如图5．3所示，傅里叶频谱的A区和D区关于坐标(M/2,N/2)对称； 傅里叶频谱的B区和C区关于坐标(M/2,N/2)对称。把图5．1(b)的原点在(0，0)的傅里叶频谱平移到以(M/2,N/2)为原点，也就是把图5．1(b)按高(height)和宽(width)折半分成四块，然后把左上角的一块移到右下角位置，把右下角的一块移到左上角位置； 同理把右上角的一块移到左下角位置，把左下角的一块移到右上角位置，就可得到如图5．5(a)所示的、原点在(M/2,N/2)的傅里叶频谱。同理，由图5．2(b)中原点在(0，0)的傅里叶频谱可得如图5．5(b)所示的、原点在(M/2,N/2)的傅里叶频谱。


图5．5把傅里叶频谱坐标原点平移到(M/2,N/2)后的频谱图


显然，图5．5(a)和图5．5(b)的中心清楚地显示出了图像的低频分量的变化情况，不仅具有可视化特点，而且能简化图像的滤波过程等(详见5．3节和5．4节)。
有关傅里叶频谱的低频分量主要集中在中心的事实，将在5．3节的例5．1中对其进行进一步的说明。
3． 对图像进行傅里叶变换的意义
利用傅里叶变换把空间域的图像变换到频率域，就可得到该图像所含频率分量分布情况的傅里叶频谱，进而可以利用频率域的图像处理方法对傅里叶频谱进行处理，其意义在于以下几点。
(1) 简化计算。在空间域中处理图像时所进行的复杂的卷积运算，等同于在频率域中简单的乘积运算。
(2) 由于在用频谱图表示的频率域图像中，中心部位是能量集中的低频特征，反映的是图像的平滑部分； 随着不断远离频谱图的中心位置，对应于空间图像中变化越来越快的细节、边缘、结构复杂区域、突变部位和噪声等高频成分逐渐加强。所以，在频率域中滤波的概念更为直观，更容易理解； 即某些在空间域中难于处理或处理起来比较复杂的问题，在频率域却比较容易处理。
(3) 某些特定应用需求只能在频率域处理，如频率域图像特征提取、数据压缩、纹理分析、水印嵌入等。
5．1．4离散傅里叶变换的实现
1．  快速离散傅里叶变换的实现思路

在数字图像处理中，当M×N图像阵列的M和N较大时，直接利用离散傅里叶变换的定义式进行计算的计算量非常大，以至于在实际中是无法实现的。快速离散傅里叶变换算法的出现，才使得傅里叶变换用于实际的图像处理成为可能。截至目前已经出现了各种各样的快速离散傅里叶变换算法，比如按时间提取的基2一维快速离散傅里叶变换算法和按时间提取的基4一维快速离散傅里叶变换算法等。由于矩阵处理在空间资源上的要求比较高，所以这些快速算法都是一维的，对于二维图像的处理是分别通过按行和按列执行一维算法实现的。由于篇幅所限，本书略去了作者曾设计实现的“按时间提取的基2一维快速离散傅里叶变换算法”内容，感兴趣的读者可参阅参考文献［9］。
2．  串行计算二维离散傅里叶变换的方法
设f(x,y)是N×N的二维实序列，为表述方便，把f(x,y)看作是N×N的图像像素阵列，称

F(u,v)=1N∑N－1x=0∑N－1y=0f(x,y)WuxNWvyN(u,v=0,1,…,N-1)(5．23)

为图像像素矩阵的二维离散傅里叶变换(2DDFT)。其逆变换(2DIDFT)为

f(x,y)=1N∑N－1u=0∑N－1v=0F(u,v)W－uxNW－vyN(x,y=0,1,…,N-1)(5．24)

其中，WuxN=e－j2π/N·ux，WvyN=e－j2π/N·vy。WuxNWvyN是变换核。F(u,v)为图像的频谱。
根据二维傅里叶变换的可分离性，正变换公式(5．23)可改写成:

F(u,v)=1N∑N－1x=0∑N－1y=0f(x,y)WuyNWvxN(u,v=0,1,…,N-1)(5．25)


式（5.25）的表示形式说明，对于二维离散傅里叶变换，可先对图像像素矩阵的各列分别进行列傅里叶变换(简称列变换)，然后再对变换结果的各行分别进行行傅里叶变换(简称行变换)，这样就可以利用一维快速离散傅里叶变换算法串行计算二维离散傅里叶变换。其计算和变换过程如图5．6所示。


图5．6利用一维快速离散傅里叶变换算法串行计算二维离散傅里叶变换


为了简化程序,可把列变换后的结果进行转置，这样在进行行变换时就可应用列变换的程序，最后再把行变换后的结果进行一次转置即为变换结果。二维正变换的流程可简要描述为

f(x,y)→DFT列［f(x,y)］=F(x,v)转置FT(v,x)→DFT列［FT(v,x)］

=FT(v,u)转置F(u,v)

其中，DFT［f(x,y)］表示对f(x,y)进行傅里叶正变换。
图5．7是一个傅里叶变换的图例。其中图5．7(b)是对图5．7(a)的图像进行傅里叶正变换所得的频谱图像，图5．7(c)是进行傅里叶反变换恢复的图像。频率域图像处理过程是对正变换获得的频谱进行进一步处理，如低通滤波(去噪声)、高通滤波(图像增强)等，然后再对处理后的傅里叶进行反变换，获得处理后的空间域图像。


图5．7傅里叶变换图例


5．2频率域图像处理的基本实现思路
由5.1节可知，二维离散傅里叶变换很好地描述了二维离散信号的空间域与频域之间的关系，所以对于那些在空间域中表述起来比较困难、甚至是不太可能实现的图像处理问题，可以先通过对图像进行离散傅里叶变换，把图像变换到频率域，然后利用适当的频率域图像处理方式对图像进行处理，处理完后再把它变换回空间域中，从而解决那些在空间域不便于解决的图像处理问题。
5．2．1基本实现思想
由傅里叶频谱的特性可知，u和v同时为0时的频率成分对应于图像的平均灰度级。当从傅里叶频谱的原点离开时，低频对应着图像的缓慢变化分量，如一幅图像中较平坦的区域； 当进一步离开原点时，较高的频率开始对应图像中变化越来越快的灰度级，它们反映了一幅图像中物体的边缘和灰度级突发改变(如噪声)部分的图像成分。频率域图像增强正是基于这种机理，通过对图像的傅里叶频谱进行低通滤波(使低频通过，使高频衰减)来滤除噪声； 或通过对图像的傅里叶频谱进行高通滤波(使高频通过，使低频衰减)突出图像中的边缘和轮廓来实现图像的增强。
设f(x,y)为M×N的图像，F(u,v)为f(x,y)的原点在(0，0)的傅里叶频谱，F＇（u,v）
为原点平移到(M/2,N/2)的傅里叶频谱，H(u,v)为转移函数(也称为滤波函数)，G＇（u,v）
为对F＇（u,v）
进行频率域滤波后的输出结果，G(u,v)为原点平移回(0，0)的傅里叶频谱，g(x,y)为经频率域滤波后的输出图像，则有


G＇（u,v）=F＇（u,v）H(u,v)(5．26)


将G′（u,v）平移回原点（0,0）可得G（u,v），则有


g(x,y)=F－1［G(u,v)］(5．27)


在式(5．26)中，H和F＇
的相乘定义为二维函数逐元素的相乘，即H的第1个元素乘以F＇
的第1个元素，H的第2个元素乘以F＇
的第2个元素，以此类推。被滤波的图像可以由G(u,v)的傅里叶反变换(F－1)结果得到。通常情况下转移函数都为实函数，所以其傅里叶反变换的虚部为0。
综上所述，频率域图像处理的步骤如下。
(1) 对图像f(x,y)进行二维傅里叶变换，得到原点在(0，0)的傅里叶频谱F(u,v)。
(2) 把原点在(0，0)的傅里叶频谱F(u,v)按其高和宽折半分成四块，通过左上角块与右下角块的对称平移位置互换和右上角块与左下角块的对称平移位置互换，得到原点在(M/2,N/2)的傅里叶频谱F＇（u,v）
。
(3) 进行频率域滤波，即用设计的转移函数H(u,v)乘以原点在(M/2,N/2)的傅里叶频谱F＇（u,v）
，根据式(5．26)可得G＇(u,v)
。
(4) 将G＇(u,v)
的原点平移回(0，0)，可得G(u,v)。
(5) 对G(u,v)进行傅里叶反变换，即计算F－1［G(u,v)］，并取计算结果的实部，可得通过频率域滤波后的图像g(x,y)。
以上过程可简要地描述为图5．8。


图5．8频率域图像增强步骤


在图5．8中，前处理和后处理可能包括将输入图像向最接近的偶数维数转换(以便图像有合适的变换中心)、灰度级标定、输入向浮点数/双精度数的转换、输出向8比特整数格式的转换等。
5．2．2转移函数的设计
假设原图像f(x,y)经傅里叶变换为F(u,v)，频率域图像处理就是选择合适的滤波器函数H(u,v)对F(u,v)的频谱成分进行调整（让低频通过，使高频衰减； 让高频通过，使低频衰减； 等等），然后经傅里叶反变换得到空间域（频率域增强、频率域去噪等）的图像g(x,y)。因此，频率域图像处理的关键是设计合适的滤波转移函数H(u,v)。
关于转移函数的设计方法，一是先凭直观感觉选择一个理想的滤波器模型，然后通过反复的滤波实验和参数修正来逼近并设计出实际的滤波器； 二是利用频率成分和图像外观之间的对应关系选择频率域滤波器； 三是基于数学和统计准则设计频率域滤波器。
另外，对于大小为M×N的函数f(x,y)和h(x,y)，其卷积表示形式为

f(x,y)*h(x,y)=1MN∑Mm=0∑Nn=0f(m,n)h(x－m,y－n)(5．28)

其中，F(u,v)和H(u,v)分别为f(x,y)和h(x,y)的傅里叶变换。
根据傅里叶变换对

f(x,y)*h(x,y)F(u,v)H(u,v)(5．29)
F(u,v)*H(u,v)f(x,y)h(x,y)(5．30)

即空间域的卷积在频率域简化为相乘，频率域的卷积在空间域简化为相乘； 有时也可以将频率域的滤波器函数变换到空间域，然后在空间域对图像进行滤波运算。
一般来说，如果两个域中的滤波器尺寸（即滤波窗口大小）相同，通常在频率域中进行滤波计算更有效，空间域更适合于尺寸较小的滤波器。
5．3基于频率域的图像噪声消除——频率域低通滤波
图像中的噪声和边缘对应于傅里叶频谱的高频部分，选择能使低频通过、使高频衰减的转移函数H(u,v)，就可以根据式(5．26)和式(5．27)实现低通滤波，达到滤除噪声的目的。
下面分别介绍理想的低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器三种滤波器的有关概念和滤波原理。
5．3．1理想的低通滤波器
一个理想的低通滤波器的转移函数定义为

H(u,v)=1D(u,v)≤D0


0D(u,v)>D0(5．31)

其中，D0是一个非负整数，D(u,v)为频率平面从原点到点(u,v)的距离。
设已经将傅里叶频谱的原点平移到(M/2,N/2)，则点(u,v)到频率平面原点(M/2,N/2)的距离为

D(u,v)=［(u－M/2)2+(v－N/2)2］1/2(5．32)


理想的低通滤波器的含义为，在半径为D0的圆内的所有频率没有衰减地通过该滤波器，而在此圆之外的所有频率完全被衰减掉，所以D0称为截止频率。也就是说，对于低通滤波来说，当某个频率之上的幅度统统被忽略为0时，那个频率就称为截止频率。实现上，在进行低通滤波时，先要设定(或经多次实验得到)一个截止频率D0，设D0=π/3； 将高于π/3
的频率乘以0，让其全部被阻止； 将低于π/3的频率乘以1，让其全部通过。此时的频率值π/3
就称为截止频率。
图5．9(a)给出了理想低通滤波器的转移函数H的横截面图，其中横轴用于表示离原点的径向距离。通过将横截面绕原点旋转360°即可得到完整的理想低通滤波器转移函数，即图5．9(b)所示的转移函数H的透视图。该透视图的含义是： 只有那些位于该圆柱体内的频率范围的信号才能通过，而位于圆柱体外的频率成分都将被滤除。


图5．9理想低通滤波器的转移函数横截面图和透视图


在图5．9(a)及本节的以下内容中，均假设D0的频率平面原点为(M/2,N/2)。
需要说明的是，理想低通滤波器的数学意义是十分清楚的，利用计算机对其进行模拟也是可行的，但在实际中用电子元件实现直上直下的理想低通滤波器是不可能的，所以才将其称为“理想”低通滤波器。
【例5．1】频率域理想低通滤波器的滤波效果及低频特性分析。
解： 若一般地设R为截止频率的圆周半径，EB为圆周内能量(图像功率)与原图像总能量(总功率)的百分比，根据图像信号能量在频率域上的分布有

EB=∑u∈R∑v∈RP(u,v)∑M－1u=0∑N－1v=0P(u,v)×100%(5．33)


图5．10(a)是一幅包含了全部细节的原始图像； 图5．10 (b)是它的傅里叶频谱图，利用其截止频率半径D分别为10、20、40和80确定的理想低通滤波器对原图像进行低通滤波，所得的图像分别为图5．10(c)~图5．10(f)。根据式(5．33)计算可知，图5．10(c)~图5．10(f)分别包含了原图像中95．5%、97．9%、99．0%和99．6%的能量。


图5．10频率域理想低通滤波器的滤波效果及低频特性分析


图5．10(c)~图5．10(f)的结果说明:
 傅里叶频谱图的低频分量主要集中在中心; 
 指明了以截止频率为半径的圆内的图像功率与图像总功率的量级关系； 
 说明了高频成分对于表现图像的轮廓和细节是十分重要的。
图5．10(c)仅滤除掉占总能量的4．5%的高频分量，图像就变得十分模糊，并有明显的振铃效应； 图5．10(d)仅滤除掉占总能量的2．1%的高频分量，图像仍存在着一定程度的模糊和振铃效应； 当图5．10(e)仅滤除掉占总能量的1．0%的高频分量时，图像视觉效果才变得尚可接受； 而当图5．10(f)仅滤除掉占总能量的0．4%的高频分量时，图像才变得与原图像几乎一样。
5．3．2巴特沃斯低通滤波器
一个n阶的巴特沃斯(Butterworth)低通滤波器的转移函数定义为

H(u,v)=11+［D(u,v)/D0］2n(5．34)

其中，D0为截止频率，D(u,v)是频率平面从原点到点(u,v)的距离，如式(5．32)所示。
图5．11(a)给出了阶数分别为1、2和3的巴特沃斯低通滤波器的转移函数H的横截面图。巴特沃斯低通滤波器的转移函数H的透视图如图5．11(b)所示，该透视图的含义是： 只有那些位于该草帽形体内的频率范围的信号才能通过，而位于草帽形体外的频率成分都将被滤除掉。由图中可见，巴特沃斯低通滤波器在高低频率间的过渡比较平滑。


图5．11巴特沃斯低通滤波器的转移函数横截面图和透视图


一般情况下取当H的最大值降到某个百分比时对应的频率为截止频率。在式(5．34)中，当D(u,v)=D0时，H(u,v)=0．5(从最大值1降到50％)。一个常用的方法是取H降到最大值的1/2时的频率为截止频率。
图5．12是利用巴特沃斯低通滤波器进行图像去噪的实验结果。图5．12(b)是加入了噪声密度为0．02的椒盐噪声图像，图5．12(c)~图5．12(h)取不同的截止频率(10和30)和不同的滤波器阶数(1、2和3)的去噪效果图像。分析可知，滤波器阶数n的值越大，去噪结果图像越模糊(去噪效果越差)； 截止频率D0的值越大，去噪结果图像越清晰(去噪效果越好)。


图5．12利用巴特沃斯低通滤波器进行图像去噪实验结果


5．3．3高斯低通滤波器
由于高斯函数的傅里叶变换和反变换均为高斯函数，并常常用来帮助寻找空间域与频率域之间的联系，所以基于高斯函数的滤波具有特殊的重要意义。
一个二维的高斯低通滤波器的转移函数定义为

H(u,v)=e－D2(u,v)/2σ2(5．35)

其中，D(u,v)是频率平面从原点到点(u,v)的距离，如式(5．32)所示。σ表示高斯曲线扩展的程度，当σ=D0时，可得到高斯低通滤波器的一种更为标准的表示形式

H(u,v)=e－D2(u,v)/2D20(5．36)

其中，D0是截止频率； 当D(u,v)=D0时，H下降到它的最大值的0．607处。
图5．13(a)给出了D0=10、D0=20和D0=40的高斯低通滤波器的转移函数H的横截面图。高斯低通滤波器的转移函数H的透视图如图5．13(b)所示。该透视图的含义是： 只有那些位于该草帽形体内的频率范围的信号才能通过，而位于草帽形体外的频率成分都将被滤除。


图5．13高斯低通滤波器的转移函数横截面图和透视图


与巴特沃斯低通滤波器相比，高斯低通滤波器没有振铃现象。另外，在需要严格控制低频和高频之间截止频率的过渡的情况下，选择高斯低通滤波器更合适一些。
在本节结束前需要说明的是，在频率域中，滤波器越窄，滤除掉的高频成分就越多，滤波后的图像就越模糊。这一特性正好对应于在空间域中，滤波器越宽(模板尺寸越大)，平滑后的图像就越模糊的现象。
5．4基于频率域的图像增强——频率域高通滤波
图像中的边缘和灰度的陡峭变化对应于傅里叶频谱的高频部分，选择能使高频通过、使低频衰减的转移函数H(u,v)，就可以根据式(5．26)和式(5．27)实现高通滤波，达到突出图像的高频边缘成分，实现图像增强的效果。
5．4．1理想的高通滤波器
一个理想的频率域高通滤波器的转移函数定义为

H(u,v)=0D(u,v)≤D0


1D(u,v)>D0(5．37)

其中，D0为截止频率，D(u,v)为频率平面从原点到点(u,v)的距离，如式(5．32)所示。对于高通滤波来说，当某个频率之下的幅度统统被忽略为0时，那个频率就称为截止频率。
理想高通滤波器的含义为，将以半径为D0的圆周内的所有频率置零，而让圆周外的所有频率毫不衰减地通过。
图5．14(a)给出了理想高通滤波器的转移函数H的横截面图，图5．14(b)给出了转移函数H的透视图。该透视图的含义是： 只有那些位于该圆柱形体外的频率范围的信号才能通过，而位于圆柱形体内的频率成分都将被滤除。


图5．14理想高通滤波器的转移函数横截面图和透视图


与理想低通滤波器一样，在实际中用电子元件实现直上直下的理想高通滤波器是不现实的，所以才将其称为“理想”高通滤波器。
5．4．2巴特沃斯高通滤波器
一个n阶的巴特沃斯高通滤波器的转移函数定义为

H(u,v)=11+［D0/D(u,v)］2n(5．38)

其中，D0为截止频率，D(u,v)是频率平面从原点到点(u,v)的距离，如式(5．32)所示。
图5．15(a)给出了阶数为1的巴特沃斯高通滤波器的转移函数H的横截面图。巴特沃斯高通滤波器的转移函数H的透视图如图5．15(b)所示，该透视图的含义是： 只有那些位于该倒立形草帽体外的频率范围的信号才能通过，而位于倒立形草帽体内的频率成分都将被滤除。与巴特沃斯低通滤波器一样，巴特沃斯高通滤波器在高低频率间的过渡比较平滑。


图5．15巴特沃斯高通滤波器的转移函数横截面图和透视图


与巴特沃斯低通滤波器一样，一般情况下取当H的最大值降到某个百分比时对应的频率为截止频率。

5．4．3高斯高通滤波器
一个截止频率距离原点为D0的高斯高通滤波器的转移函数定义为

H(u,v)=1-e－D2(u,v)/2D20(5．39)

其中，D(u,v)是频率平面从原点到点(u,v)的距离，如式(5．32)所示。
图5．16(a)给出了典型的高斯高通滤波器的转移函数H的横截面图，图5．16(b)给出了该高斯高通滤波器转移函数H的透视图。该透视图的含义是： 只有那些位于该倒立形草帽体外的频率范围的信号才能通过，而位于倒立形草帽体内的频率成分都将被滤除。


图5．16高斯高通滤波器的转移函数横截面图和透视图


图5．17给出了用高斯高通滤波器进行图像增强时取不同截止频率时的实验结果。分析可知： 随着截止频率D0值的增大，增强效果越来越明显，即使对于细线条的头发丝，肉眼可以看出其细微的增强效果变化。但同时，随着D0的增大，低频信息减少较多，图像中灰度剧烈变化的区域也出现了微弱的振铃现象(即图像灰度剧烈变化处产生的振荡，就好像钟被敲击后产生的空气振荡)。


图5．17利用高斯高通滤波器进行图像增强的实验结果


5．5带阻滤波和带通滤波
带阻滤波器(bandpass filter)与带通滤波器(bandstop filter)用于对某些区域的某一频率范围内的频率分量抑制其通过或让其通过。
5．5．1带阻滤波器
在某些应用中，图像的质量可能受到带有一定规律的结构性噪声的影响。比如，图像上叠加有正弦干扰图案就是这类噪声的一个典型情况。当正弦干扰图案比较明显时，会在图像的频谱平面上出现两个比较明显的对称点(由于傅里叶变换的共轭对称性所致)。这种用于消除以某点为对称中心的给定区域内的频率，或用于阻止以原点为对称中心的一定频率范围内信号通过的问题，就可以用带阻滤波器实现。
一个用于消除以某点为对称中心、以D为半径的圆上的带阻滤波器，可以通过把以原点为中心的高通滤波器平移到该点得到，设该带阻滤波器的中心为点(u0,v0)，半径为D0，则其传递函数定义为

H(u,v)=0D(u,v)≤D0


1D(u,v)>D0(5．40)

其中

D(u,v)=(u－u0)2+(v－v0)212(5．41)


由于傅里叶变换的共轭对称性，要求带阻滤波器必须成对出现，所以一个用于消除以(u0,v0)为中心、以D0为半径的对称区域内的所有频率的理想带阻滤波器的转移函数定义为

H(u,v)=0D1(u,v)≤D0或D2(u,v)≤D0


1其他(5．42)

其中

D1(u,v)=［(u－u0)2+(v－v0)2］1/2(5．43)
D2(u,v)=［(u+u0)2+(v+v0)2］1/2(5．44)


利用上述构建带阻滤波器的思路，还可以把前面所讲的几种高通滤波器转变为带阻滤波器。例如一种n阶径向对称的巴特沃斯带阻滤波器的传递函数可定义为

H(u,v)=11+D(u,v)WD2(u,v)－D202n(5．45)

其中，W为阻带带宽，D0为阻带中心半径。
图5．18给出了一个典型的带阻滤波器的转移函数H的透视图。该透视图的含义是： 只有那些位于两个圆柱体外的频率范围的信号才能通过，而位于两个圆柱体内的频率成分都将被滤除。


图5．18理想带阻滤波器的转移函数的透视图


5．5．2带通滤波器
带通滤波器与带阻滤波器相反，它允许以原点为对称中心的一定频率范围内的信号通过，而使其他频率衰减或抑制。理想的带通滤波器的转移函数可定义为

H(u,v)=1D1(u,v)≤D0或D2(u,v)≤D0


0其他(5．46)


带通滤波器也可以通过对相应的带阻滤波器进行“翻转”获得。若设H′(u,v)为带阻滤波器的传递函数，则对应的带通滤波器的传递函数H(u,v)可定义为

H(u,v)=1－H′(u,v)(5．47)


图5．19是一个典型的带通滤波器的转移函数H的透视图。该透视图的含义是： 只有那些位于两个圆柱体内的频率范围的信号会通过，而位于两个圆柱体外的频率成分都将被滤除。


图5．19理想带通滤波器的转移函数的透视图


习题5
5．1解释下列术语。
(1) 频率域图像增强(2) 低通滤波
(3) 高通滤波(4) 截止频率
5．2已知有数学关系式exp(x)=ex和eix=cosx+isinx，请写出下列公式的展开式。
(1) exp-j2π(xu+yv)N
(2) expj2π(xu+yv)N
5．3一般来说，进行图像傅里叶变换的目的是什么？
5．4设N×N的图像f(x,y)的傅里叶变换为F(u,v)，请写出反映该图像的灰度平均值与傅里叶变换的直流(DC)分量的关系的公式。
5．5回答下列问题： 
(1) 对图像f(x,y)进行傅里叶变换后得到的频谱图，低频位于其中心处还是四角处？
(2) 当把对图像f(x,y)进行傅里叶变换后得到的频谱图的坐标原点移到中心对称时，低频位于其中心处还是四角处？
(3) 频率域低通滤波实现的图像处理功能是什么？
(4) 频率域高通滤波实现的图像处理功能是什么？
5．6编写一个利用巴特沃斯低通滤波器进行图像去噪的MATLAB程序。
5．7编写一个利用高斯高通滤波器进行图像增强的MATLAB程序。