第25章 DSP变换运算-快速傅里叶变换原理（FFT）
在数字信号处理中常常需要用到离散傅立叶变换(DFT)，以获取信号的频域特征。尽管传统的DFT算法能够获取信号频域特征，但是算法计算量大，耗时长，不利于计算机实时对信号进行处理。因此导致DFT被发现以来，在很长的一段时间内都不能被应用到实际工程项目中，直到一种快速的离散傅立叶计算方法——FFT被发现，离散是傅立叶变换才在实际的工程中得到广泛应用。需要强调的是，FFT并不是一种新的频域特征获取方式，而是DFT的一种快速实现算法。
特别声明：FFT原理的讲解来自网络和书籍。
25.1 FFT由来
25.3 直接计算DFT的问题及改进路径
25.3 改善DFT运算效率的基本途径
25.4 按时间抽选的基2-FFT算法
25.5 按频率抽选的基2-FFT算法
25.6 总结
25.1 初学者重要提示
1、 为大家推荐西电的国家级精品课程信号与系统，含课堂视频，书籍和课堂PPT
http://www.armbbs.cn/forum.php?mod=viewthread&tid=94886 。
2、 非常好的DSP基础知识普及书籍：
http://www.armbbs.cn/forum.php?mod=viewthread&tid=97312 。
25.2 FFT的由来
离散傅里叶变换（Discrete Fourier Transform，DFT）是数字信号处理最重要的基石之一，也是对信号进行分析处理时，最常用的工具之一。在200多年前，法国数学家、物理学家傅里叶提出以他的名字命名的傅里叶级数之后，用DFT这个工具来分析信号就已经被人们所知。
历史上最伟大的数学家之一，欧拉是第一个使用“函数”一词来描述包含各种参数的表达式，例如：y = f(x)。他是把微积分应用于物理学的先驱者之一，给出了一个用实变量函数表示傅立叶系数的方程，用三角级数来描述离散声音在媒介中传播，发现某些函数可以通过余弦函数之和来表达。 但在很长时间内，这种分析方法并没有引起更多的重视，最主要的原因在于这种方法运算量比较大。直到1965年，Cooley和Tukey在《计算机科学 》发表著名的《机器计算傅立叶级数的一种算法》论文，FFT才开始大规模应用。
那个年代，有个肯尼迪总统科学咨询委员会，其中有项研究主题是对苏联核测试进行检测，Tukey就是其中一员。美国/苏联核测试提案的批准，主要取决于不实地访问核测试设施而做出检测方法。其中一个想法是，分析离海岸的地震情况，这种计算需要快速算法来计算DFT。其它应用是国家安全，如用声学探测远距离的核潜艇。所以在军事上，迫切需要一种快速的傅立叶变换算法，这也促进了FFT的正式提出。
FFT充分利用了DFT运算中的对称性和周期性，从而将DFT运算量从N2减少到 。当N比较小时，FFT优势并不明显。但当N大于32开始，点数越大，FFT对运算量的改善越明显。比如当N为1024时，FFT的运算效率比DFT提高了100倍。在库利和图基提出的FFT算法中，其基本原理是先将一个N点时域序列的DFT分解为N个1点序列的DFT，然后将这样计算出来的N个1点序列DFT的结果进行组合，得到最初的N点时域序列的DFT值。实际上，这种基本的思想很早就由德国伟大的数学家高斯提出过，在某种情况下，天文学计算（也是现在FFT应用的领域之一）与等距观察的有限集中的行星轨道的内插值有关。由于当时计算都是靠手工，所以产生一种快速算法的迫切需要。 而且，更少的计算量同时也代表着错误的机会更少，正确性更高。高斯发现，一个富氏级数有宽度N=N1*N2，可以分成几个部分。计算N2子样本DFT的N1长度和N1子样本DFT的N2长度。只是由于当时尚欠东风——计算机还没发明。在20世纪60年代，伴随着计算机的发展和成熟，库利和图基的成果掀起了数字信号处理的革命，因而FFT发明者的桂冠才落在他们头上。
之后，桑德（G.Sand）-图基等快速算法相继出现，几经改进，很快形成了一套高效运算方法，这就是现在的快速傅立叶变换（FFT）。这种算法使DFT的运算效率提高1到2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了良好的条件，大大推进了数学信号处理技术。1984年，法国的杜哈梅（P.Dohamel）和霍尔曼（H.Hollamann）提出的分裂基块快速算法，使运算效率进一步提高。
库利和图基的FFT算法的最基本运算为蝶形运算，每个蝶形运算包括两个输入点，因而也称为基-2算法。在这之后，又有一些新的算法，进一步提高了FFT的运算效率，比如基-4算法，分裂基算法等。这些新算法对FFT运算效率的提高一般在50%以内，远远不如FFT对DFT运算的提高幅度。从这个意义上说，FFT算法是里程碑式的。可以说，正是计算机技术的发展和FFT的出现，才使得数字信号处理迎来了一个崭新的时代。除了运算效率的大幅度提高外，FFT还大大降低了DFT运算带来的累计量化误差，这点常为人们所忽略。
25.3 直接计算DFT的问题及改进路径
25.3.1 问题的提出
设有限长序列x(n)，非零值长度为N，若对x(n)进行一次DFT运行，共需要多大的运算工作量。
25.3.2 DFT的运算量
DFT和IDFT的变换式：
下面以DFT为例说明计算量：
计算机运算时（编程实现）：
由上面的结算可得DFT的计算量如下：
复数乘法的计算量：(a+jb)(c+jd)=(ab-bd)+j(bc+ad)
下面通过两个实例来说明计算量：
例一：计算一个N点DFT，共需 次复乘。以做一次复乘1 计算，若N=4096，所需时间为
例二：石油勘探，有24个通道的记录，每通道波形记录长度为5秒，若每秒抽样500点/秒。
- 每道总抽样点数：500*5 = 2500点。
- 24道总抽样点数：24*2500=6万点。
由于计算量大，且要求相当大的内存，难以实现实时处理，限制了DFT的应用，人们一直在寻求一种能提高DFT运算速度的方法。
FFT便是Cooley和Tukey在1995年提出来的快速算法，它可以使运算速度提高几百倍，从而使数字信号处理成为一个新兴的应用学科。
25.4 改善DFT运算效率的基本途径
1、利用DFT运算的系数 的固有对称性和周期性，改善DFT的运算效率。
1）对称性
2）周期性
3）可约性
2、将长序列DFT利用对称性和周期性分解为短序列DFT的思路
因为DFT的运算量与N2成正比，如果一个大点数N的DFT能分解为若干小点数DFT的组合，则显然可以达到减少运算工作量的效果。
FFT算法的基本思想：
- 利用DFT系数的特性，合并DFT运算中的某些项。
- 把长序列DFTà短序列DFT，从而减少运算量。
FFT算法分类：
时间抽选法
DIT: Decimation-In-Time
频率抽选法
DIF: Decimation-In-Frequency
25.5 按时间抽选的基2-FFT算法
25.5.1 算法原理
设输入序列长度为N = 2M（M为正整数），将该序列按时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的FFT算法。也称为Coolkey-Tukey算法。
其中基2表示：N = 2M，M为整数。若不满足这个条件，可以人为地加上若干零值（加零补长）使其达到N = 2,M。
25.5.2 算法步骤
- 分组，变量置换
- 分组，变量置换
其中k = 0， 1，…. N/2 – 1。 和 只有N/2个点，以N/2为周期；而X(k)却有N个点，以N为周期。要用x1(k)和x2(k)表达全部的X(k)值。还必须利用
系数的周期特性。
有了上面的计算结果后，我们可以得到如下的蝶形运算流图符号：
关于这个蝶形运算流图符号说明如下：
- 1个蝶形运算需要1次复乘，2次复加。
- 左边两路为输入。
- 右边两路为输出。
- 中间以一个小圆表示加减运算（右上路为相加输出，右下路为相减输出）。
- 分解后的运算量
运算量减少了近一半。
例子：求N=23=8点FFT变化。按N=8àN/2=4，做4点的DFT，先将N=8点的DFT分解成2个4点的DFT:
可知： 时域上 x(0)，x(2)，x(4)，x(6)为偶子序列。
x(1)，x(3)，x(5)，x(7)为奇子序列。
频域上 X(0) 到X(3) 由X(k)给出
X(4) 到X(7) 由X(k+N/2)给出
此外，还有4个蝶形结，每个蝶形结需要1次复乘，2次复加。一共是：复乘4次，复加8次。
用分解的方法的得到X(k)需要：
复乘：32+4=36次
复加：24+8=32次
N = 23 = 8 按时间抽取的DFT分解过程：
因为4点DFT还是比较麻烦，所以再继续分解。
若将N/2(4点)子序列按奇/偶分解成两个N/4点(2点)子序列。即对将x1(r)和x2(r)分解成奇、偶两个N/4点(2点)的子序列。
因此可以对两个N/2点的DFT再分别作进一步的分解。将一个８点的DFT可以分解成四个２点的DFT,直到最后得到两两点的DFT为止。
由于这种方法每一步分解都是按输入序列是属于偶数还是奇数来抽取的，所以称为“按时间抽取的FFT算法”。
下图是由4个两点DFT组成的8点DFT：
下图是按8点抽取的FFT运算流图：
这里注意观察蝶形图的系数
25.5.3 FFT算法和直接计算DFT运算量的比较
FFT算法与直接计算DFT所需乘法次数的比较曲线
25.6 按频率抽选的基2-FFT算法
在基2快速算法中，频域抽取法FFT也是一种常用的快速算法，简称DIF-FFT。
鉴于网上和课本中关于FFT原理已经讲解非常详细了，在这里就不再赘述了。有兴趣的查阅相关书籍进行学习即可。
25.7 总结
本章节主要讲解了FFT的基2算法实现原理，讲解稍显枯燥，不过还是希望初学的同学认真学习，搞懂一种快速傅里叶算法的实现即可。
微信公众号：armfly_com
