第二章傅里叶变换

1、数字图像处理技术数字图像处理技术Digital Image Processing Technique吴昊天吴昊天E-mail: 2第二章第二章 图像变换技术图像变换技术要点：要点：1. 主要介绍图像处理重要的工具主要介绍图像处理重要的工具 傅里叶变换傅里叶变换.2. 傅立叶变换在图象处理中的意义是什么傅立叶变换在图象处理中的意义是什么?3. 什么是高频、中频和低频成分，它们分别对应空间域图像的哪些部分什么是高频、中频和低频成分，它们分别对应空间域图像的哪些部分?4. 什么是卷积定理，它在图象处理中的作用是什么什么是卷积定理，它在图象处理中的作用是什么?5. 傅立叶变换的性质。傅立叶变换的性质。

2、31822年，傅立叶（年，傅立叶（Fourier）发表了发表了“热传导解析理论热传导解析理论”，提出了傅，提出了傅立叶变换。它本质上提出了一种与空间思维不同的频域思维方法。立叶变换。它本质上提出了一种与空间思维不同的频域思维方法。傅立叶变换是十九世纪数学界和工程界最辉煌的成果之一。它一直是傅立叶变换是十九世纪数学界和工程界最辉煌的成果之一。它一直是信号处理领域中最完美、应用最广泛、效果最好的一种分析手段。它信号处理领域中最完美、应用最广泛、效果最好的一种分析手段。它也是线性系统分析的有利工具。也是线性系统分析的有利工具。傅立叶变换能使我们从空间域（或时域）与频率域两个不同的角度来傅立叶变换能使

3、我们从空间域（或时域）与频率域两个不同的角度来看待信号或图象的问题。有时在时域无法解决的问题，在频域却是显看待信号或图象的问题。有时在时域无法解决的问题，在频域却是显而易见的。而易见的。 4n傅里叶分析中最重要的结论就是几乎傅里叶分析中最重要的结论就是几乎“所有所有”的函数的函数(信号信号)都可以表都可以表示为示为(分解成分解成)简单的简单的(加权加权)正弦波和余弦波之和。从而提供了一种具正弦波和余弦波之和。从而提供了一种具有物理意义的函数表达方式。有物理意义的函数表达方式。设设: f(x)是以是以T为周期的函数为周期的函数, 满足一定的条件满足一定的条件, 例如绝对可积例如绝对可积, 则有则

4、有xTjkkkeaxf2 )(TxTjkkdxexfTa02)(156特别注意特别注意n傅傅氏级数中基底的物理意义非常明确氏级数中基底的物理意义非常明确, 每一个基函数都是一个单频谐每一个基函数都是一个单频谐波波, 而相应的系数而相应的系数(频谱频谱)表明了原函数对这种频率成份贡献的大小表明了原函数对这种频率成份贡献的大小(原原函数在这个谐波上的投影函数在这个谐波上的投影), 或者说原函数中某种频率成分的多少或者说原函数中某种频率成分的多少.n从图像从图像(信号信号)处理的角度处理的角度, 利用谐波的物理性质可以通过对系数的处理利用谐波的物理性质可以通过对系数的处理达到对图像的处理达到对图像的

5、处理, 如增强、压缩等等如增强、压缩等等.nf(x)傅傅氏系数氏系数ak的计算的计算, 需要用到函数在整个空间需要用到函数在整个空间(或时间或时间)上的分布情上的分布情况况.TxTjkkdxexfTa02)(17一维傅里叶变换及反变换一维傅里叶变换及反变换考虑定义在无穷区间连续函数的傅里叶变换公式考虑定义在无穷区间连续函数的傅里叶变换公式(通常函数要满足一通常函数要满足一定的条件才能保证傅里叶变换的存在性和收敛性定的条件才能保证傅里叶变换的存在性和收敛性):dxexfuFuxj2)()(dueuFxfuxj2)()(8连续情形的傅里叶变换比较方便用于公式推导和定理证明连续情形的傅里叶变换比较方

6、便用于公式推导和定理证明, 但在实际但在实际应用中应用中, 面临更多也更实际的是离散的情况面临更多也更实际的是离散的情况. 定义离散情形的傅里叶变换定义离散情形的傅里叶变换(DFT)公式公式: f(x)为离散函数为离散函数, 其中其中x = 0, 1, , M-1.1, 1 , 0 )(1)(10/2MuexfMuFMxMuxj1, 1 , 0 )()(10/2MxeuFxfMxMuxj离散傅里叶变换和它对应的反变换总是存在的离散傅里叶变换和它对应的反变换总是存在的, 不必特地关心分析不必特地关心分析各项的意义各项的意义.9频率域的概念频率域的概念:利用欧拉公式利用欧拉公式: ej = cos

7、 + jsin , 有有101( )( )cos(2/)sin(2/)MxF uf xux Mjux MM其中其中u = 0, 1, 2, , M-1.变量变量u确定了变换的频率成分确定了变换的频率成分 u的取值范围称为频率域的取值范围称为频率域(给定一个给定一个u 上述公式可以计算出离散信号中包含了上述公式可以计算出离散信号中包含了“多少多少”这个频率的谐波这个频率的谐波).对每一对每一个个u, F(u)称为变换的频率分量称为变换的频率分量(也叫振幅也叫振幅). F(u)可以看作f(x)在谐波上的投影，即f(x)在频率为u的谐波上占有的成份。10谱的概念谱的概念:注意到傅里叶变换后的函数是在

8、复数域内注意到傅里叶变换后的函数是在复数域内, 也可以表示为也可以表示为F(u) = R(u) + iI(u)或极坐标的形式或极坐标的形式: F(u) = |F(u)|ej (u).我们把量我们把量|F(u)| = R2(u) + I2(u)1/2称为傅里叶变换的幅度称为傅里叶变换的幅度(Magnitude)或者谱或者谱(Spectrum). 这是在图像处理中要经常用到的量这是在图像处理中要经常用到的量. 谱可以表示原函谱可以表示原函数数(或图像或图像)对某一频谱分量的贡献对某一频谱分量的贡献.)()(arctan)(uRuIu称为变换的相角或者相位谱称为变换的相角或者相位谱, 用来表示原函数

9、中某一频谱分量的起用来表示原函数中某一频谱分量的起始位置始位置*.另外另外, 一个重要的量是功率谱一个重要的量是功率谱(有时也叫能量谱、谱密度有时也叫能量谱、谱密度)P(u) = |F(u)|2 = R2(u) + I2(u)11例例4.1 两个简单一维函数的傅里叶谱两个简单一维函数的傅里叶谱特征特征: (1) 当曲线下的面积在当曲线下的面积在x域加倍时域加倍时, 频率谱的高度也加倍频率谱的高度也加倍;(2) 当函数的长度加倍时当函数的长度加倍时, 相同长度区域内的零点数量也加倍相同长度区域内的零点数量也加倍. 极限情况极限情况? 12解释解释: 图图a函数的傅里叶变换为函数的傅里叶变换为:1

10、010/2)(1)(KxxuKxMuxjrMAAeMuF)(/2Mjer易见易见, 当当u = 0时时, ru = 1, 故而故而MAKMAuFKx 1)(10若若u 0, 则则ru 1, 对对u = 1, 2, , M-1,uuKKxxrrMArMAuF11 )()(10u利用欧拉公式可知利用欧拉公式可知: r = cos(2 /M) - jsin(2 /M).所以所以, 当当uK是是M的倍数时的倍数时, 就有就有ruK = 1(当然这时也有当然这时也有ru2K = 1 ), 从而从而F(u) = 0. 如果图如果图a中函数中函数f(x)非零点的个数是非零点的个数是K时时, F(u) = 0

11、的点数是的点数是n个个, 那么那么, 当当f(x)的非零点数是的非零点数是2K时时, F(u) = 0的点数应该是的点数应该是2n个个. 13关于变量的说明关于变量的说明:书中的记号书中的记号 f(x)(x = 0, 1, , M-1)表示从连续函数中取表示从连续函数中取M个样点个样点, 这些这些点不一定选取为区间点不一定选取为区间0, M-1中的整数点中的整数点. 通常用通常用x0(任意位置的任意位置的)表示第一表示第一个取样点个取样点, x是取样间隔是取样间隔. 所以所以, f(x)理解为理解为其中其中: u = 0, 1, , M-1.值得注意的是值得注意的是, 当当M固定时固定时, x

12、和和 u之间有如下的反比关系之间有如下的反比关系:)( )(0 xxxfxf其中其中: x = 0, 1, , M-1.同理同理, 变量变量u有相似的解释有相似的解释, 但序列通常总是从但序列通常总是从0频率开始频率开始. 因此因此, u的取的取值序列为值序列为u = 0, u, 2 u, , (M-1) u. F(u)理解为理解为:)( )(uuFuFxMu114二维二维DFT及反变换及反变换 离散情形离散情形完全类似完全类似.设设f(x, y)是一幅尺寸为是一幅尺寸为MN的图象函数的图象函数, 相应的相应的离散傅里叶变换离散傅里叶变换(DFT)可以表示为可以表示为: dxdyeyxfvuF

13、vyuxj)(2),(),( dudvevuFyxfvyuxj)(2),(),(1010)/(2),(1),(MxNyNvyMuxjeyxfMNvuF傅里叶谱傅里叶谱: |F(u, v)| = R2(u, v) + I2(u, v)1/2二维傅里叶变换本质上是一维情形向两个方向的简单扩展二维傅里叶变换本质上是一维情形向两个方向的简单扩展. 1010)/(2),(),(MuNvNvyMuxjevuFyxf15相角相角: (u, v) = arctan(I(u, v)/R(u, v)功率谱功率谱: P(u, v) = |F(u, v)|2 = R2(u, v) + I2(u, v)注意注意:为了把

14、变换后的中心移到图像的中心为了把变换后的中心移到图像的中心(M/2, N/2), 通常在变换通常在变换之前都要在函数上乘以之前都要在函数上乘以(-1)x+y. 这是由于这是由于傅傅里叶变换的平移性质决定的里叶变换的平移性质决定的, 以一维为例以一维为例:f(x)ej2 u0 x/M F(u - u0)f(x - x0) F(u)e-j2 x0u/M当当x0 = M/2或或u0 = M/2时时f(x)(-1)x F(u - M/2)f(x - M/2) F(u)(-1)u16对频谱图像的认识对频谱图像的认识?171010),(1)0 , 0(MxNyyxfMNF易见易见是图像的平均灰度是图像的平

15、均灰度. 因为在原点处两个方向的频率都为零因为在原点处两个方向的频率都为零, 所以所以, 这这个量经常被称为频谱的直流分量个量经常被称为频谱的直流分量(DC). DC就是电子工程领域中的直流就是电子工程领域中的直流, 也就是零频率的电流也就是零频率的电流.和一维情形和一维情形, 同样也有同样也有xMu1yNv118例例4.2 一个简单函数的频谱一个简单函数的频谱(已经做过中心化处理已经做过中心化处理).图像是图像是512 512的黑色背景上叠加一个的黑色背景上叠加一个20 40 象素的白色矩形象素的白色矩形. 频谱的显示频谱的显示经过了对数变换处理以加强灰度级细节经过了对数变换处理以加强灰度级

16、细节, 并适当调整了灰度强度并适当调整了灰度强度.可以看出可以看出, u方方向谱的零点分隔恰好是向谱的零点分隔恰好是v方向零点分隔的两倍方向零点分隔的两倍, 在不同方向上符合了原图中在不同方向上符合了原图中1:2的的矩形尺寸比例矩形尺寸比例. 这和一维情形完全类似这和一维情形完全类似.极限情况、能量分布情况极限情况、能量分布情况? 19频率域滤波频率域滤波频率域的基本性质频率域的基本性质傅里叶变换每个频谱分量都要涉及到图像空间中的每个傅里叶变换每个频谱分量都要涉及到图像空间中的每个像像素素, 所以所以一般来说一般来说, 频谱信息中很难看出空间的信息频谱信息中很难看出空间的信息。但由于频率反映的

17、是空间但由于频率反映的是空间强度的变化率强度的变化率, 如低频对应着图像变化慢的部分如低频对应着图像变化慢的部分, 高频对应着图像变化快高频对应着图像变化快的部分的部分。所以所以, 在某种意义上两者之间仍然有不可分割的联系在某种意义上两者之间仍然有不可分割的联系, 尽管这些尽管这些联系是联系是“总体总体”的的.1010)/(2),(1),(MxNyNvyMuxjeyxfMNvuF 1010)/(2),(),(MuNvNvyMuxjevuFyxf(4.2.16)(4.2.17)20例例4.3 一幅图像和显示某些重要特征的傅利叶谱一幅图像和显示某些重要特征的傅利叶谱集成电路的扫描电子显微镜图像集成

18、电路的扫描电子显微镜图像, 放大放大2500倍倍.注意注意: 45 角的角的两个强边缘和热感应两个强边缘和热感应不足引起的两个白色不足引起的两个白色氧化突起氧化突起.高频和低频部分高频和低频部分, 能量分布的一般情况能量分布的一般情况.21频率域滤波基础频率域滤波基础 f(x, y) F(u, v)步骤步骤:(1) 用用(-1)x+y乘以输入图像乘以输入图像, 做频谱中心化处理做频谱中心化处理;(2) 计算计算(1)结果的结果的DFT, 即即F(u, v);(3) 用滤波器函数用滤波器函数H(u, v)乘以乘以F(u, v)(在频谱域处理图像在频谱域处理图像)滤波器函数滤波器函数下面讨论下面讨

19、论;(4) 计算计算(3)中结果的反中结果的反DFT;(5) 得到得到(4)结果中的实部结果中的实部;(6) 用用(-1)x+y乘以乘以(5)中的结果中的结果滤波和滤波器滤波和滤波器: 滤波顾名思义就是阻止或减少信号或图像中的某些频滤波顾名思义就是阻止或减少信号或图像中的某些频率成分率成分. 滤波器滤波器(函数函数)就是能起到这样作用的函数就是能起到这样作用的函数. 一般表达式一般表达式:G(u, v) = H(u, v)F(u, v)滤波后的结果图像可以从滤波后的结果图像可以从G(u, v)的反傅里叶变换得到的反傅里叶变换得到.g(x, y) = -1G(u, v)22频域滤波的基本步骤频域

20、滤波的基本步骤(包括前处理和后处理包括前处理和后处理):陷波滤波器陷波滤波器(Notch Filter): 使得处理后的图像平均值为零使得处理后的图像平均值为零, 从而图从而图像的整体灰度降低像的整体灰度降低.others 1)2/ , 2/(),( 0),(NMvuvuH23对图对图4.4a使用陷波滤波器使用陷波滤波器, 将傅里叶变换的原点设置为将傅里叶变换的原点设置为0.24低通滤波器和高通滤波器低通滤波器和高通滤波器25高频增强滤波高频增强滤波26空间域滤波器和频率域滤波器之间的对应关系空间域滤波器和频率域滤波器之间的对应关系 结论结论: 空间域的滤波器空间域的滤波器, 和频率域的滤波器

21、组成了一组傅里叶变换对和频率域的滤波器组成了一组傅里叶变换对.也就是说也就是说, 给出在频率域的滤波器给出在频率域的滤波器, 通过傅里叶反变换就可以得到在空间通过傅里叶反变换就可以得到在空间域相应的滤波器域相应的滤波器, 反之亦然反之亦然. 这个最基本的联系是由著名的卷积定理建立起来的这个最基本的联系是由著名的卷积定理建立起来的. 两个两个M N的离散函数的离散函数f(x, y)和和h(x, y)的卷积定义为的卷积定义为1010),(),(),(),(MmNnnymxhnmfyxhyxf空间域滤波的本质上就是用选好的掩模空间域滤波的本质上就是用选好的掩模(mask), 经过一定的处理经过一定的

22、处理, 与与给定的图像作卷积给定的图像作卷积.例如例如: 平滑滤波的掩模平滑滤波的掩模27卷积定理卷积定理F(u, v) f(x, y)的傅里叶变换的傅里叶变换H(u, v) h(x, y)的傅里叶变换的傅里叶变换则有则有f(x, y)*h(x, y) F(u, v)H(u, v)f(x, y)h(x, y) F(u, v)*H(u, v)28卷积定理的证明提示：卷积定理的证明提示： （51） （52） 证明：由定义： )()()()(sGsFtgtfF)()()()(tgtfsGsF1F)()()()()()()()()()()(2)(22)(22sGsFutvduufdvevgdudteu

23、tgufdtdueutguftgtfsujvsjsujutsjstj 令F29定义定义: 在坐标在坐标(x0, y0)处强度为处强度为A的冲激的冲激(脉冲脉冲)函数函数A (x x0, y y0)满足满足一些结论：一些结论：),(),(),(00101000yxAsyyxxAyxsMxNy1。)0 , 0(),(),(1010AsyxAyxsMxNy2。MNeyxMNvuFNvyMuxjMxNy1),(1),()/(21010 3。),(1),(),(1),(*),(1010yxhMNnymxhnmMNyxhyxfMmNn 30这个结论的指导意义这个结论的指导意义: 在频率域选择滤波器更为直观

24、在频率域选择滤波器更为直观, 物理意义比较物理意义比较明确明确. 在空间域用较小的模板进行滤波计算则比较经济在空间域用较小的模板进行滤波计算则比较经济. 如果两个滤波器如果两个滤波器的大小一样的大小一样, 则在频率域进行滤波计算比较方便则在频率域进行滤波计算比较方便. 更令人感兴趣的是更令人感兴趣的是, 在在频率域找到一个有意义的滤波器频率域找到一个有意义的滤波器, 然后作傅里叶反变换然后作傅里叶反变换, 并以此为依据并以此为依据, 尽可能在空间域建造尽可能小的滤波模板尽可能在空间域建造尽可能小的滤波模板.高斯滤波器高斯滤波器以一维为例以一维为例, 设高斯滤波器设高斯滤波器(频率域频率域)为为

25、H(u) = Ae-u2/2 2, 其中其中: 为高斯为高斯曲线的标准差曲线的标准差. 而相关的空间滤波器则为而相关的空间滤波器则为22222)(xAexh注意高斯函数有两个重要特性注意高斯函数有两个重要特性: 它本身和它的傅里叶变换都是实高斯函数它本身和它的傅里叶变换都是实高斯函数; 这一对傅利叶变换相互间有某种反比特性这一对傅利叶变换相互间有某种反比特性, 即当即当H(u)有较宽的轮有较宽的轮廓廓(大的大的 值值)时时, h(x)有较窄的轮廓有较窄的轮廓, 反之亦然反之亦然.当接近无限时当接近无限时, H(u)趋于常数量趋于常数量, h(x)则趋于冲激函数则趋于冲激函数.31高斯函数是低高

26、斯函数是低通滤波器通滤波器(?), 但也但也可以通过组合构造可以通过组合构造出更复杂的滤波器出更复杂的滤波器, 如高通和带通等如高通和带通等. 2222122/2/)(uuBeAeuH22222212222122)(xxBeAexh高频增强高频增强空间域空间域频域滤波器频域滤波器空间域空间域其中其中:21B，A32在频率域技术的发展过程中在频率域技术的发展过程中, 经常被问到的问题就是计算复杂性经常被问到的问题就是计算复杂性. 为为什么在空间域用小模板就可以做到的事情什么在空间域用小模板就可以做到的事情, 却要在频率域里去考虑却要在频率域里去考虑? 这个这个问题可以从两方面来回答问题可以从两方

27、面来回答:一方面是因为频率域直观一方面是因为频率域直观, 可以很容易凭经验指定滤波器可以很容易凭经验指定滤波器;另一方面的答案取决于空间模板的大小另一方面的答案取决于空间模板的大小, 并可以用对比计算复杂性的并可以用对比计算复杂性的方式来回答方式来回答.一个用于计算复杂性对比的基准是计算卷积一个用于计算复杂性对比的基准是计算卷积. 空间域的卷积由公式空间域的卷积由公式(4.2.30)给出给出. 而由卷积定理可知而由卷积定理可知, 对两个函数傅里叶变换的乘积做反变换对两个函数傅里叶变换的乘积做反变换, 也可以得到同样的结果也可以得到同样的结果.假定假定: 在同一台微机上用软件实现这两种算法在同一

28、台微机上用软件实现这两种算法(频率域使用频率域使用4.6.6节讨节讨论的论的FFT), 会发现对较小的会发现对较小的M和和N值值, 频率域的计算会非常快频率域的计算会非常快.频率域可以看成是一个实验室频率域可以看成是一个实验室, 从中可以利用频率成分和图像外观之从中可以利用频率成分和图像外观之间的对应关系间的对应关系.后面的许多例子将反复证明后面的许多例子将反复证明: 一些在空间域中很难表述、甚至不可能一些在空间域中很难表述、甚至不可能的增强任务的增强任务, 在频率域会变得非常简单在频率域会变得非常简单.例如周期性背景噪声例如周期性背景噪声.33这一节是本章的重要内容之一这一节是本章的重要内容

29、之一. 目的是解决实际计算傅里叶变换时目的是解决实际计算傅里叶变换时出现的问题出现的问题, 包括快速傅里叶变换包括快速傅里叶变换. 假设假设: f(x, y)是是M N的二维图像的二维图像, F(u, v)是其傅里叶变换是其傅里叶变换, 有有1010)/(2),(1),(MxNyNvyMuxjeyxfMNvuF1010)/(2),(),(MuNvNvyMuxjevuFyxf平移性质平移性质f(x, y)ej2 (u0 x/M+v0y/N) F(u u0, v v0)f(x x0, y y0) F(u, v)e-j2 (x0u/M+x0v/N)二维傅里叶变换的性质二维傅里叶变换的性质34注意注意

30、: 当当u0 = M/2和和v0 = N/2时时, 有有ej2 (u0 x/M+v0y/N) = ej (x+y) = (-1)x+yf(x, y)(-1)x+y F(u M/2, v N/2)提供了频谱中心化的方法提供了频谱中心化的方法.分配性和比例变换性分配性和比例变换性旋转性质旋转性质: 如果引入极坐标如果引入极坐标:x = rcos ,y = rsin ,u = cos , v = sin f(r, + 0) F( , + 0)此式表明当原图像此式表明当原图像f(x, y)以某一角度旋转时以某一角度旋转时, 其谱平面也将转过同样其谱平面也将转过同样的角度的角度. 注注: 对离散情形对离

31、散情形, 若有若有M = N, 证明易见证明易见. ),(),(),(),(2121yxfbyxfayxbfyxaf)/,/(1),(bvauFabbyaxff(x, y)ej2 (u0 x/M+v0y/N) F(u u0, v v0)352021-7-2635(a)原始图像 (b) DFT变换 (c) 原始图像旋转45 (d) 旋转之后DFT变换结果 36*周期性和对称性周期性和对称性直接验证可知直接验证可知, 离散傅里叶变换具有周期性离散傅里叶变换具有周期性F(u, v) = F(u+M, v) = F(u, v+N) = F(u+M, v+N)注意注意e-j2 x = 1, 故故e-j2

32、 (u+M)x/M+vy/N) = e-j2 (ux/M+vy/N)同理同理, 反变换也具有周期性反变换也具有周期性f(x, y) = f(x+M, y) = f(x, y+N) = f(x+M, y+N)除此之外除此之外, 还有共轭对称性还有共轭对称性F(u, v) = F*(-u, -v)共轭的定义为共轭的定义为: 若一个复数若一个复数Z = x + i y, 则其共轭为则其共轭为Z* = x - i y, 且且 |Z| = |Z*|.37由傅里叶变换的共轭对称性质可以得到结论由傅里叶变换的共轭对称性质可以得到结论: 频谱是关于原点对频谱是关于原点对称的称的.|F(u, v)| = |F(

33、-u, -v)|上述这些性质说明在这里给出的离散傅里叶变换是周期函数上述这些性质说明在这里给出的离散傅里叶变换是周期函数, 并并关于原点对称关于原点对称(经过平移即为中心对称经过平移即为中心对称). 38可分性可分性10/21010/2/2),(1 ),(1),(MxMuxjMxNyNvyjMuxjevxFMNeyxfeMNvuF还可交换顺序还可交换顺序.因而在实际计算中因而在实际计算中, 可逐列或逐行先计算一维的傅里可逐列或逐行先计算一维的傅里叶变换叶变换, 然后对另外一个方向实行计算然后对另外一个方向实行计算 反变换的计算类似。反变换的计算类似。39用前向变换计算傅里叶反变换用前向变换计算

34、傅里叶反变换上式两边取复共轭并除以上式两边取复共轭并除以M, 有有上式的形式和前向变换一样上式的形式和前向变换一样, 左边再取复共轭并乘以左边再取复共轭并乘以M就是对就是对F(u)反反变换的函数变换的函数. )(1)(10/2MxMuxjexfMuF )()(0/2MxMuxjeuFxf )(*1)(*10/2MxMuxjeuFMxfM主要思路是基于傅里叶变换的共轭性质主要思路是基于傅里叶变换的共轭性质, 利用正变换公式计算反变利用正变换公式计算反变换换. 傅里叶变换傅里叶变换40我们知道我们知道, 傅里叶变换可以简化一些运算傅里叶变换可以简化一些运算, 例如例如: 卷积、微分等卷积、微分等.

35、 但由但由于离散傅里叶变换于离散傅里叶变换(反变换反变换)自动将原函数做了周期化的处理自动将原函数做了周期化的处理, 因此在利用因此在利用傅里叶变换处理函数傅里叶变换处理函数(或图象或图象)时时, 一定要注意这种周期化对运算的影响一定要注意这种周期化对运算的影响. 以卷积为例以卷积为例(特别注意在频谱域滤波就相当于在空间域的卷积特别注意在频谱域滤波就相当于在空间域的卷积).以一维为例以一维为例: 如果如果f(x)和函数和函数h(x)各有各有400个点个点(如图如图4.36a和和4.36b), 那么那么, 10)()()()()(Mmmxhmfxhxfxg(4.3.20) 关于周期性的更多讨论关

36、于周期性的更多讨论:l 直接按公式直接按公式(4.6.20)计算计算(注意在定义域外的注意在定义域外的点不参与运算点不参与运算), 它们的卷积函数它们的卷积函数g(x)应该是应该是有有800个点个点, 从从0到到799, 结果函数如图结果函数如图4.36e所示所示(这个是正确的解这个是正确的解).41l 如果把函数如果把函数f(x)和和h(x)先作周期化处理先作周期化处理(为为什么什么?)l 然后利用公式然后利用公式(4.3.20)计算计算, 这时它们的卷这时它们的卷积函数积函数g(x)只有只有400个点个点, 并且在开头有数并且在开头有数据误差据误差, 结尾处出现数据丢失结尾处出现数据丢失4

37、2l 利用离散傅里叶变换的性质利用离散傅里叶变换的性质, 两个函数的卷积可以由它们傅里叶变换两个函数的卷积可以由它们傅里叶变换的乘积做反变换得到的乘积做反变换得到由于傅里叶变由于傅里叶变换换(反变换反变换)附加的附加的周期性周期性(周期的长度周期的长度等于函数的长度等于函数的长度), 如果不谨慎处理周如果不谨慎处理周期问题期问题, 会得到和会得到和上面类似的结果上面类似的结果.43如何处理周期问题如何处理周期问题? 处理的方法是对要处理的函数作零延拓处理的方法是对要处理的函数作零延拓(Padding).假设函数假设函数f(x)和和h(x)分别有分别有A个点和个点和B个点组成个点组成, 对两个函

38、数同时添加零对两个函数同时添加零, 使它们具有相同的周使它们具有相同的周期期, 用用P表示表示PxAAxxffe 010 )(PxBBxxhhe 010 )(只要只要P A+B-1, 将将fe(x)和和he(x)作周期化处理作周期化处理, 然后再用然后再用(4.6.20)作卷积作卷积, 就得到图就得到图4.37e的结果的结果, 这个结果在这个结果在0-799的范围内的范围内, 和原来的卷积结果和原来的卷积结果(图图4.36e)是相同的是相同的. 另一方面另一方面, 在作了适当的函数延拓后在作了适当的函数延拓后, 我们就可以利用傅我们就可以利用傅里叶变换和反变换来计算相应的卷积里叶变换和反变换来

39、计算相应的卷积(或者相应的滤波或者相应的滤波)了了. 4445在频率域计算卷积的正确步骤在频率域计算卷积的正确步骤(1) 对两个函数作零延拓至适当的长度对两个函数作零延拓至适当的长度;(2) 做两个序列的傅里叶变换做两个序列的傅里叶变换(长度和延拓后的一致长度和延拓后的一致);(3) 将两个序列的傅里叶变换相乘将两个序列的傅里叶变换相乘;(4) 计算乘积的傅里叶反变换计算乘积的傅里叶反变换.二维图像滤波问题二维图像滤波问题同样的道理同样的道理, 在处理二维图像在处理二维图像, 例如滤波等时例如滤波等时, 也应遵循同样的理念也应遵循同样的理念.假设假设f(x, y)和和h(x, y)是两幅图像是

40、两幅图像(可把可把h(x, y)看成是空间滤波器看成是空间滤波器), 大小分别大小分别为为A B和和C D, 选选: P A+C-1, Q B+D-1.延拓函数延拓函数f(x, y)和和h(x, y).QxBPxAByAxyxfyxfeor 00 and 10 ),(),(QxDPxCDyCxyxhyxheor 00 and 10 ),(),(46在频谱域作卷积的正确步骤和一维类似在频谱域作卷积的正确步骤和一维类似. 图图4.38给出了几种滤波方式给出了几种滤波方式的结果的结果, 这里这里h(x, y)是滤波器是滤波器H(u, v)的傅里叶反变换的傅里叶反变换. 4.38a是图像未作延是图像未

41、作延拓就实行傅里叶变换拓就实行傅里叶变换, 4.38b是正确的函数延拓是正确的函数延拓, 4.38c是延拓后的图像正是延拓后的图像正确滤波的结果确滤波的结果.注意在这个过程中注意在这个过程中用到了频率域滤波函数用到了频率域滤波函数的傅里叶反变换的傅里叶反变换, 对反对反变换作零延拓变换作零延拓, 然后再然后再进行正变换进行正变换.同样要注意同样要注意反变换的实部和虚部反变换的实部和虚部, 不能在傅里叶计算的中不能在傅里叶计算的中间过程中忽略虚部间过程中忽略虚部, 而而且延拓要同时包括实部且延拓要同时包括实部和虚部和虚部.47在空间域延拓低通滤波器在空间域延拓低通滤波器滤波结果滤波结果48众所周

42、知，两个二维函数众所周知，两个二维函数f(x, y)和和h(x, y)的卷积定义为的卷积定义为:1010),(),(),(),(MmNnnymxhnmfyxhyxf它们之间的傅里叶变换有如下关系它们之间的傅里叶变换有如下关系f(x, y)*h(x, y) F(u, v)H(u, v)f(x, y)h(x, y) F(u, v)*H(u, v)两个二维函数两个二维函数f(x, y)和和h(x, y)的相关函数定义为的相关函数定义为:(没有镜像变换没有镜像变换)其中其中: f*表示表示f的复共轭的复共轭. 注意相关函数和卷积有非常类似的结构注意相关函数和卷积有非常类似的结构, 因此因此f(x, y

43、) h(x, y) F*(u, v)H(u, v)f*(x, y)h(x, y) F(u, v) H(u, v)1010),(),(*),(),(),(MmNnnymxhnmfyxhyxfyxg卷积和相关性理论卷积和相关性理论(引入一个新的概念引入一个新的概念, 相关性相关性)49相关的重要应用是匹配相关的重要应用是匹配.假设假设f(x, y)是一幅包含物体是一幅包含物体或者区域的图像或者区域的图像, 如果想知道中如果想知道中是否包含有感兴趣的物体或区是否包含有感兴趣的物体或区域域, 就让就让h(x, y)作为那个区域或作为那个区域或物体物体(模板模板), 然后求两者的相关然后求两者的相关函数

44、函数.如果存在匹配如果存在匹配(即即f中有感中有感兴趣的物体或区域兴趣的物体或区域), 则相关函则相关函数会在数会在h和和f对应的位置达到最对应的位置达到最大值大值.卷积和相关性理论卷积和相关性理论50二维傅里叶变换的一些基本性质二维傅里叶变换的一些基本性质51二维傅里叶变换的一些基本性质二维傅里叶变换的一些基本性质52532021-7-2653快速傅立叶变换快速傅立叶变换(FFT)n离散傅立叶变换已成为数字信号处理的重要工具离散傅立叶变换已成为数字信号处理的重要工具，然而，然而，它的计算量大，运算时间长它的计算量大，运算时间长，在某种程，在某种程度上却限制了它的使用范围。度上却限制了它的使用

45、范围。n快速算法快速算法大大提高了运算速度大大提高了运算速度，在某些应用场合，在某些应用场合已已能作实时处理能作实时处理，并且应用在控制系统中，并且应用在控制系统中n快速傅立叶变换快速傅立叶变换不是一种新的变换不是一种新的变换，它是离散傅，它是离散傅立叶变换的一种算法，它是在分析离散傅立叶变立叶变换的一种算法，它是在分析离散傅立叶变换中的换中的多余运算的基础多余运算的基础上，进而消除这些重复工上，进而消除这些重复工作的思想指导下得到的作的思想指导下得到的542021-7-2654其原理：其原理：对于一个有限长序列对于一个有限长序列f(x)(0 x N-1)，它的傅立叶变换由下它的傅立叶变换由下

46、式表示：式表示：令令NjNjeWeW212傅立叶变换对可写为：傅立叶变换对可写为：(1)(2) 10101NuxuNxxuWuFNxfWxfuF1, 1 , 0/2exp)()(10NxNuxjxfuFNx55将正变换将正变换(1)展开得到：展开得到：矩矩阵阵表表示示为：为：) 1)(1(1 ) 1(0) 1() 1(22120) 1( 11110) 1(00100) 1() 1 ()0() 1() 1() 1 ()0()2() 1() 1 ()0() 1 () 1() 1 ()0()0(NNNNNNNWNfWfWfNFWNfWfWfFWNfWfWfFWNfWfWfF) 1() 1 () 0(

47、) 1() 1 () 0() 1)(1(1 ) 1(0) 1() 1( 11110) 1(00100NfffWWWWWWWWWNFFFNNNNNN从上式可以看出，要得到每一个频率分量，需进行从上式可以看出，要得到每一个频率分量，需进行N次乘法和次乘法和N-1次加次加法运算。要完成整个变换需要法运算。要完成整个变换需要N2次乘法和次乘法和N(N-1)次加法运算。当序列次加法运算。当序列较长时，必然要花费大量的时间较长时，必然要花费大量的时间562021-7-2656观察上面的系数矩阵，发现观察上面的系数矩阵，发现Wmn是以是以N为周期的，即为周期的，即mnhNnlNmWW)(当当N=8时，其周期

48、性如图所示，由于时，其周期性如图所示，由于NjNeWNj2sin2cos2所以当所以当N=8时，可得：时，可得：jWjWWWNNNN434211W0W1W2W3W4W5W6W7)1)(1(1)1(0)1()1(11110)1(00100NNNNNNWWWWWWWWW572021-7-2657可见，离散傅立叶变换中的乘法运算有许多重复内容。可见，离散傅立叶变换中的乘法运算有许多重复内容。1965年年库利库利-图基提出原始的图基提出原始的N点序列依次分解成一系列短序列，然后，点序列依次分解成一系列短序列，然后，求出这些短序列的离散傅立叶变换，以此来减少乘法运算，例求出这些短序列的离散傅立叶变换，以

49、此来减少乘法运算，例如，设：如，设：由此，离散傅立叶变换可写成下面的形式：由此，离散傅立叶变换可写成下面的形式：12/, 1 , 0) 12()(12/, 1 , 0)2()(21NxxfxfNxxfxf12/0)12(12/0)2(12/0212/0110) 12()2()()()()(NxuxNNxuxNNxxuNNxxuNNxxuNWxfWxfWxfWxfWxfuF582021-7-2658因为：因为：22kNkNWW所以：所以：F1(u)和和F2(u)分别是分别是f1(x)和和f2(x)的的N/2点点的傅立叶变换的傅立叶变换 )()() 12()2() 12()2()(2112/02/

50、12/02/12/02/12/02/uFWuFWxfWWxfWWxfWxfuFuNNxxuNuNNxxuNNxuNxuNNxxuN由于由于F1(u)和和F2(u)均是以均是以N/2为周期，所以为周期，所以 )()2/()()2/(2211uFNuFuFNuF这说明当这说明当mN/2时，上式也是成立的，因此下式成立：时，上式也是成立的，因此下式成立：1, 1 , 0)()()(21NuuFWuFuFuN592021-7-2659由上面的分析可见，一个由上面的分析可见，一个N点的离散傅立叶变换可由两个点的离散傅立叶变换可由两个N/2点的傅立叶变换得到。点的傅立叶变换得到。离散傅立叶变换的计算时间主

51、要由乘法决定，分解后所需的乘法次数离散傅立叶变换的计算时间主要由乘法决定，分解后所需的乘法次数大为减少。第一项为大为减少。第一项为(N/2)2次，第二项为次，第二项为(N/2)2N次，总共为次，总共为2(N/2)2N次运算即可完成，而原来需要次运算即可完成，而原来需要N2次运算，可见分解后的次运算，可见分解后的乘法计算次数减少了近一半。乘法计算次数减少了近一半。当当N为为2的整数幂时，则上式中的的整数幂时，则上式中的F1(u)和和F2(u)还可以再分成两个更短还可以再分成两个更短的序列，因此计算时间会更短。的序列，因此计算时间会更短。由此可见，利用由此可见，利用WNm的周期性和分解运算，从而减少乘法运算次数是的周期性和分解运算，从而减少乘法运算次数是实现快速运算的关键实现快速运算的关键602021-7-2660实例实例61快速傅里叶变换快速傅里叶变换一直是重要的研究课题一直是重要的研究课题.先看看直接利用定义公式计算傅里叶变换的