第 03章 离散傅里叶变换及其快速算法邹江
zoujiang@public.wh.hb.cn
3,5 快速傅里叶变换 (FFT)
3.5.1 DFT的计算量离散傅里叶变换在实际应用中是非常重要的,利用它可以计算信号的频谱、功率谱和线性卷积等。但是,如果使用定义式 (3.20)来直接计算 DFT,当 N很大时,即使使用高速计算机,所花的时间也太多。因此,如何提高计算 DFT的速度,便成了重要的研究课题。
1965年库利 (Cooley)和图基 (Tukey)在前人的研究成果的基础上提出了快速计算 DFT的算法,之后,又出现了各种各样快速计算 DFT的方法,这些方法统称为快速傅里叶变换 (Fast Fourier Transform),简称为
FFT。 FFT的出现,使计算 DFT的计算量减少了两个数量级,从而成为数字信号处理强有力的工具。
快速傅里叶变换 (FFT)是离散傅里叶变换 (DFT)的快速算法。它是
DSP领域中的一项重大突破,它考虑了计算机和数字硬件实现的约束条件、研究了有利于机器操作的运算结构,使 DFT的计算时间缩短了 1~2个数量级,还有效地减少了计算所需的存储容量,FFT技术的应用极大地推动了 DSP的理论和技术的发展。
DFT的定义在导出 FFT算法之前,首先来估计一下直接计算 DFT所需的计算量。
其中将 DFT定义式展开成方程组将方程组写成矩阵形式用向量表示从矩阵形式表示可以看出,由于计算一个 X(k)值需要 N次复乘法和
(N-1)次复数加法,因而计算 N个 X(k)值,共需 N2次复乘法和 N(N-1)次复加法。每次复乘法包括 4次实数乘法和 2次实数加法,每次复加法包括 2次实数加法,因此计算 N点的 DFT共需要 4N2次实数乘法和
(2N2+2N·(N-1))次实数加法。当 N很大时,这是一个非常大的计算量。
FFT算法主要利用了 WNk的两个性质:
(1)对称性,即
(2)周期性,即用复数表示:
r为任意整数。
FFT算法是基于可以将一个长度为 N的序列的离散傅里叶变换逐次分解为较短的离散傅里叶变换来计算这一基本原理的。这一原理产生了许多不同的算法,但它们在计算速度上均取得了大致相当的改善。
在本章中我们集中讨论两类基本的 FFT算法。
第一类 称为按时间抽取 (Decimation-in-Time)的基 2FFT算法,它的命名来自如下事实:在把原计算安排成较短变换的过程中,
序列 x(n)(通常看作是一个时间序列 )可逐次分解为较短的子序列。
第二类称为按频率抽取 (Decimation-in-Frequency)的基 2FFT算法,
在这类算法中是将离散傅里叶变换系数序列 X(k)分解为较短的子序列。
前面两种算法特别适用于 N等于 2的幂的情况。
对于 N为合数的情况,本章也将介绍两种处理方法。
3,5,2 时间抽选基 2FFT算法 (库里 — 图基算法 )
这种算法简称为时间抽选 FFT算法,其基本出发点是,利用旋转因子 WNk的对称性和周期性,将一个大的 DFT分解成一些逐次变小的 DFT来计算。
分解过程遵循两条规则:
①对时间进行偶奇分解;
②对频率进行前后分解。
设 N= 2M,M为正整数。为了推导方便,取 N= 23= 8,即离散时间信号为按照规则 (1),将序列 x(n)分为奇偶两组,一组序号为偶数,另一组序号为奇数,即分别表示为:
根据 DFT的定义因为 WN2=WN/21,所以上式变为上式中的 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)保存,等等。这样一直到最后一级的最后一个蝶形运算完成。
可以看出,每一级的蝶形的输入与输出在运算前后可以存储在同一地址 (原来位置上 )的存储单元中,这种同址运算的优点是可以节省存储单元,从而降低对计算机存储量的要求或降低硬件实现的成本。
蝶形运算的特点是,首先每一个蝶形运算都需要两个输入数据,计算结果也是两个数据,与其它结点的数据无关,其它蝶形运算也与这两结点的数据无关、因此,一个蝶形 运算一旦计算完毕,原输入数据便失效了。这就意味着输出数据可以立即使用原输 人数据结点所占用的内存。原来的数据也就消失了。输出、输人数据利用同一内存单 元的这种蝶形计算称为同位 (址 )计算。
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)经过变址计算后变成了码位倒置的排列顺序,便可进入第一级的蝶形运算。
最后介绍一下时间抽选 FFT算法的另外一些形式的流程图。对于任何流程图,只要保持 各节点所连支路及其传输系数不变,则不论节点位置怎样排列,所得到的流程图总是等效的,因而都能得到
DFT的正确结果,只是数据的提取和存储次序不同而已。
把 图 3,19中与 x(4)水平相邻的所有节点和与 x(1)水平相邻的所有节点交换,把与 x(6)水平相邻的所有节点和与 x(3)水平相邻的所有节点交换,而与 x(2),x(5)和 x(7)水平相邻各节点位置不变,就可以从图 3,
19得到图 3.22。图 3.22与图 3.19的区别只是节点的排列不同,而支路传输比,即 WN的 各次幂保持不变。显然图 3.22所示流程图的输入是正序 (自然顺序 )排列的,输出是码位倒置 排列的,所以输出要进行变址计算。图 3,22所示的流程图相当于最初由库利和图基给出的时间抽选算法。
另一种形式的流程图是将节点排列成输入 和输出两者都是正序排列,但这类流程图不能进行同址计算,因而需要两列长度为 N的复数存储器。
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个蝶形。因此,
总计算量为频率抽选的 FFT算法的计算量与时间抽选 FFT算法的计算量是一样。
与时间抽选算法一样,频率抽选 FFT算法也具有同址 (原位 )计算的优点。但是,与时间抽选不同的是,频率抽选 FFT算法的信号输入为正序排列,输出为码位倒置排列,因此输出要进行变址计算。
3,5,5 IFFT的计算方法
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流程图。