《数字信号处理》课件1第4章

第4章快速傅里叶变换(FFT)4.1直接计算DFT的运算量及改进途径

4.2时间抽取法(DIT)基-2FFT算法

4.3频域抽取法(DIF)基-2FFT算法

4.4快速傅里叶逆变换(IFFT)算法

4.5

FFT算法的工程实现考虑

习题

4.1直接计算DFT的运算量及改进途径

4.1.1直接计算DFT的运算量

设x(n)是长度为N的有限长序列，其N点DFT为

离散傅里叶逆变换(IDFT)为(4.1.2)(4.1.1)上述两式的差别只在于WN上指数符号的不同，以及一个常数因子，因此，DFT和IDFT运算量基本上是相同的，为

简便，下面只讨论DFT的运算量。

由于实数序列是复数序列的特例，为分析简便，这里统一考虑x(n)为复数序列的情况，而通常也是复数，因此，这里以复数乘法和复数加法的次数作为运算量衡量指标。显然，计算一个X(k)值需要N次复数乘法，N－1次复数加法。而X(k)总共有N个值，计算所有X(k)值的运算量为复数乘法N2次，复数加法N(N－1)次。对于复数乘法和加法，实际上是转化为实数进行运算的，式(4.1.1)可以写成(4.1.3)可见，1次复数乘法相当于4次实数乘法和2次实数加法，1次复数加法相当于2次实数加法。因此，计算一个X(k)需要4N次实数乘法和2N+2(N－1)=4N－2次复数加法，整个DFT运算需要4N2次实数乘法和N(4N－2)次实数加法。

当N

1时，N(N－1)≈N2，N(4N－2)≈4N2，因此，不管以复数还是实数进行统计，直接计算N点DFT的乘法或加法次数都与N2成正比，随着N的增加，运算量增加越来越快，特别是N很大时，运算量将非常可观。例如，N=8时，次数为N2=64；N=1024时，次数为N2=1048576，超过100万次。对于各种实时性很强的信号处理应用来说，要求计算速度特别快，因此必须改进DFT的计算方法，减少运算量。＞＞4.1.2减少运算量的途径

由于N点DFT的运算量随N2快速增长，当N增加1倍时，N2增加到4倍。如果能够将N点DFT分解为几个较短的DFT，运算量将会大大减少。例如，分解为2个N/2点DFT，复数乘法次数为 ，运算量减少1倍；若分解为4个点DFT，复数乘法次数为 ，运算量将减少为原来的。可以说，正是DFT分解思想形成了快速计算DFT的基本途径。以N点DFT分解为2个N/2点DFT为例，假定N点序列x(n)，n=0，1，2，…，N－1，N为偶数，那么将x(n)分解为2个N/2点序列的方法归纳起来有两个：奇偶分解和前后分解。

(1)奇偶分解：x(n)分解为偶序列和奇序列，分别表示为

偶序列：x(2r)，r=0，1，2，…，－1

奇序列：x(2r+1)，r=0，1，2，…，－1

(2)前后分解：x(n)前后对半分解为两部分，分别为

前半部分：x(r)，r=0，1，2，…，－1

后半部分： ，r=0，1，2，…，－1快速傅里叶变换(FFT)通过不断地将长序列的DFT分解为短序列的DFT，并利用的特性，来达到减少DFT运算量的目的。为了便于对N点DFT进行分解，通常N取为2的幂级数形式，即N=2M，M为整数。此时的快速傅里叶变换称为基-2FFT算法。若序列长度不满足条件N=2M，可以对序列补0，使之达到这一条件。

FFT算法推导过程中，主要利用的对称性和周期性，在序列分解后，对DFT计算式(4.1.1)中的某些项进行合并，从而减小乘法和加法的次数。称为旋转因子，其对称性和周期性表现为

对称性：

周期性：(4.1.4)(4.1.5)(4.1.6)此外，的特性还有

下面两节中将学习两类基-2FFT算法，分别由上述两种分解方式推导得到，其中，奇偶分解对应于时间抽取法FFT(Decimation-In-TimeFFT)，简称DIT-FFT，前后分解对应于频率抽取法FFT(Decimation-In-FrequencyFFT)，简称DIF-FFT。(4.1.7)

4.2时间抽取法(DIT)基-2FFT算法

4.2.1

DIT-FFT算法原理

设序列x(n)长度N=2M，N点DIT-FFT算法对应序列奇偶分解，令x1(r)、x2(r)分别表示偶序列和奇序列，则有

(4.2.2)(4.2.1)根据DFT的定义式(4.1.1)，有(4.2.3)由于

所以(4.2.4)观察上式：右边前一项是序列x1(r)的N/2点DFT，后一项求和部分是序列x2(r)的N/2点DFT，即

对于N/2点DFTX1(k)和X2(k)，变量k的取值只有X(k)中k的取值的一半。因此，对于X(k)表达式(4.2.4)，(4.2.6)(4.2.5)

(1) ：确定X(k)前半部分。

(2) ：确定X(k)的后半部分。(4.2.7)为表述方便，X(k)后半部分表示为 。由于点DFTX1(k)和X2(k)具有周期性，且周期均为N/2，即

而 ，因此，X(k)的后半部分为(4.2.8)由此可见：N点DFT可以分解为两个N/2点DFT，按照式(4.2.7)和式(4.2.8)又可组合成N点DFT。因此求X(k)时，只要求出k=0，1，…，N/2－1时的X1(k)和X2(k)值，即可得到所有的X(k)值，即k=0，1，…，N－1，从而大大节省了运算量，这也是快速傅里叶变换的特点和好处所在。

式(4.2.7)和式(4.2.8)的运算可以用图4.2.1所示的蝶形运算流图符号表示。图中，左侧为两个输入节点，右侧为两个输出节点，左下支路上标注系数，没有标注时系数默认为1，右上支路默认为相加运算，右下支路为相减运算。图4.2.1

DIT-FFT的蝶形运算流图符号从图4.2.1中可以看出，一个蝶形运算需要1次复数乘法和2次复数加法。

通过采用蝶形运算流图符号，DIT-FFT经过1次DFT分解之后，整个分解过程可用图4.2.2所示的运算流图表示，其中DFT点数N=23=8，蝶形输出值X(0)～X(3)由式(4.2.7)确定，X(4)～X(7)由式(4.2.8)确定。由于1个蝶形包括两个输入和两个输出，总共有N/2个蝶形。图4.2.2

DIT-FFT的一次分解运算流图(N=8)从图4.2.2可以看出，一个N点DFT可以由分解的两个N/2点DFT通过N/2个蝶形进行合成得到，总运算量包括两个N/2点DFT以及N/2个蝶形的计算。N/2个蝶形运算需要N/2次复数乘法和2×N/2=N次复数加法。

对于两个N/2点DFT，如果直接计算，需要复数乘法次数为复数加法次数为

因此，经过一次分解后的N点DFT运算量为

复数乘法：

复数加法：与直接计算N点DFT的运算量相比，一次DFT分解后的运算量减少了一半左右。这也充分说明：通过DFT分解可以有效地减少DFT的计算量。如果针对N/2点DFT也采用分解措施，将一个N/2点DFT分解为两个N/4点DFT，那么可以进一步减少运算量。这就是下面讨论的DFT的二次分解过程。

不妨以N/2点DFTX1(k)为例，将N/2点序列x1(r)按照奇偶分解方式进行分解，得到两个N/4点序列，分别为(4.2.10)(4.2.9)则有(4.2.11)上式右边前一部分、后一部分求和项分别是序列x3(l)和x4(l)的N/4点DFT，即(4.2.12)(4.2.13)因此，当 时，式(4.2.11)可直接表述为

当

时，利用X3(k)、X4(k)周期性

以及

性质，X1(k)表示为(4.2.14)(4.2.15)图4.2.3给出了N=8时，一个N/2点DFT分解为两个N/4点DFT的蝶形运算流图，即X1(k)分解为X3(k)和X4(k)，同时由X3(k)和X4(k)通过N/4个蝶形也可以合成X1(k)。图4.2.3

N/2点DFT分解的蝶形运算流图(N=8)同理，对于N/2点DFTX2(k)也可以分解为两个N/4点DFT：(4.2.16)(4.2.17)其中，X5(k)、X6(k)分别为x2(r)分解的偶序列和奇序列对应的N/4点DFT。(4.2.18)(4.2.19)图4.2.4给出了N=8点DFT经过两次分解后的蝶形运算流图。与第一次分解后的运算量相比，利用4个N/4点DFT及两级蝶形来计算N点DFT的运算量又降低了。图4.2.4

DIT-FFT的二次分解运算流图(N=8)可以看出，当N=8时，通过两次DFT分解后，N/4点DFT实际上为2点DFT，即X3(k)、X4(k)、X5(k)、X6(k)，k=0，1。此时，2点DFT可以直接进行计算。以X3(k)计算式(4.2.12)为例，可得

而 ， ，因此有(4.2.20)(4.2.21)(4.2.22)式(4.2.21)和式(4.2.22)表明：2点DFT仅涉及加减法运算，不需要乘法运算；X4(k)、X5(k)、X6(k)也具有类似特点，它们都可用一个简单的2点蝶形运算表示。

依此类推，一个N=2M点DFT通过M－1次分解后，可以分解为N/2个2点DFT，得到M级蝶形运算。对于8点DFT，通过二次分解后，可以得到三次蝶形运算，图4.2.5给出了完整的8点DIT-FFT蝶形运算流图。图4.2.5

DIT-FFT的蝶形运算流图(N=8)图4.2.5中，旋转因子采用的表示形式；输出为顺序排列，但输入并不是顺序排列，而是在每一次分解过程中，将输入序列按照时间上的偶数和奇数次序分解为两个短序列，相当于在时间上进行抽取。最后得到的输入序列也是非常有规律的，将在后面进行介绍。所以，具有这种时间抽取关系的快速傅里叶变换称为“时间抽取法FFT(DIT-FFT)”。4.2.2

DIT-FFT运算量分析与比较

根据时间抽取法FFT算法的蝶形运算流图可知，当N=2M时，共有M级蝶形运算，每级均有N/2个蝶形，而每个蝶形运算包含1次复数乘法和2次复数加法。因此，每一级蝶形都需要N/2次复数乘法和N次复数加法。M级蝶形的总运算量为

复数乘法：

复数加法：(4.2.23)(4.2.24)

由于旋转因子存在一些特例，如

，

，

，

，与这几个系数相乘实际上不需要乘法运算，这种情况在直接计算DFT时也存在。但是当N较大时，这种特例相对较少。为了便于统一比较运算量，这里不考虑这些特殊情况。

表4.2.1列出了N点DFT直接计算和DIT-FFT计算的运算量，两者复数乘法和复数加法之比分别为

复数乘法之比： (4.2.25)复数加法之比：(4.2.26)表4.2.1

N点DFT直接计算和DIT-FFT的运算量比较当N

1时，N

log2N，有N2

N·log2N，N(N－1)

N·log2N，因此，与直接计算DFT相比，DIT-FFT的运算量大大减少。表4.2.2列出了不同N值条件下直接计算DFT与DIT-FFT的复数乘法次数及比例关系。可以看出，随着N增大，复数乘法次数的比值越大，DIT-FFT的优势越来越明显。

但是也要注意到：当N较小时(如N≤16)，比值也相对较小，考虑到DIT-FFT实际编程时的复杂性和指令开销，DIT-FFT的整体运算量不一定小于直接计算DFT。因此，在实际计算DFT时，需要根据N的大小，在直接计算和FFT之间灵活选择。＞＞＞＞＞＞＞＞表4.2.2直接计算DFT与DIT-FFT的复数乘法次数的比较4.2.3

DIT-FFT运算规律

为了更好地理解和掌握时间抽取法FFT算法，为算法的实际编程和硬件实现打下良好的基础，下面对DIT-FFT的运算规律和特点进行分析和讨论。

1.原位运算

从图4.2.5所示的DIT-FFT蝶形运算流图中可以看出：N=2M点FFT共有M级蝶形运算，每级由N/2个蝶形构成。在同一级中，每个蝶形的输入和输出都位于同一水平线上，并且每个输入只参与本蝶形运算，与其它蝶形无关。该特性意味着蝶形的输出可以直接存入输入所占用的存储单元，这就是原位计算的特点。通过原位计算，每一级的N/2蝶形运算完成后，所有输出存入原输入的存储位置，然后开始下一级的蝶形运算，只不过蝶形运算的组合关系有所不同。这种原位计算结构只需要N个存储单元，节省了存储开销，降低了设备成本。

2.位码倒序

观察图4.2.5所示的蝶形运算流图的输入输出可以看出，输出序列是按照X(0)，X(1)，…，X(7)的顺序排列，而输入序列次序是x(0)，x(4)，…，x(7)，看起来似乎很乱，但实际上是有规律的，这种规律称为位码倒序。首先看看输入序列是如何形成x(0)，x(4)，…，x(7)排列的。造成这种排列关系的原因是序列x(n)不断地按照n的偶奇特性进行分解。假设n用二进制数表示为(n2n1n0)，那么，第一次分解是按照n0=0和n0=1分解为偶数序列和奇数序列，第二次分解是分别针对偶数序列和奇数序列，按照n1=0和n1=1进行偶奇分解，最后得到的2点序列是按照n2=0和n2=1排列的。这种不断分解为偶数序列和奇数序列的过程可用图

4.2.6表示。图4.2.6形成位码倒序的树状图若DIT-FFT输入顺序编号0，1，2，…，7用二进制码(n2n1n0)表示，那么图4.2.6中DIT-FFT输入序列可表示为x(n0n1n2)，其序号(n0n1n2)实际上是二进制码(n2n1n0)的比特左右反转结果，两者形成倒序关系。因此称为位码倒序。

表4.2.3列出了N=8时顺序二进制数及对应的倒序二进制数。给定顺序的输入序列x(n)，计算DIT-FFT时，将位码倒序十进制数作为序号来选择x(n)。表4.2.3顺序二进制数与倒序二进制数的对照表在实际编程实现DIT-FFT算法的过程中，可以采取以下方式：

(1)若采用通用计算机编程，可按照表4.2.3，依次将十进制顺序转换为二进制数、位码倒序二进制数以及最后的位码倒序十进制数，并依据位码倒序十进制数来选择x(n)作为DIT-FFT的输入。

(2)若通用计算机的存储空间足够，可将顺序与位码倒序十进制数的对应关系以表的形式存储起来，通过查表方式来选择x(n)。此时，只需要一个数组存储位码倒序十进制数，数组位置表示顺序号，数组的值代表位码倒序十进制数。

(3)若采用数字信号处理器(DSP)编程，可利用DSP自身的位码倒序寻址专用指令来完成转换。以美国TI公司的TMS320C54系列DSP为例，假定N=8，辅助寄存器AR2指向x(0)的存储单元，辅助寄存器AR0设置为FFT点数的一半，即AR0=4，那么，位码倒序寻址的专用指令为

*AR2+0B

该指令表示用反向进位的方式将AR0加至AR2上，即加法按比特从高位向低位进位，然后再赋值给AR2，*表示AR2指向地址的数值。注意到AR2初始地址低3位必须为零，以便进行反向进位。以低3位运算为例，初始值为0，以反向进位方式依次加4，可以得到4、2、6、1、5、3、7。运算完毕后，AR0=4始终固定不变，AR2则按照位码倒序的方式依次指向x(0)，x(4)，…，x(7)。

3.蝶形运算规律

从图4.2.5中N=8点DIT-FFT蝶形运算流图可以看出，从左至右，第一级蝶形对应2点FFT，输入数据相距1点，或者说蝶形张口大小为1；第二级蝶形输入数据相距2点，蝶形张口大小为2；第三级蝶形输入数据相距4点，蝶形张口大小为4。依此类推，对于N=2M点DIT-FFT，从左至右第m级蝶形输入数据相距2m－1点，蝶形张口大小也为2m－1。利用蝶形张口大小的特点，便于从前一级输出中选择相应数据作为输入，来进行本级的蝶形运算。与蝶形运算密切相关的有旋转因子，观察图4.2.5可知，旋转因子的个数与蝶形级数有关，第m级蝶形的旋转因子有2m－1个，可以表示为

对于最后一级蝶形，m=M，旋转因子有2M－1=N/2个。由于

因此，最后一级蝶形的旋转因子均包含着前面M－1级蝶形的旋转因子，所有旋转因子以集合形式表示为

。(4.2.28)(4.2.27)表4.2.4给出了不同蝶形级数下的蝶形张口大小和旋转因子。在实际DSP编程实现DIT-FFT时，可以将旋转因子 制作成表的形式，然后根据式(4.2.27)中蝶形级数与WN指数k·2M－m的关系，查表得到当前蝶形运算所需要的旋转因子。这样可以避免直接计算复指数，能够减少FFT运算量。表4.2.4蝶形张口大小和旋转因子与蝶形级数的关系4.2.4

DIT-FFT其它形式流图

对于时间抽取法FFT算法，图4.2.5的蝶形运算流图并不是唯一的，只要能够保持各节点间的支路及其传输系数不变，不论如何改变输入节点、输出节点以及中间节点的排列顺序，所得到的运算流图是等效的。这样，通过对图4.2.5进行变形，就可以得到DIT-FFT其它形式的运算流图。

图4.2.5中，蝶形运算流图的输入为倒序，输出为顺序。通过变形，图4.2.7给出了输入为顺序、输出为倒序的运算流图形式，该流图同样具有原位计算的特点，其旋转因子、运算量也与图4.2.5相同，只是在蝶形张口大小次序和旋转因子排列上有所差别。若从左至右来看蝶形张口，图4.2.5中是由小变大，而图4.2.7中是由大变小。图4.2.7

DIT-FFT的变形运算流图(输入顺序、输出倒序)对于旋转因子，图4.2.5中最后一级按照 顺序排列，而图4.2.7中最后一级是按照 排列的，并且前一级的旋转因子正好是本级上一半蝶形运算的旋转因子，顺序也不变。这种流图形式就是最初由库利和图基给出的时间抽取法FFT。

如果要获得输入和输出均是顺序排列的运算流图，可以对图4.2.7的最后一级蝶形输出进行调整，得到图4.2.8所示的运算流图。该流图的旋转因子、运算量均与图4.2.7相同，但是在最后一级上不能采用原位计算。图4.2.8

DIT-FFT的变形运算流图(输入顺序、输出顺序) 4.3频域抽取法(DIF)基-2FFT算法

4.3.1

DIF-FFT算法原理

设序列x(n)长度N=2M，N点DIF-FFT算法对应着x(n)前后对半分解为两部分，即前半部分x(n)和后半部分 。根据DFT的定义有(4.3.1)

由于 ，故

需要说明的是，上式中旋转因子项是，而不是，因此，上式并不是一个N/2点DFT。式(4.3.2)要根据k为偶数或者奇数，将X(k)分为两部分进行讨论。(4.3.2)

(1)k为偶数。令k=2r，

，则有(4.3.3)

(1)k为奇数。令k=2r+1，

，则有(4.3.4)可以看出，式(4.3.3)和式(4.3.4)都是N/2点DFT表达式，其中，式(4.3.3)变换对象是x(n)前半部分和后半部分相加形成的序列，式(4.3.4)变换对象则是x(n)前半部分和后半部分相减后再乘以形成的序列。定义两个序列为(4.3.5)(4.3.6)将上述两式分别代入式(4.3.3)和式(4.3.4)，可得

式(4.3.7)和式(4.3.8)表明：N点DFT按照k的奇偶特性，可以分解为两个N/2点DFT。具体方法是将x(n)前后对半分解为两部分，合成两个新的N/2点序列，再进行N/2点DFT。合成序列x1(n)、x2(n)与x(n)的关系可用图4.3.1所示的蝶形运算流图符号表示。(4.3.7)(4.3.8)

利用上述蝶形运算流图符号，N=8点DFT经过一次分解后，得到的运算流图如图4.3.2所示。图4.3.1

DIF-FFT的蝶形运算流图符号图4.3.2

DIF-FFT一次分解的运算流图(N=8)与时间抽取法的推导过程一样，由于N=2M，N/2仍为偶数，在图4.3.2的基础上，可以将每个N/2点DFT进一步分解为两个N/4点DFT。这相当于分别将x1(n)和x2(n)进行前后对半分解后，通过蝶形运算，合成为4个N/4点序列，再进行DFT。图4.3.3给出了N=8点DFT经过二次分解后的运算流图。

当N=8时，经过两次分解得到的N/4点DFT即为2点DFT，可以直接进行计算，相当于一个基本的蝶形运算符号。一个N=2M点DFT通过M－1次分解后，最后可分解为N/2个2点DFT，形成M级蝶形运算。图4.3.4给出了完整的8点DIF-FFT蝶形运算流图。图4.3.3

DIF-FFT二次分解的运算流图(N=8)图4.3.4

DIF-FFT的蝶形运算流图(N=8)观察图4.3.4可知，DIF-FFT的蝶形运算流图仍具有原位计算的特点，其输入序列是顺序的，而输出是倒序的。这是由于每一级分解的输出都按照k的奇偶次序分成为两部分，这相当于在频率上进行抽取，最后得到位码倒序的输出。因此，具有这种频率抽取关系的快速傅里叶变换称为“频率抽取法FFT(DIF-FFT)”。4.3.2

DIF-FFT与DIT-FFT的比较

比较图4.2.5中DIT-FFT蝶形运算流图和图4.3.4中DIF-FFT的蝶形运算流图，两者相同点如下：

(1)原位运算。对于DIT-FFT和DIF-FFT，每个蝶形的输入和输出都位于同一水平线上，并且每个输入只参与本蝶形运算，蝶形的输出可直接存入输入所占用的存储单元。

(2)运算量相同。当N=2M时，DIT-FFT和DIF-FFT都有M级蝶形运算，每级均有个蝶形，复数乘法总数为log2N次，复数加法总数为N·log2N次。

DIT-FFT和DIF-FFT的差异如下：

(1)输入和输出排列次序不同。DIT-FFT输入为倒序、输出为顺序，DIF-FFT输入为顺序、输出为倒序。DIT-FFT的输入序列的倒序由于序列不断进行奇偶分解所致，而DIF-FFT输出的倒序是由于序列前后对半分解后，其合成子序列正好对应着频率的奇偶部分。DIT和DIF在名称上也体现了这种不同点。

(2)蝶形张口大小和旋转因子次序的不同。从左至右来看，DIT-FFT蝶形张口由小到大，旋转因子由少到多；而DIF-FFT正好相反，蝶形张口由大到小，旋转因子由多到少。

(3)基本蝶形运算符号不同。图4.2.1中DIT-FFT蝶形不同于图4.3.1中DIF-FFT蝶形，DIT蝶形运算在频域进行，先乘旋转因子，后加减法；而DIF蝶形运算在时域进行，先加减法，后乘旋转因子。基本蝶形的不同才是两种FFT算法本质上的不同。

DIF-FFT和DIT-FFT的相互关系如下：

(1)基本蝶形运算符号的转置关系。若将图4.2.1中的DIT的基本蝶形进行转置，这里转置包括蝶形180°左右翻转、支路方向反向以及输入输出交换，就可以得到图4.3.1DIF的基本蝶形；同理，将DIF的基本蝶形加以转置，也可得到DIT的基本蝶形。

(2)蝶形运算流图的转置关系。对比图4.2.5DIT和图4.3.4DIF的两种蝶形运算流图，可以互相转置。这种特性有助于加深对两种FFT算法的理解和把握。 4.4快速傅里叶逆变换(IFFT)算法

本节研究离散傅里叶逆变换(IDFT)的快速算法，即快速傅里叶逆变换(InverseFastFourierTransform，IFFT)。比较IDFT与DFT的计算公式：(4.4.2)(4.4.1)可以看出，两者的差异在于变换对象(X(k)、x(n))、旋转因子(、)以及有无修正因子(1/N)的不同。若从公式计算角度来看，只需要DFT公式中的旋转因子换成，最后乘以1/N，就可以计算IDFT。

根据上述对IDFT和DFT的两个层面比较分析，IFFT的计算可以有两种方式：

(1)利用FFT蝶形运算流图。在FFT流图的基础上，按照变换对象、旋转因子以及有无修正因子三个方面进行适当修改，得到IFFT的蝶形运算流图。这部分内容在下面详细讨论。

(2)直接调用FFT程序。在用DFT公式计算IDFT的基础上，通过FFT程序来计算IFFT，这样，可以沿用FFT的蝶形运算流图。对IDFT表达式(4.4.1)两边取复共轭，可得

因此(4.4.3)(4.4.4)式(4.4.4)表明，利用FFT计算IDFT的过程如下：先将X(k)取共轭，然后直接调用FFT程序，计算结果取共轭后再乘以1/N。该方法虽然需要两次取共轭，但FFT和IFFT算法可以共有一个程序，使用很方便，在实际应用中大多采用这种方法。

下面重点讨论IFFT的第一种计算方式——蝶形运算流图。若基于图4.2.5所示的DIT-FFT蝶形运算流图，输入x(n)要换成X(k)，旋转因子变为，输出换成x(n)后再乘以1/N，这样得到图4.4.1所示的IFFT蝶形运算图。可以看出，原DIT-FFT的时间抽取变为IFFT的频率抽取，因此，该IFFT算法称为频率抽取法IFFT(DIF-IFFT)。图4.4.1

DIF-IFFT的蝶形运算流图(N=8)若基于图4.3.4所示的DIF-FFT蝶形运算流图，同样也需要修改三个方面：输入x(n)换成X(k)，旋转因子变为，输出换成x(n)后再乘以1/N，这样得到的IFFT蝶形运算图如图4.4.2所示的。图中，原DIF-FFT的频率抽取变为IFFT的时间抽取，因此，该IFFT算法称为时间抽取法IFFT(DIT-IFFT)。图4.4.2

DIT-IFFT的蝶形运算流图(N=8)在实际应用中，为了防止IFFT算法运行过程出现溢出，可以将分摊到每一级蝶形中。由于，因此，正好M级蝶形中每个蝶形输出均乘以。以图4.4.2的DIF-IFFT为例，经过防溢出处理后的蝶形运算图如图4.4.3所示。图4.4.3

DIT-IFFT的蝶形运算流图(N=8，防止溢出)下面关于溢出问题展开进一步讨论。

由于在实际DSP编程实现过程中，数值的表示存在着位数的限制，如定点DSP为16位(bit)，最高位表示符号位，可表示的整数值介于－32768～32767之间，浮点DSP则通常为32位。在进行FFT或者IFFT运算时，难免会出现溢出的问题。如定点DSP对一个单音进行频谱分析，其频域上的单频成分就非常高，很容易溢出。在这种情况下，可以考虑在某几级蝶形运算中再乘以1/2，避免数值溢出；同时，可以设置DSP的防溢出控制位，来限定最大值和最小值，避免数值溢出后出现正负数颠倒。

4.5

FFT算法的工程实现考虑

前面几节讨论的FFT和IFFT算法由于结构非常有规律，算法编程效率高，在实际数字信号处理中有着广泛的应用。下面将讨论FFT算法在工程实现时需要考虑的问题。

4.5.1旋转因子的生成

在FFT算法中，如何有效且快速地生成旋转因子是一个关键，直接涉及到蝶形中的乘法运算。以为例，k=0，1，2，…，－1，共有个复指数值，可展开为实部和虚部的形式：

可以看出，的生成包括余弦值和正弦值的计算，在编程实现时如何快速产生直接影响运算速度。

旋转因子的生成通常有两种方式：预先计算和直接计算。

1.预先计算

基本思想是预先计算出来所有N/2个值，以表的形式存储起来，供查找使用。该方法称为查表法，适用于FFT点数N已知的情况，且N不宜过大；否则，占用内存过多。从相位上看，相当于将2π分成N等份，取前N/2个弧度值2πk/N(k=0，1，2，…，N/2－1)。(4.5.1)由于式(4.5.1)包括余弦值和正弦值，按照常理可以分开制成余弦表和正弦表。但是，考虑到余弦和正弦在相位上相差π/2，即sin(α+π/2)=cosα，可以将余弦表和正弦表合成一张表，相位覆盖［0，2π］，共N个弧度值2πk/N(k=0，1，2，…，N－1)，只计算N个正弦值，构成一张正弦表。查表时，首先查找正弦值，然后向前搜索π/2，即N/4个弧度值，即可得到对应的余弦值。合成一张［0，2π］正弦表的好处是该表可同时用于除计算旋转因子之外的其它场合，如查表计算任意相位的余弦值或正弦值、产生数字频率等，这在数字信号处理领域特别是数字通信中常常出现。下面讨论如何查表得到任意相位的正弦值。假设相位α∈［0，2π)，单位为弧度，在N个弧度值2πk/N中α对应的位置可表示为

其中，［·］表示向下取整，即不超过的最大整数。当然，也可以采用四舍五入。通过查正弦表中nα位置，其对应的值即为正弦值。通过查表计算正弦值的实质是查找相邻相位及其正弦值来进行逼近，其精度高低与N的大小密切相关。若N值比较大，逼近误差较小，精度较高，但同时内存开销较大；若N值比较小，逼近误差较大，精度有限。(4.5.2)在实际中N值通常是确定的，为了进一步提高计算精度，基于给定的正弦表，可以通过内插获得更高的精度。如N=512，相邻相位差约为0.7°，比较小，连续正弦波可以用所有离散正弦值的直线连接来逼近。此时，利用α的左右相邻相位进行线性内插即可得到所需要的正弦值，其计算公式为(4.5.3)

2.直接计算

基本思想是运用Taylor级数展开式，直接计算余弦值和正弦值。余弦和正弦的Tylor公式为

直接计算的特点是精度很高，但是运算量比较大，适合于非实时性处理。(4.5.4)(4.5.5)4.5.2旋转因子的使用

一旦旋转因子生成好了，就可以直接参与蝶形乘法运算，按照FFT蝶形运算流图，逐级完成整个FFT计算。那么，是否还有必要专门讨论如何使用旋转因子呢？实际上还是非常有必要的，因为在实现FFT或DFT时，可能使用DSP、FPGA或者通用计算机，不同的实现方式有不同的特点，对旋转因子的使用和处理方式也不同。

就旋转因子本身而言，在计算FFT和DFT时，有很多特殊值，如

这样，可以不需要采用乘法甚至加减法，只需要简单的操作即可。鉴于通常情况下乘法的运算量要大于加法，针对这些特殊值作特别的编程处理，似乎可以降低运算量。但是，特殊编程处理增加了程序的复杂性，是否值得要取决于乘法和加法有多大差异，这与实现FFT或DFT的器件密切相关。

(1)若采用数字信号处理芯片(DSP)，由于DSP通常拥有专门的乘法指令，其指令运算时间与加减法指令一样，因此，完全可以进行乘法运算，没有必要针对旋转因子作特殊编程处理。

(2)若采用FPGA或通用计算机，由于乘法运算单元占用的硬件和运算资源通常要大大超过加减法，因此，根据实际FPGA型号或者计算机的资源状况，可在直接进行乘法运算和特殊编程处理之间进行合理选择。4.5.3实序列的FFT计算

在实际工作中，序列x(n)通常为实数序列，如模拟信号经过A/D采样后为实数字信号。如果直接按照FFT蝶形流图进行计算，需要将x(n)看做虚部为零的复序列，而由DFT的共轭对称性可知，实序列的FFT具有共轭对称性。如果能够利用实序列及其FFT的特点，可以降低FFT的运算量。

(1)一个N点FFT计算两个N点实序列的FFT。

基本做法是将两个N点实序列分别作为实部和虚部，构成N点复序列，再进行FFT。根据DFT的共轭对称性可知，实部FFT对应到复序列FFT的共轭对称部分，j和虚部的FFT对应到复序列FFT的共轭反对称部分。具体过程参见第3章中DFT共轭对称性的应用——两个实序列的