1、,第四章 快速傅里叶变换 （FFT）,第十二讲,本章内容：,介绍傅里叶变换的一些快速算法快速算法的思想 根据原是变换定义的运算规律，及其中某些算子的特殊性，找出减少乘法和加法运算次数的有效途径，实现原始变换的各种高效算法。,本讲学习目标,了解直接计算N点DFT的运算量了解减少运算量的基本途径理解按时间抽取的基-2FFT算法的算法原理、运算流图、所需计算量和算法特点了解按时间抽取的基-2FFT算法的编程思想,本章作业练习,P127: 4.14.24.44.5,第四章 快速傅里叶变换,FFT: Fast Fourier Transform1965年，Cooley, Tukey机器计算傅里叶级数的一
2、种算法4.2 基-2FFT算法,4.2.1直接计算DFT的特点及减少运算量的基本途径,1. DFT的运算量及特点,运算量,FFT算法分类:,时间抽选法DIT: Decimation-In-Time频率抽选法DIF: Decimation-In-Frequency,4.2.2、时间抽取法基-2FFT算法基本思想 （基-2 Decimation-In-Time FFT）,1、算法原理设序列点数 N = 2M，M 为整数。 若不满足，则补零,将序列x(n)按n的奇偶分成两组：,N为2的整数幂的FFT算法称基-2FFT算法。,n为偶数时: n为奇数时:,则x(n)的DFT:,其中,2.两点结论: (1
3、) X (k),X (k)均为N/2点的DFT。 (2) X(k)=X (k)+W X (k)只能确定出 X(k)的k= 个；即前一半的结果。,1 2,1 2,k,N,同理, 这就是说,X1(k),X2(k)的后一半,分别 等于其前一半的值。,3.X(k)的后一半的确定,由于 （周期性）,所以:,可见,X(k)的后一半，也完全由X1(k), X2 (k)的前一半所确定。 * N点的DFT可由两个N/2点的DFT来计算。,又由于 ,所以,实现上式运算的流图称作蝶形运算,前一半,4.蝶形运算,后一半,（N/2个蝶形),(前一半),(后一半),1 1,1,1,-1,由X1(k)、X 2(k)表示X(
4、k)的运算是一种特殊的运算-碟形运算,图4.2.1 蝶形运算符号,x1(0)=x(0) x1(1)=x(2) N/2点 x1(2)=x(4) DFT x1(3)=x(6) x2(0)=x(1) x2(1)=x(3) N/2点 x2(2)=x(5) DFT x2(3)=x(7), ,X1(0),X1(1),X1(2),X1(3),X2(0),X2(1),X2(2),X2(3),W,N,2,W,N,1,W,N,0,W,N,3,-1,-1,-1,-1,X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7),5.对X1 (k)和X2 (k)进行蝶形运算,前半部为 X(0)- X
5、(3),后半部分为X(4)- X(7) 整个过程如下图所示:,6. N/2分解后的运算量：,运算量减少了近一半,由于N=2 M ,所以 N/2仍为偶数,可以进 一步把每个N/2点的序列再按其奇偶部分 分解为两个N/4的子序列。例如，n为偶数时的 N/2点分解为:,(二) N/4点DFT,进行N/4点的DFT,得到,(偶中偶),(偶中奇),同理可将奇数序号组成的N/2点序列进行分解得：,其中：,X2的偶数序列,X2的偶数序列,这样逐级分解，直到2点DFT当N = 8时，即分解到X3(k)，X4(k)，X5(k)，X6(k)，k = 0, 1,不再做DFT,(0) (4) (2) (6) (1)
6、(5) (3) (7),W,N,0,W,N,0,W,N,0,W,0,N,-1,-1,-1,-1,-1,-1,-1,-1,N,0,N,1,N,2,N,3,-1,-1,-1,-1,X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7),因此,8点DFT的FFT的运算流图如下,4.2.3 基2DIT-FFT算法与直接DFT运算量的比较,当N = 2M时，共有M级蝶形，每级N / 2个蝶形，每个蝶形有1次复数乘法，2次复数加法。,复数乘法：,复数加法：,与DFT 比较,4.2.4 DIF-FFT的运算规律 及编程思想,1）原位计算 DIT-FFT的运算过程很有规律，共进行M级运
7、算，每级由N/2个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，与其它蝶形运算无关。 这样,蝶形运算的两个输出值仍可放回蝶形运算的两个输入所在的存储器中,这种利用同一存储单元存储蝶形计算输入、输出的方法即为原位运算。每一级(列)有N/2个蝶形运算,所以只需N个存储单元，可以节省存储单元。,运算规律,(0)=A0(0) A1(0) A2(0) A3(0)=X(0) (4)=A0(1) A1(1) A2(1) A3(1)=X(1) (2)=A0(2) A3(2)=X(2) (6)=A0(3) A3(3)=X(3) (1)=A0(4) A1(4) A2(4) A3(4)=X(4)
8、 (5)=A0(5) A3(5)=X(5) (3)=A0(6) A3(6)=X(6) (7)=A0(7) A1(7) A2(7) A3(7)=X(7),W,W,W,W,N,0,N,0,N,0,N,0,-1,-1,-1,-1,W,W,W,W,N,0,N,2,N,0,N,2,-1,-1,-1,-1,W,W,W,W,N,N,N,N,0,1,2,3,.,.,.,.,原位运算的图示,输入数据、中间运算结果和最后输出均用同一存储器。,L=1,L=2,L=3,开始时，输入序列的N个数放于此N个存储器内，倒序重排后仍存于这N个存储器中，每一次迭代运算后的结果也仍然存于这N个存储器中，整个运算完成后，这N个存储
9、器中即为所求的频谱序列X(k) （k = 0、1、.、 N-1）。这就是所谓的同址计算，这样可以大大节约存储元件。,N点DITFFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子WpN，称其为旋转因子，p称为旋转因子的指数。,2）旋转引子的变化规律,3）蝶形运算规律,对N = 2M点FFT，输入倒位序，输出自然顺序，设第L级运算每个蝶形的两节点距离为 B行，则第L级运算：,旋转因子的确定仅与指数P有关，当L一定时J可以确定，进而可以确定指数P。对于同一P，参与本级蝶形运算的次数为 ，每级第一次调用 的蝶形结第一个输入节点的位置为J，第二次调用 的蝶形结第一个输入节点的位置为J+2L，即
10、以2L为步进，搜索下一蝶形运算第一个输入节点的位置。以此类推，直至做完 个蝶形运算。,同一级中，同一P的蝶形计算完成后，计算下一个P值对应的蝶形运算，直至完成本级所有蝶形运算。,DIT-FFT原位计算步骤,1.确定L；2.求出第L级的 个 ；3.对每一个 完成其所参与的 个蝶形运算，4. 重复步骤13完成M级蝶形运算。,以2L为步进参与本级同一旋转因子对应的其他蝶形运算,作为练习请写出L=3时的运算过程。,4） 编程思想及程序框图在一般情况下，进行FFT运算的序列至少都有几百点的长度，因此需要编制FFT算法程序以便能够利用计算机来快速进行计算。输入倒位序，输出顺序的DIT-FFT的编程思想，
11、N必须等于2的正整数幂，FFT的计算程序可以分为两部分：一部分是倒序重排，另一部分是用三层嵌套的循环来完成M=log2N次迭代。,三层循环的功能是：最里的一层循环完成相同WNP 的蝶形运算，中间的一层循环完成因子WNP的变化，而最外的一层循环则是完成M次迭代过程。,循环L,确定输入数据所在的位置,循环J,L、P、B确定后进行蝶形计算,5） 序列的倒序 由DIT-FFT的规律可知,输出X(k)按正常顺序排列在存储单元,而输入是按顺序: 这种顺序称作倒位序,即二进制数倒位。在编程时需完成倒位序，才能执行原位计算。,倒位序由奇偶分组造成,以N=8为例 说明如下:,0 0 0 0 0 0 0 0 1
12、0 0 1 1 0 0 4 2 0 1 0 0 1 0 2 3 0 1 1 1 1 0 6 4 1 0 0 0 0 1 1 5 1 0 1 1 0 1 5 6 1 1 0 0 1 1 3 7 1 1 1 1 1 1 7,自然顺序n 二进制n n n 倒位序二进制n n n 倒位顺序n,2 1 0 0 1 2,权重：,X(I),X(J),最低位加1,若最高位为0，二进制最高位加1，对应的十进数加N/2=4,若最高位为1，最高位置0，同时J=J-N/2，若次高位为0次高位加1,J=J+N/4,若次高位为1，次高位置0同时J=J-N/4，次次高位加1,J=J+N/8,顺序I起始及终止序号为：16倒序
13、J起始序号为：N/2=4 当IJ 时，A(I)和A(J)的内容调换,形成倒序J后，将原存储器放的输入序列重新按倒序排列。,X(I),X(J),I=J时不需要交换,图4.2.9 倒序程序框图,确定J的起始序号及I的终止序号,K=N/2,最高位为0否？,J=J+N/2,最高位为1，将最高位置0，J=J-N/2；次高位加1，K=N/4，,倒序重排的程序是一段经典程序，它以巧妙的构思、简单的语句用高级编程语言来完成了倒序重排的功能。下面是一段用FORTRAN语言编写的倒序重排程序。,N=2*M （表示N=2M ，M是输入的正整数） LH=N/2 （LH是一个整数变量） N1=N-2 （N1也是一个整数
14、变量） J=1 （对变量J赋初值） 100 DO 7 I=1,N1(循环开始，到语句7结束； 循环变量I从1开始，到N1结束，步长为1),IF (I.GE.J) GOTO 5（如果IJ，就转移到语句5） （将输入序列中序号互为倒序 的两个数值交换位置）,5 K=N2 6 IF(K.GE.J) GOTO 7 (如果KJ，就转移到语句7) J=J-kK=K/2 GOTO 6 （转移到语句6）7 J=J+K ,附：DIT算法的其他形式流图,输入倒位序输出自然序输入自然序输出倒位序输入输出均自然序相同几何形状输入倒位序输出自然序输入自然序输出倒位序,用Matlab方法计算信号的DFT时，是用函数fft(x,N) 和ifft(x.N)。对于这两个函数，如果N为2的正整数幂，则可以得到本章中介绍的基2 FFT快速算法；如果N既不是2的正整数幂，也不是质数，则函数将N分解成质数，得到较慢的混合基 FFT算法；如果 N 为质数，则fft函数采用原来的 DFT 算法。,附：利用Matlab计算FFT,
