3. 5 快速傅里叶变换(FFT) 

3.5.1 DFT的计算量

离散傅里叶变换在实际应用中是非常重要的，利用它可以计算信号的频谱、功率谱和线性卷积等。但是，如果使用定义式(3.20)来直接计算DFT，当N很大时，即使使用高速计算机，所花的时间也太多。因此，如何提高计算DFT的速度，便成了重要的研究课题。1965年库利(Cooley)和图基(Tukey)在前人的研究成果的基础上提出了快速计算DFT的算法，之后，又出现了各种各样快速计算DFT的方法，这些方法统称为快速傅里叶变换(FastFourier Transform)，简称为FFT。FFT的出现，使计算DFT的计算量减少了两个数量级，从而成为数字信号处理强有力的工具。

快速傅里叶变换(FFT)是离散傅里叶变换(DFT)的快速算法。它是DSP领域中的一项重大突破，它考虑了计算机和数字硬件实现的约束条件、研究了有利于机器操作的运算结构，使DFT的计算时间缩短了1~2个数量级，还有效地减少了计算所需的存储容量，FFT技术的应用极大地推动了DSP的理论和技术的发展。

DFT的定义

在导出FFT算法之前，首先来估计一下直接计算DFT所需的计算量。

其中

将DFT定义式展开成方程组

将方程组写成矩阵形式

用向量表示

DFT复数运算量

.

计算一个X(k)(一个频率成分)值，运算量为

例k=1则

要进行N次复数乘法+(N-1)次复数加法

所以，要完成整个DFT运算，其计算量为：

N*N次复数相乘和N*(N-1)次复数加法

一次复数乘法换算成实数运算量

.

复数运算要比加法运算复杂，需要的运算时间长。

.

一个复数乘法包括4个实数乘法和2个实数加法。

(a+jb)(c+jd)=(ac-bd)+j(bc+ad)

因此计算N点的DFT共需要4N2次实数乘法和(2N2+2N·(N-1))次实数加法。当N很大时，这是一个非常大的计算量。

4次实数乘法

2次实数加法

例子

.

例1:当N=1024点时，直接计算DFT需要：

N2=220=1048576次，即一百多万次的复乘运算

这对实时性很强的信号处理(如雷达信号处理)来讲，它对计算速度有十分苛刻的要求-->迫切需要改进DFT的计算方法，以减少总的运算次数。

.

例2：石油勘探，24道记录，每道波形记录长度5秒，若每秒抽样500点/秒，

每道总抽样点数=500*5=2500点

24道总抽样点数=24*2500=6万点

DFT运算时间=N2=(60000)2=36*108次

FFT算法主要利用了WNk的两个性质：

(1)对称性，即

(2)周期性，即

r为任意整数。

例子

.

例:

利用以上特性,得到改进DFT直接算法的方法.

FFT算法是基于可以将一个长度为N的序列的离散傅里叶变换逐次分解为较短的离散傅里叶变换来计算这一基本原理的。这一原理产生了许多不同的算法，但它们在计算速度上均取得了大致相当的改善。

在本章中我们集中讨论两类基本的FFT算法。

第一类称为按时间抽取(Decimation-in-Time)的基2FFT算法，它的命名来自如下事实：在把原计算安排成较短变换的过程中，序列x(n)(通常看作是一个时间序列)可逐次分解为较短的子序列。

第二类称为按频率抽取(Decimation-in-Frequency)的基2FFT算法，在这类算法中是将离散傅里叶变换系数序列X(k)分解为较短的子序列。

前面两种算法特别适用于N等于2的幂的情况。

对于N为合数的情况，本章也将介绍两种处理方法。

将长序列DFT利用对称性和周期性分解为短序列DFT--方法

把N点数据分成二半：

其运算量为：

再分二半：

这样一直分下去，剩下两点的变换。

3. 5. 2时间抽选基2FFT算法(库里—图基算法) 

这种算法简称为时间抽选FFT算法，其基本出发点是，利用旋转因子WNk的对称性和周期性，将一个大的DFT分解成一些逐次变小的DFT来计算。

分解过程遵循两条规则：

①对时间进行偶奇分解；

②对频率进行前后分解。

设N＝2M，M为正整数。为了推导方便，取N＝23＝8，即离散时间信号为

按照规则(1)，将序列x(n)分为奇偶两组，一组序号为偶数，另一组序号为奇数，即

分别表示为：

根据DFT的定义

因为，所以上式变为

上式中的G(k)和H(k)都是N/2点的DFT。

(3.64)

按照规则(2)，将X(k)分成前后两组，即

由(3.64)表示的是N/2点DFT，前4个k值的DFT可表示为

后4个k值的X(k)表示为：

因为

所以

(3.66)

(3.65)

按照式(3.65)和式(3.66)可画出图3．15所示的信号流程图。

式(3.65)和式(3.66)把原来N点DFT的计算分解成两个N/2点DFT的计算。照此可进一步把每个N/2点DFT的计算再各分解成两个N/4点DFT的计算。具体说来，是把{x(0)，x(2)，x(4)，x(6)}和{x(1)，x(3)，x(5)，x(7)}分为{x(0)，x(4) | x(2)，x(6)}和{x(1)，x(5) | x(3)，x(7)}。这样，原信号序列被分成{x(0)，x(4) | x(2)，x(6) I x(1)，x(5) I x(3)，x(7)}4个2项信号。G(k)和H(k)分别计算如下:

(3.67)

(3.68)

(3.69)

(3.70)

这样，用式(3.67)～(3.70)4个公式就可计算图3.15中的两组N/2点DFT。图3.16所示的是其中一组G(k)的计算。

将图3.16与图3.15所示的信号流程图合并，便得到图3.17所示的信号流程图。

因为N=8，所以上图中N/4点的DFT就是2点的DFT，不能再分解了。

将2点DFT的信号流程图移入图3.17，得到图3.19所示的8点时间抽选的完整的FFT流程图。

返回

从图3.19中可看出几个特点：

(1)流程图的每一级的基本计算单元都是一个蝶形；

(2)输入x(n)不按自然顺序排列，称为“混序”排列，而输出X(k)按自然顺序排列，称为“正序”排列，因而要对输入进行“变址”；

(3)由于流程图的基本运算单元为蝶形，所以可以进行“同址”或“原位”计算。

3. 5. 3蝶形、同址和变址计算

1.  蝶形计算

任何一个N为2的整数幂(即N=2M)的DFT，都可以通过M次分解，最后成为2点的DFT来计算。M次分解构成了从x(n)到X(k)的M级迭代计算，每级由N/2个蝶形组成。图3.20表示了蝶形的一般形式表示。

其输入和输出之间满足下列关系：

从上式可以看出完成一个蝶形计算需一次复数乘法和两次复数加法。因此，完成N点的时间抽选FFT计算的总运算量为

大多数情况下复数乘法所花的时间最多，因此下面仅以复数乘法的计算次数为例来与直接计算进行比较。

直接计算DFT需要的乘法次数为αD=N2，于是有

例如，当N=1024时，则：205，即直接计算DFT所需复数乘法次数约为FFT的205倍。显然，N越大，FFT的速度优势越大。

表3. 2列出了不同N值所对应的两种计算方法的复数乘法次数和它们的比值。

2．同址(原位)计算

图3. 19包含log2N级迭代运算，每级由N/2个蝶形计算构成。蝶形计算的优点是可以进行所谓同址或原位计算。

现在来考察第一级的计算规律。设将输入x(0)，x(4)，x(2)，x(6)，x(1)，x(5)，x(3)，x(7)分别存入计算机的存储单元M(1), M(2), M(3),…，M(7)和M(8)中。首先，存储单元M(1)和M(2)中的数据x(0)和x(4)进入运算器并进行蝶形运算，流图中各蝶形的输入量或输出量是互不相重的，任何一个蝶形的二个输入量经蝶形运算后，便失去了利用价值，不再需要保存。这样，蝶形运算后的结果便可以送到M(1)和M(2)存储起来。类似地，M(3)和M(4)中的x(2)和x(6)进入运算器进行蝶形运算后的结果也被送回到M(3)和M(4)保存，等等。第二级运算与第一级类似，不过，M(1)和M(3)存储单元中的数据进行蝶形运算后的结果被送回M(1)和M(3)保存，M(2)和M(4)中的数据进行蝶形运算后送回M(2)和M(4)保存，等等。这样一直到最后一级的最后一个蝶形运算完成。

可以看出，每一级的蝶形的输入与输出在运算前后可以存储在同一地址(原来位置上)的存储单元中，这种同址运算的优点是可以节省存储单元，从而降低对计算机存储量的要求或降低硬件实现的成本。

蝶形运算的特点是，首先每一个蝶形运算都需要两个输入数据，计算结果也是两个数据，与其它结点的数据无关，其它蝶形运算也与这两结点的数据无关、因此，一个蝶形运算一旦计算完毕，原输入数据便失效了。这就意味着输出数据可以立即使用原输人数据结点所占用的内存。原来的数据也就消失了。输出、输人数据利用同一内存单元的这种蝶形计算称为同位(址)计算。

蝶形结

即蝶式计算结构也即为蝶式信号流图

上面频域中前/后半部分表示式可以用蝶形信号流图表示。

作图要素：

(1)左边两路为输入

(2)右边两路为输出

(3)中间以一个小圆表示加、减运算（右上路为相加输出、右下路为相减输出)

(4)如果在某一支路上信号需要进行相乘运算，则在该支路上标以箭头，将相乘的系数标在箭头旁。

(5)当支路上没有箭头及系数时，则该支路的传输比为1。

看出：用原位运算结构运算后，A(0)…A(7)正好顺序存放X(0)…X(7),可以直接顺序输出。

因子的分布

结论：每由后向前（m由M-1-->0级）推进一级，则此系数为后级系数中偶数序号的那一半。

3．变址计算

从图3. 19所示的流程图看出，同址计算要求输入x(n)是“混序”排列的。所谓输入为“混序”，并不是说输入是杂乱无章的，实际上它是有规律的。如果输入x(n)的序号用二进制码来表示，就可以发现输入的顺序恰好是正序输入的“码位倒置”，表3. 3列出了这种规律。

在实际运算中，按码位倒置顺序输入数据x(n)，特别当N较大时，是很不方便的。因此，数据总是按自然顺序输入存储，然后通过“变址”运算将自然顺序转换成码位倒置顺序存储。实现这种转换的程序可用图3. 21来说明。

图中用n表示自然顺序的标号，用l表示码位倒置的标号。当l=n时，x(n)和x(l)不必互相调换。当l≠n时，必须将x(l)和x(n)互相调换，但只能调换一次，为此必须规定每当l>n时，要将x(l)和x(n)相互调换，即把原来存放x(n)的存储单元中的数据调入存储x(l)的存储单元中，而把原来存储x(l)的存储单元中的数据调入到存储x(n)的存储单元中。

这样，按自然序输入的数据x(n)经过变址计算后变成了码位倒置的排列顺序，便可进入第一级的蝶形运算。

例子

以N=8为例：

3.5.4频率抽选基2FFT算法

频率抽选基2FFT算法简称为频率抽选，它的推导过程遵循两个规则：①对时间前后分；②对频率偶奇分。

设N＝2M，M为正整数。为推导方便，取N=23＝8。

首先，根据规则(1)，将x(n)分成前一半和后一半，即

x(n)＝{x(0), x(1), x(2), x(3), I x(4), x(5), x(6), x(7)} 

这样

(3. 72)

式(3. 72)虽然包含两个N/2点求和公式，但是在每个求和公式中出现的是WNkn，而不是WN/2kn，因此这两个求和公式都不是N/2点的DFT。如果合并这两个求和公式，并利用

则得：

其次，按规则(2)，将频率偶奇分，即

X(k)={X(0), X(2), X(4), X(6), | X(1), X(3), X(5), X(7)}

如果用X(2r)和X(2r+1)分别表示频率的偶数项和奇数项，则有

令

得到

上面两式所表示的是N/2的DFT。

在实际计算中，首先形成序列g(n)和h(n)，然后计算h(n)WNn，最后分别计算g(n) 和h(n)WNn的N/2点DFT，便得到偶数输出点和奇数输出点的DFT。计算流程图如图3. 24所示。

由于N是2的整数幂，所以N/2仍然是偶数。这样，可以将N/2点DFT的输出再分为偶数组和奇数组，也就是将N/2点的DFT计算进一步分解为两个N/4点的DFT计算，其推导过程如下。

将g(n)分为前后两组，得到

因为

所以

对频率再进行偶奇分，则得频率的偶数项为

频率的奇数项为

通过类似的推导可得

上面4式所表示的计算都是N/4点的DFT计算，从而得到图3.25所示的形式。

与时间抽选的FFT算法一样，图3.27所示的流程图的基本运算也是蝶形运算，但是与时间抽选中的蝶形单元(图3.20)有所不同。

图3. 28所示的是频率抽选FFT算法的蝶形单元的一般化的流程图

显然，一个蝶形的计算包括1次复数乘法和2次复数加法。图3. 27所示的整个流程图共有log2N级迭代运算，每级有N/2个蝶形。因此，总计算量为

DIF与DIT比较1

.

相同之处：

.

（1）DIF与DIT两种算法均为原位运算。

.

（2）DIF与DIT运算量相同。

.

它们都需要

DIF与DIT比较2

.

不同之处：

（1）DIF与DIT两种算法结构倒过来。

DIF为输入顺序，输出乱序。运算完毕再运行“二进制倒读”程序。

DIT为输入乱序，输出顺序。先运行“二进制倒读”程序，再进行求DFT。

（2）DIF与DIT根本区别：在于蝶形结不同。

DIT的复数相乘出现在减法之前。

DIF的复数相乘出现在减法之后。

3. 5. 5IFFT的计算方法

FFT算法同样可以应用于IDFT的计算，称为快速傅里叶反变换，简写为IFFT。前述DFT和IDFT公式为

比较上面两式，可以看出，只要把DFT公式中的系数改为，并乘以系数1/N，就可用FFT算法来计算IDFT，这就得到了IFFT的算法。

当把时间抽选FFT算法用于IFFT计算时，由于原来输入的时间序列x(n)现在变为频率序列X(k)，原来是将x(n)偶奇分的，而现在变成对X(k)进行偶奇分了，因此这种算法改称为频率抽选IFFT算法。类似地，当把频率抽选FFT算法用于计算IFFT时，应该称为时间抽选IFFT算法。

在IFFT计算中经常把常量1/N分解成M个1/2连乘，即1/N＝(1/2)M，并且在M级的迭代运算中，每级的运算都分别乘上一个1/2因子。图3.29表示的是时间抽选IFFT流程图。

IFFT的基本蝶形运算