第4章FFT.ppt

1、第四章 快速傅里叶变换(FFT),Fast Fourier Transform,内容简介,基2FFT算法 时域抽选基2FFT 频域抽取法FFT 进一步减少运算量的措施,4.2 基2FFT算法,4.2.1 DFT的计算量,DFT的定义,计算量：复数乘法 N2 次，复数加法 N(N-1) 实数乘法 4N2 次，实数加法 4N(N-1) x(n)也可为复数。Wkn可事先计算好,直接计算DFT的计算量为N2数量级，需要改进算法提高计算效率,算法改进的途径： 把长序列分解成短序列计算，然后再组合出结果； 利用计算中的对称性，减少计算量。,r为任意整数。,典型的FFT算法,1965年,库利(J.W.Coo

2、ley)和图基(J.W.Tukey)在Mathematics of Computation杂志发表论文机器计算傅立叶级数的一种算法，首先提出FFT算法,这就是库利图基算法（时间抽取基2FFT算法）。对于N点DFT,仅需(N/2)log2N 次复数乘法运算.例如N=1024=210 时，需要(1024/2)log2 210=512*10=5120次。5120/1048576=4.88% ,速度提高20倍。 在此基础上有法长出许多改进算法如:桑德图基算法（ 频率抽取基2FFT算法）；WFTA算法（小N素数组合算法）；CZT算法（线性调频Z变换算法）；ZFFT算法（局部频谱细化算法）等等。,基本出发

3、点：利用旋转因子WNk的对称性和周期性，将一个大的DFT分解成一些逐次变小的DFT来计算。 分解过程遵循两条规则： 对时间进行偶奇分解； 对频率进行前后分解。,设N2M，M为正整数。为了推导方便，取N238，即离散时间信号为,4.2.2时间抽选基2FFT算法(库里图基算法),按照规则(1)，将序列x(n)分为奇偶两组，一组序号为偶数，另一组序号为奇数，即,分别表示为：,根据DFT的定义,因为 ,所以上式变为,上式中的X1(k)和X2(k)分别为x1(r)和x2(r) 的N/2点的DFT。,(4.2.4),按照规则(2)，将X(k)分成前后两组，即,由(4.2.4)表示的是N/2点DFT，前4个

4、k值的DFT可表示为,因为X1(k)和X2(k)均以N/2为周期，后4个k值的X(k)表示为：,因为,所以,(4.2.8),(4.2.7),图4.2.1 蝶形运算符号,N点DFT分解为两个N/2点DFT和式(4.2.7)和(4.2.8)的蝶形运算，如图所示：,按照式(4.2.7)和式(4.2.8)可画出N点DFT经过一次奇偶抽取分解后的信号流程图。,图4.2.2 N点DFT的一次时域抽取分解图(N=8),一次奇偶抽取分解后的运算量分析： 一个蝶形运算：一次复数乘法和两次复数加法运算 根据图4.2.2，一次分解后，1个N点DFT需计算两个N/2点DFT和N/2个蝶形运算。总的复数乘法次数：,复数

5、加法次数为,因此，一次分解后，运算量减少近一半。且N=2，N/2仍为偶数，可再进一步分解。,将x1(r)和x2(r)按奇偶分解成两个N/4长的子序列x3(l)和x4(l)，即把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项信号。X1(k)和X2(k)分别计算如下:,(4.2.9),式中,同理，由X3(k)和X4(k)的周期性和Wm N/2的对称 性 Wk

6、+N/4 N/2= - Wk N/2 最后得到：,(4.2.10),用同样的方法可计算出,(4.2.11),其中,图4.2.3 N点DFT的第二次时域抽取分解图(N=8),将2点DFT的信号流程图移入图4.2.3，得到图4.2.4所示的8点时间抽选的完整的FFT流程图。,图4.2.4 N点DITFFT运算流图(N=8),从图4.2.4中可看出几个特点： (1) 流程图的每一级的基本计算单元都是一个蝶形； (2) 输入x(n)不按自然顺序排列，称为“混序”排列，而输出 X(k)按自然顺序排列，称为“正序”排列，因而要对输入进行“变址”； (3) 由于流程图的基本运算单元为蝶形，所以可以进行“同址

7、”或“原位”计算。,4.2.3DIT-FFT算法与直接计算DFT运算量的比较,由DIT-FFT算法的分解过程，当N=2M时，其运算流图有M级蝶形，每一级都由N/2个蝶形运算构成，每一级运算都需要N/2次复数乘法和N次复数加法。M级运算总的复数乘法次数为,复数加法次数：,而直接计算DFT的复数乘法为N2次，复数加为N(N-1)次，当N1时，DIT-FFT算法比直接计算DFT的运算次数大大减少。,图4.2.5 FFT算法与直接计算DFT所需乘法次数的比较曲线,1同址(原位)计算,图4.2.4包含log2N级迭代运算，每级由N/2个蝶形计算构成。蝶形计算的优点是可以进行所谓同址或原位计算。,现在来考

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

9、0)和A(2)存储单元中的数 据进行蝶形运算后的结果被送回A(0)和A(2)保存，A(1)和A(3)中的数据进行蝶形运算 后送回A(1)和A(3)保存，等等。这样一直到最后一级的最后一个蝶形运算完成。,4.2.4DIT-FFT的运算规律及编程思想,可以看出， 每一级的蝶形的输入与输出在运算前后可以存储在同一地址(原来位置上)的存储单元中，这种同址运算的优点是可以节省存储单元，从而降低对计算机存储量的要求或降低硬件实现的成本。,蝶形运算的特点是，首先每一个蝶形运算都需要两个输入数据，计算结果也是两个数据，与其它结点的数据无关，其它蝶形运算也与这两结点的数据无关、因此，一个蝶形 运算一旦计算完毕，

10、原输入数据便失效了。这就意味着输出数据可以立即使用原输 人数据结点所占用的内存。原来的数据也就消失了。输出、输人数据利用同一内存单 元的这种蝶形计算称为同位(址)计算。,2. 旋转因子的变化规律,称为旋转因子，p称为旋转因子的指数。用L表示从左到右的运算级数，L=1,2,，M ，则第L级共有2L-1个不同的旋转因子。 N=8时各级旋转因子表示：,对N=2M的一般情况，第L级的旋转因子为,(4.2.12),(4.2.13),因为,此即确定第L级运算的旋转因子。,3. 蝶形运算规律,观察图4.2.4，序列x(n)经时域抽选（倒序）后存入数组A，若蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶

11、形运算可表示为：,其中,4. 编程思想及程序框图,图4.2.6 DITFFT运算和程序框图,流程：从输入端开始，逐级进行，共M级。 第L级运算中，依次求出B个不同的旋转因子，每求出一个旋转因子，计算完它对应的所有2M-L个蝶形，采用三重循环程序实现。,5序列的倒序,从图4.2.4所示的流程图看出，同址计算要求输入x(n)是“混序”排列的。所谓输入为“混 序”，并不是说输入是杂乱无章的，实际上它是有规律的。如果输入x(n)的序号用二进制码来 表示，就可以发现输入的顺序恰好是正序输入的“码位倒置”，表3. 3列出了这种规律。,在实际运算中，按码位倒置顺序输入数据x(n)，特别当N较大时，是很不方便

12、的。因此，数据总是按自然顺序输入存储，然后通过“变址”运算将自然顺序转换成码位倒置顺序存储。实现这种转换的程序可用图4.2.8来说明。,说明：自然顺序数I增加1，是在顺序数的二进制数最低位加1，逢2向高位进位，而倒序数则是在M位二进制数最高位加1，逢2向低位进位。,x(n),x(l),图中用n表示自然顺序的标号，用l表示码位倒置的标号。当l=n时，x(n)和x(l)不必互相调换。当ln时， 必须将x(l)和x(n)互相调换，但只能调换一次，为此必须规定每当ln时，要将x(l)和x(n)相互调换，即把原来存放x(n)的存储单元中的数据调入存储x(l)的存储单元中，而把原来存储x(l)的存储单元中

13、的数据调入到存储x(n)的存储单元中。,这样，按自然序输入的数据x(n)经过变址计算后变成了码位倒置的排列顺序，便可进入第一级的蝶形运算。,图4.2.9 倒序程序框图,时间抽选FFT算法的另外一些形式的流程图。对于任何流程图，只要保持 各节点所连支路及其传输系数不变，则不论节点位置怎样排列，所得到的流程图总是等效的，因而都能得到DFT的正确结果，只是数据的提取和存储次序不同而已。,把图4.2.4中与x(4)水平相邻的所有节点和与x(1)水平相邻的所有节点交换，把与x(6)水平相邻的所有节点和与 x(3)水平相邻的所有节点交换，而与x(2)、x(5)和x(7)水平相邻各节点位置不变，就可以从图4

14、.2.4得到图3.22。图3.22与图4.2.4的区别只是节点的排列不同，而支路传输比，即WN的 各次幂保持不变。显然图3.22所示流程图的输入是正序(自然顺序)排列的，输出是码位倒置 排列的，所以输出要进行变址计算。图3. 22所示的流程图相当于最初由库利和图基给出的时 间抽选算法。,另一种形式的流程图是将节点排列成输入 和输出两者都是正序排列，但这类流程图不能进行同址计算，因而需要两列 长度为N的复数存储器。,另一种常用的快速算法频率抽取法FFT的推导过程遵循两个规则：对时间前后分；对频率偶奇分。 设N2M，M为正整数。为推导方便，取N=238。 首先，根据规则(1)，将x(n)分成前一半

15、和后一半，即 x(n)x(0), x(1), x(2), x(3), I x(4), x(5), x(6), x(7) 这样,4.2.5频率抽取法基2FFT(DIF-FFT),利用,则得：,再按规则(2)，将频率偶奇分，即,X(k)=X(0), X(2), X(4), X(6), | X(1), X(3), X(5), X(7),如果用X(2m)和X(2m+1)分别表示频率的偶数项和奇数项，则有,(4.2.14),令,(4.2.15),将x1(n)和x2(n)分别代入(4.2.14)和(4.2.15)式，可得,结论： X(k)按奇偶k值分为两组，其偶数组是x1(n)的N/2点DFT，奇数组则是

16、x2(n)的N/2点DFT，x1(n)、x2(n)与x(n)的关系可用图4.2.10所示的蝶形运算流图符合表示。,图4.2.10 DIFFFT蝶形运算流图符号,图4.2.11 DIFFFT一次分解运算流图(N=8),由于N是2的整数幂，所以N/2仍然是偶数。这样，可以将N/2点DFT的输出再分为偶数组和奇数组，也就是将N/2点的DFT计算进一步分解为两个N/4点的DFT计算，其推导过程如下。,将x1(n)分为前后两组，得到,因为,所以,对频率再进行偶奇分，则得频率的偶数项为,频率的奇数项为,通过类似的推导可得,上面4式所表示的计算都是N/4点的DFT计算，从而得到图4.2.12所示的形式。,图

17、4.2.12 DIFFFT二次分解运算流图(N=8),图4.2.13 DIFFFT运算流图(N=8),比较频率抽选和时间抽选的蝶形运算单元，频率抽选的蝶形为,显然，一个蝶形的计算包括1次复数乘法和2次复数加法。图4.2.13 所示的整个流程图共有log2N级迭代运算，每级有N/2个蝶形。因此，总计算量为,频率抽选的FFT算法的计算量与时间抽选FFT算法的计算量是一样。,与时间抽选算法一样，频率抽选FFT算法也具有同址(原位)计算的优点。但是，与时间抽选不同的是，频率抽选FFT算法的信号输入为正序排列，输出为码位倒置排列，因此输出要进行变址计算。,FFT算法同样可以应用于IDFT的计算，称为快速

18、傅里叶反变换，简写为IFFT。前述DFT和IDFT公式为,比较上面两式，可以看出，只要把DFT公式中的系数 改为 ，并乘以系数1/N，就可用FFT算法来计算IDFT，这就得到了IFFT的算法。 当把时间抽选FFT算法用于 IFFT计算时，由于原来输入的时间序列x(n)现在变为频率序列X(k)，原来是将x(n)偶奇分的，而现在变成对X(k)进行偶奇分了，因此这种算法改称为频率抽选IFFT算法。类似地，当把频率抽选FFT算法用于计算IFFT时，应该称为时间抽选IFFT算法。,4.2.6IFFT的计算方法,在IFFT计算中经常把常量1/N分解成M个1/2连乘，即1/N(1/2)M，并且在M级的迭代运

19、算中，每级的运算都分别乘 上一个1/2因子。图4.2.14表示的是时间抽选IFFT流程图。,图4.2.14 DITIFFT运算流图(防止溢出),4.3 进一步减少运算量的措施,4.3.1多类蝶形单元运算,DITFFT运算量：N=2M点FFT共需要MN/2次复数乘法。,对N=2M的一般情况，第L级的旋转因子为,当L=1时，只有一种旋转因子W0N=1，第一级不需要乘法运算。 当L=2时，有两个旋转因子，W0N =1，WN/4N = - j，第二级也不需乘法运算 无关紧要的旋转因子：值为1，-1，j，-j的选择因子。,a) 除去第一、二两级后，所需复数乘法次数应是,b) 再考虑无关紧要的旋转因子，因同一旋转因子对应2M-L=N/2L个蝶形运算，从L=3至L=M共减少复数乘法次数为,DIT-FFT的复乘法次数减为：,(4.3.