本节书摘来自异步社区《精通Matlab数字图像处理与识别》一书中的第6章，第6.3节，作者 张铮 , 倪红霞 , 苑春苗 , 杨立红，更多章节内容可以访问云栖社区“异步社区”公众号查看
精通Matlab数字图像处理与识别6.2节介绍了离散傅立叶变换（DFT）的原理，但并没有涉及其实现问题，这主要是因为DFT的直接实现效率较低。在工程实践中，我们迫切地需要一种能够快速计算离散傅立叶变换的高效算法，快速傅立叶变换（FFT）便应运而生。本节将给出快速傅立叶变换算法的原理及其实现细节。
之所以提出快速傅立叶变换（FFT）方法，是因为在计算离散域上的傅立叶变换时，对于N点序列，它的DFT变换与反变换对定义为
于是不难发现，计算每个u值对应的F(u)需要N次复数乘法和N-1次复数加法。因此，为了计算长度为N的序列的快速傅里叶变换，共需要执行N2次复数乘法和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 偶(u) 和F 奇(u) 分别表示f(2r)和f(2r+1)的N/2点DFT。
而且，根据DFT序列的周期性特点，还可得到如下式子成立。
这是一个递推公式，它就是FFT蝶形运算的理论依据。该公式表明，一个偶数长度序列的傅立叶变换可以通过它的奇数项和偶数项的傅立叶变换得到，从而可以将输入序列分成两部分分别计算，并按公式相加/相减。而在这个运算过程中，实际上只需要计算
因此，一个8点按时间抽取的FFT算法的第一步骤如图6.7所示。
图6.7所示是根据式（6-39）绘制的，这一算法也可以用图6.8抽象地表示出来。
由于我们讨论的是基-2的FFT算法，N/2 一般应是偶数，因此得到的序列还可以继续分解，分解过程可以一直持续到每个序列只需要2点的DFT。这样只需要如下的运算即可计算这一DFT值，这一运算是FFT的基本运算，称为蝶形运算。
这一基础单元是对初始输入序列进行傅立叶变换操作的第一步，即2点时的FFT。把这个基本的DFT运算和上面的抽象化蝶形运算比较，可以发现它们的基本结构是完全一致的。在蝶形算法中，我们可以只计算一次 W_N^k F 2，而后分别与F 1相加和相减，从而每一次蝶形算法只需1次复数乘法和2次复数加法（从复杂度分析的角度，相减当然也可看作是一次加法）。并且，注意到 W_N^K =1，因此可以进一步简化计算。尤其第一级蝶形运算更是可以完全简化为单纯的复数加减法。
一个8点FFT的完整计算过程如图6.10所示，请思考这个过程与DFT过程的区别，以及这个过程所需的算法复杂度和存储空间问题。稍后我们将讨论这个问题。
用基-2的时间抽取FFT算法比直接计算DFT的效率高得多。在计算长度为N=2L序列的FFT时，在不对复数乘法进行额外优化的情况下，所需运算量分析如下。
对于每一个蝶形运算，我们需要进行1次复数乘法和2次复数加法。而FFT运算的每一级都含有N/2 = 2L-1个蝶形运算单元。因此，完成L级FFT运算共需要的复数乘法次数M cm和复数加法M ca数目分别为
因此，在N或L取值增大时，FFT运算的优势更加明显。例如，当N=210时，C(N)=102.4，即FFT算法的速度是DFT的102.4倍。
此外，从占用的存储空间看，按时间抽取的FFT算法也远比DFT算法节约。一对复数进行完蝶形运算后，就没有必要再次保留输入的复数对。因此，输出对可以和输入对放在相同的存储单元中。所以，只需要和输入序列大小相等的存储单元即可。也就是一种“原位运算”。
但是，经过观察上面的8点FFT运算全过程，我们发现，如果要使用这种“原位”运算，输入序列就必须按照倒序存储。由于f(x)是逐次抽取的，所以必须对原输入码列倒转位序，得到的次序相当于是原序列编号的二进制码位倒置。也即将原序列编号按照二进制表示，并且将二进制的所有位次序颠倒，就得到了在实际的输入序列中应该使用的排序位置。下面同样以8点FFT为例说明码位倒置的方法，如表6.1所示。
按照表6.1中的顺序排列输入数据，就可以方便地进行原位运算，以节约内存空间。
离散反傅立叶变换的形式与离散傅立叶变换很相似,首先比较它们的公式形式。
离散反傅立叶变换IDFT的公式为
因此，我们只需先将F(u)取共轭，就可以直接使用FFT算法计算IFFT了。
6.3.5 N维快速傅立叶变换N维快速傅立叶变换用于对高维信号矩阵执行傅立叶频谱分析操作。其中二维的快速傅立叶变换常常用于数字图像处理。
N维快速傅立叶变换是由一维FFT组合而成的，其运算实质就是在给定二维或多维数组的每个维度上依次执行一维FFT，并且使用“原位运算”的方法。在开始之前，算法将输入直接复制到输出上，所以之后在每个维度上执行FFT的原位操作都不会改变原本的输入数组，同时也使这个算法输出的数组和输入的数组拥有同样的大小和维度。也就是说，如果对一幅灰度图像执行二维快速傅立叶变换操作，得到的结果也将是一个二维数组。
Matlab中提供了fft2和ifft2函数分别计算二维傅立叶变换和反变换，它们都经过了优化，运算速度非常快；另一个与傅立叶变换密切相关的函数是fftshift，常需要利用它来将傅立叶频谱图中的零频点移动到频谱图的中心位置。
下面分别介绍这3个函数。
1．fft2函数该函数用于执行二维快速傅立叶操作，因此可以直接用于数字图像处理。调用语法为
Y = fft2(X) Y = fft2(X,m,n)参数说明
X为输入图像。m和n分别用于将X的第一和第二维规整到指定的长度。当m和n均为2的整数次幂时，算法的执行速度要比m和n均为素数时更快。返回值
Y是计算得到的傅立叶频谱，是一个复数矩阵。提示
计算abs(Y )可得到幅度谱，计算angle(Y )可得到相位谱。2．fftshift函数在fft2函数输出的频谱分析数据中，是按照原始计算所得的顺序来排列频谱的，而没有以零频为中心来排列，因此造成了零频在输频谱矩阵的角上，显示幅度谱图像时表现为4个亮度较高的角（零频处的幅值较高），如图6.11!（a）所示。image
fftshift函数利用了频谱的周期性特点，将输出图像的一半平移到另一端，从而使零频被移动到图像的中间。其调用语法为
Y = fftshift(X) Y = fftshift(X,dim)参数说明
X为要平移的频谱。dim指出了在多维数组的哪个维度上执行平移操作。返回值
Y是经过平移的频谱。利用fftshift函数对图6.11（a）中的图像平移后的效果如图6.11（b）所示。
下面给出对于二维图像矩阵，fftshift函数的平移过程，如图6.12所示。
可见，输出矩阵被分为了4个部分，其中1、3两部分对换，2、4两部分对换。这样，原来在角上的零频率点（原点）位置就移动到了图像的中央位置。而dim参数则可以指定在多维数组的哪个维度上执行对换操作。例如，对于矩阵而言，dim取1和2的情形分别如图6.13所示。
3．ifft2函数该函数用于对图像（矩阵）执行逆傅立叶变换。输出矩阵的大小与输入矩阵相同。调用形式为
Y = ifft2(X) Y = ifft2(X,m,n)参数说明
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)); %对幅值做对数变换以压缩动态范围 subplot(1,2,1); imshow(temp,[]); title('FFT'); subplot(1,2,2); imshow(I1); title('Source') I2 = imread('circuit.tif'); %读入原图像 fcoef = fft2(double(I2)); %做fft变换 spectrum = fftshift(fcoef); %将零点移到中心 temp =log(1+abs(spectrum)); %对幅值做对数变换以压缩动态范围 figure; subplot(1,2,1); imshow(temp,[]); title('FFT'); subplot(1,2,2); imshow(I2); title('Source')上述程序的运行结果如下。
由图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中关于幅度谱和相位谱各自作用的讨论，您能想到这样做将会产生怎样的结果吗?
% ex6_2.m % 读取图片 A = imread('../beauty.jpg'); B = imread('../cat.jpg'); % 求傅立叶变换 Af = fft2(double(A)); Bf = fft2(double(B)); % 分别求幅度谱和相位谱 AfA = abs(Af); AfB = angle(Af); BfA = abs(Bf); BfB = angle(Bf); % 交换相位谱并重建复数矩阵 AfR = AfA .* cos(BfB) + AfA .* sin(BfB) .* i; BfR = BfA .* cos(AfB) + BfA .* sin(AfB) .* i; % 傅立叶反变换 AR = abs(ifft2(AfR)); BR = abs(ifft2(BfR)); % 显示图像 subplot(2,2,1); imshow(A); title('美女原图像'); subplot(2,2,2); imshow(B); title('猫的原图像'); subplot(2,2,3); imshow(AR, []); title('美女的幅度谱和猫的相位谱组合'); subplot(2,2,4); imshow(BR, []); title('猫的幅度谱和美女的相位谱组合');程序运行结果如图6.15所示。
通过这个示例，我们可以发现，交换相位谱之后，反变换之后得到的图像内容与其相位谱对应的图像一致，这就印证了我们之前关于相位谱决定图像结构的论断。而图像中整体灰度分布的特性，如明暗、灰度变化趋势等则在比较大的程度上取决于对应的幅度谱，因为幅度谱反映了图像整体上各个方向的频率分量的相对强度。
相关资源：精通Matlab数字图像处理与识别 (张铮等) 扫描版pdf