FFT算法介绍和编程1

,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第四章 快速傅里叶变换 （,FFT,）,第十二讲,本章内容：,介绍傅里叶变换的一些快速算法,快速算法的思想,根据原是变换定义的运算规律，及其中某些算子的特殊性，找出减少乘法和加法运算次数的有效途径，实现原始变换的各种高效算法。,本讲学习目标,了解直接计算,N,点,DFT,的运算量,了解减少运算量的基本途径,理解按时间抽取的基,-2FFT,算法的算法原理、运算流图、所需计算量和算法特点,了解,按时间抽取的基,-2FFT,算法的编程思想,本章作业练习,P127:,4.1,4.2,4.4,4.5,第四章 快速傅里叶变换,FFT:Fast Fourier Transform,1965,年，,Cooley,Tukey,机器计算傅里叶级数的一种算法,4.2,基,-2FFT,算法,4.2.1,直接计算,DFT,的特点及减少运算量的基本途径,1.DFT,的运算量及特点,运算量,复数乘法,复数加法,一个,X,(,k,),N,N,1,N,个,X,(,k,),(,N,点DFT),N,2,N,(,N,1),实数乘法,实数加法,一次复乘,4,2,一次复加,2,一个,X,(,k,),4,N,2,N,+2(,N,1)=2(2,N,1),N,个,X,(,k,),(,N,点DFT),4,N,2,2,N,(2,N,1),FFT,算法分类,:,时间抽选法,DIT:Decimation-In-Time,频率抽选法,DIF:Decimation-In-Frequency,4.2.2,、时间抽取法基,-2FFT,算法基本思想 （基,-2 Decimation-In-Time FFT,）,1,、算法原理,设序列点数,N,=2,M,，,M,为整数。,若不满足，则补零,将序列,x,(,n,),按,n,的奇偶分成两组：,N,为,2,的整数幂的,FFT,算法称基,-2FFT,算法。,n,为偶数时,:,n,为奇数时,:,则,x,(,n,)的DFT:,其中,2.,两点结论:,(1),X,(k),X,(k),均为,N/2,点的,DFT,。,(2)X(k)=X,(k)+W,X,(k),只能确定出,X(k),的,k=,个；,即,前一半,的结果。,1 2,1 2,k,N,同理,这就是说,X,1,(k),X,2,(k),的后一半,分别,等于其前一半的值。,3.X(k),的后一半的确定,由于,（周期性）,所以,:,可,见,X(k),的后一半，也完全由,X,1,(k),X,2,(k),的前一半所确定。,*,N,点的,DFT,可由两个,N/2,点的,DFT,来计算。,又由于,所以,实现上式运算的流图称作蝶形运算,前一半,4.,蝶形运算,后一半,（,N/2,个蝶形,),(,前一半,),(,后一半,),1 1,1,1,-1,由,X,1,(k),、,X,2,(k),表示,X(k),的运算是一种特殊的运算,-,碟形运算,图,4.2.1,蝶形运算符号,x,1,(0)=,x,(0),x,1,(1)=,x,(2),N/2,点,x,1,(2)=,x,(4)DFT,x,1,(3)=,x,(6),x,2,(0)=,x,(1),x,2,(1)=,x,(3),N/2,点,x,2,(2)=,x,(5),DFT,x,2,(3)=,x,(7),X,1,(0),X,1,(1),X,1,(2),X,1,(3),X,2,(0),X,2,(1),X,2,(2),X,2,(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.,对,X,1,(k),和,X,2,(k),进行蝶形运算,前半部为,X(0)-X(3),后半部分为,X(4)-X(7),整个过程如下图所示,:,6.N/2,分解后的运算量：,复数乘法,复数加法,一个,N,/2,点,DFT,(,N,/2),2,N,/2(,N,/2,1),两个,N,/2,点,DFT,N,2,/2,N,(,N,/2,1),一个蝶形,1,2,N,/2,个蝶形,N,/2,N,总计,运算量减少了近一半,由于,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,时，即分解到,X,3,(,k,),，,X,4,(,k,),，,X,5,(,k,),，,X,6,(,k,),，,k,=0,1,不再做,DFT,(0),(4),(2),(6),(1),(5),(3),(7),W,N,0,W,N,0,W,N,0,W,0,N,-1,-1,-1,-1,X,(0),X,(1),X,(0),X,(1),X,(0),X,(1),X,(0),X,(1),3,3,4,4,5,5,6,6,W,N,0,W,N,2,W,N,0,W,N,2,-1,-1,-1,-1,X,(0),X,(1),X,(2),X,(3),X,(0),X,(1),X,(2),X,(3),1,1,1,2,1,2,2,2,W,W,W,W,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,=2,M,时，共有,M,级蝶形，每级,N,/2,个蝶形，每个蝶形有,1,次复数乘法，,2,次复数加法。,复数乘法：,复数加法：,与,DFT,比较,4.2.4 DIF-FFT,的运算规律 及编程思想,1,）原位计算,DIT-FFT,的运算过程很有规律，共进行,M,级运算，每级由,N/2,个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，与其它蝶形运算无关。这样,蝶形运算的两个输出值仍可放回蝶形运算的两个输入所在的存储器中,这种利用同一存储单元存储蝶形计算输入、输出的方法即为原位运算。每一级,(,列,),有,N/2,个蝶形运算,所以只需,N,个存储单元，可以节省存储单元。,运算规律,(0)=,A,0,(0),A,1,(0),A,2,(0)A,3,(0),=X(0),(4)=,A,0,(1),A,1,(1)A,2,(1)A,3,(1),=X(1),(2)=,A,0,(2),A,3,(2),=X(2),(6)=,A,0,(3),A,3,(3),=X(3),(1)=,A,0,(4),A,1,(4)A,2,(4)A,3,(4),=X(4),(5)=,A,0,(5),A,3,(5),=X(5),(3)=,A,0,(6),A,3,(6),=X(6),(7)=,A,0,(7),A,1,(7)A,2,(7)A,3,(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,个存储器中即为所求的频谱序列,X(k,),（,k=0,、,1,、,.,、,N-1,）。,这就是所谓的同址计算，这样可以大大节约存储元件。,N,点,DITFFT,运算流图中，每级都有,N/2,个蝶形。每个蝶形都要乘以因子,W,p,N,，称其为旋转因子，,p,称为旋转因子的指数。,2,）旋转引子的变化规律,3,）蝶形运算规律,对,N,=2,M,点,FFT,，,输入倒位序，输出自然顺序，,设第,L,级运算每个蝶形的两节点距离为,B,行，则,第,L,级运算：,旋转因子的确定仅与指数,P,有关，当,L,一定时,J,可以确定，进而可以确定指数,P,。对于同一,P,，参与本级蝶形运算的次数为 ，每级第一次调用 的蝶形结第一个输入节点的位置为,J,，第二次调用 的蝶形结第一个输入节点的位置为,J+2,L,，,即以,2,L,为步进，搜索下一蝶形运算第一个输入节点的位置。以此类推，直至做完 个蝶形运算。,同一级中，同一,P,的蝶形计算完成后，计算下一个,P,值对应的蝶形运算，直至完成本级所有蝶形运算。,DIT-FFT,原位计算步骤,1.,确定,L,；,2.,求出第,L,级的 个 ；,3.,对每一个 完成其所参与的 个蝶形运算，,4.,重复步骤,13,完成,M,级蝶形运算。,以,2,L,为步进参与本级同一旋转因子对应的其他蝶形运算,作为练习请写出,L=3,时的运算过程。,4,）编程思想及程序框图,在一般情况下，进行,FFT,运算的序列至少都有几百点的长度，因此需要编制,FFT,算法程序以便能够利用计算机来快速进行计算。输入倒位序，输出顺序的,DIT-FFT,的编程思想，,N,必须等于,2,的正整数幂，,FFT,的计算程序可以分为两部分：一部分是倒序重排，另一部分是用三层嵌套的循环来完成,M=log,2,N,次迭代。,三层循环的功能是：最里的一层循环完成相同,W,N,P,的蝶形运算，中间的一层循环完成因子,W,N,P,的变化，而最外的一层循环则是完成,M,次迭代过程。,循环,L,确定输入数据所在的位置,循环,J,L,、,P,、,B,确定后进行蝶形计算,5,）,序列的,倒序,由,DIT-FFT,的规律可知,输出,X,(k),按正常顺序排列在存储单元,而输入是按顺序,:,这种顺序称作倒位序,即二进制数,倒位。在编程时需完成倒位序，才能执,行原位计算。,n=0,0,n=1,0,n=0,1,n=1,1,n=0,1,n=1,1,0,1,0,1,0,1,0,1,(n,2,),x,(000)0,x,(100)4,x,(010)2,x,(110)6,x,(001)1,x,(101)5,x,(011)3,x,(111)7,（,偶）,（,奇）,倒位序由奇偶分组造成,以,N=8,为例,说明如下,:,0 0,0,0 0,0,0 0,1 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,倒序,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=2,M,，,M,是输入的正整数）,