四章快速傅立叶变换FastFourierTransformP课件

Click to edit Master title style,Click to edit Master text styles,Second Level,Third Level,Fourth Level,Fifth Level,*,第四章 快速傅立叶变换,F,ast,F,ourier,T,ransform,第一节 直接计算,DFT,的问题及改进途径,1、问题的提出,设有限长序列,x(n),，,非零值长度为,N,，,若对,x(n),进行一次,DFT,运算，共需,多大的运算工作量,？,计算成本?,计算速度?,2.,DFT,的运算量,回忆,DFT,和,IDFT,的变换式：,1）,x(n),为,复数,， 也为,复数,。,2）,DFT,与,IDFT,的,计算量相当。,注意：,计算机运算时（编程实现）：,N,次复乘，,N-1,次复加,N,个点,以,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),运算量,(,a+,jb,)(c+,jd,)=(ac-,bd,)+j(,bc,+ad),例：计算一个,N,点,DFT,，,共需,N,2,次复乘,。以做一次,复乘1,s,计，若,N,=4096，,所需时间为,例：石油勘探，有24个通道的记录，每通道波形记,录长度为5秒，若每秒抽样500点,/秒，,1）每道总抽样点数：500*5=2500点,2）24道总抽样点数：24*2500=6万点,3）DFT,复乘运算时间：,N,2,=(60000),2,=36*10,8,次,由于计算量大，且要求,相当大的内存,，,难以实现实时处理,，限制了,DFT,的应用。长期以来，人们一直在寻求一种能,提高,DFT,运算速度,的方法。,FFT,便是,Cooley &,Tukey,在1965 年提出的的快速算法，它可以使运算速度提高几百倍，从而使数字信号处理学科成为一个新兴的应用学科。,第二节 改善,DFT,运算效率的基本途径,1、利用,DFT,运算的系数 的固有对称性和周期,性，改善,DFT,的运算效率。,1）对称性,2）,周期性,3）可约性,2、将长序列,DFT,利用对称性和周期性分解为短,序列,DFT,的,思路,因为,DFT,的运算量与,N,2,成正比的，如果一个,大点数,N,的,DFT,能分解为,若干小点数,DFT,的组合,，则显然可以达到,减少运算工作量,的效果。,N,点,DFT,N/2,点,DFT,N/2,点,DFT,N/4,点,DFT,N/4,点,DFT,N/4,点,DFT,N/4,点,DFT,.,复乘：,FFT,算法的基本思想：,利用,DFT,系数的特性，合并,DFT,运算中的某些项,把长序列,DFT,短序列,DFT，,从而减少运算量。,FFT,算法分类:,时间抽选法,DIT: Decimation-In-Time,频率抽选法,DIF: Decimation-In-Frequency,第三节 按时间抽选的基2-,FFT,算法,1、算法原理,设输入序列长度为N=2,M,(M为正整数，将该序列,按时间顺序的奇偶分解,为越来越短的子序列，称为基2按时间抽取的FFT算法。也称为,Coolkey-Tukey,算法。,其中,基2,表示：N=2,M,，M为整数.若不满足这个条件，可以人为地加上若干零值（,加零补长,）使其达到 N=2,M,。,先将,x(n),按n的奇偶分为两组，作变量置换:,当n=偶数时，令n=2r;,当n=奇数时，令n=2r+1;,分组，变量置换,2、算法步骤,得到：,带入,DFT,中,所以,由于,?,X,1,(k)、X,2,(k),只有,N/2,个点，以,N/2,为周期；而,X,(k),却有,N,个点，以,N,为周期。要用,X,1,(k)、X,2,(k),表达全部的,X,(k),值，还必须利用,W,N,系数的,周期特性,。,后半部分,前半部分,又考虑到 的对称性：,有：,后半部分,前半部分,蝶形运算流图符号,说明：,(1) 左边两路为输入,(2) 右边两路为输出,(3) 中间以一个小圆表示加、,减运算（右上路为相加,输出、右下路为相减输,出),1个蝶形运算需要1次复乘，2次复加,复数乘法,复数加法,一个,N,点,DFT,N,2,N,(,N,1),一个,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,/2 + N/2,N,2,/2,N(N/2-1) + N,N,2,/2,运算量减少了近一半,分解后的运算量：,先将N=8点的DFT分解成2个4点DFT,：,可知：,时域上：,x(0),x(2),x(4),x(6),为偶子序列,x(1),x(3),x(5),x(7),为奇子序列,频域上：,X(0)X(3)，由X(k)给出,X(4)X(7)，由X(k+N/2)给出,例子：求,N=2,3,=8,点,FFT,变换 按,N=8N/2=4，,做4点的DFT,：,N=8,点的直接DFT的计算量为：,复乘：,N,2,次,=,64次,复加：,N(N-1)次,=,8,7,=56次,此外，还有4个蝶形结，每个蝶形结需要1次复乘，2次复加。,一共是：复乘4次，复加8次。,得到,X,1,(k),和,X,2,(k),需要：,复乘：(,N,/2),2,+,(,N,/2),2,次,= 32,次,复加：,N,/2,(N,/2,-1),+,N,/2,(N,/2,-1),=12+12,=,24,次,用分解的方法得到,X,(k),需要：,复乘：32+4 = 36,次,复加：24+8,=,32,次,N,点,DFT,的一次时域抽取分解图(,N=8),4,点,DFT,4,点,DFT,x(0),x(2),x(4),x(6),x(1),x(3),x(5),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),X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7),因为4点,DFT,还是比较麻烦，所以再继续分解。,若将,N/2(4,点)子序列按奇/偶分解成两个N/4点,(,2点,),子序列。即对将x,1,(r)和x,2,(r)分解成奇、偶两个N/4点,(,2点,),点的子序列。,那么，,X,1,(k),又可表示为,X,2,(k),也可以进行相同的分解：,注意：通常我们会把 写成 。,N,点,DFT,的第二次时域抽取分解图(,N=8),2点,DFT,2点,DFT,2点,DFT,2点,DFT,x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7),X,3,(0),X,3,(1),X,4,(0),X,4,(1),X,5,(0),X,5,(1),X,6,(0),X,6,(1),X,1,(0),X,1,(1),X,1,(2),X,1,(3),X,2,(0),X,2,(1),X,2,(2),X,2,(3),X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7),4,点,DFT,4,点,DFT,x(0),x(2),x(4),x(6),x(1),x(3),x(5),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),X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7),8,8,X,3,(0),X,3,(1),x(0)=x,3,(0),x(4)=x,3,(1),N,点,DITFFT,运算流图(,N=8),x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7),X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7),3、DITFFT,算法与直接计算,DFT,运算量的比较,1),、,N=2,M,的,DFT,运算可分成,M,级，每一级有,N/2,个蝶形,，每个蝶形有一次复乘两次复加。,2)、所以,M,级共有 次复乘和 次复加。,3)、若直接计算,DFT,，,需,N,2,次复乘和,N(N-1),次复加。,显然，当,N,较大时，有：,例如，,N=2,10,=1024,时,FFT,算法与直接计算,DFT,所需乘法次数的比较曲线,4、,DITFFT,的运算规律及编程思想,FFT,的每级（列）计算都是由,N,个复数数据（输入）两两构成一个蝶型（共,N/2,个蝶形）运算而得到另外,N,个复数数据（输出）。,当数据输入到存储器以后，每一组运算的结果，,仍然存放在这同一组存储器中,直到最后输出。,例：将,x(0),放在单元,A(0),中，将,x(4),放在单元,A(1),中，,W,8,0,放在一个暂存器中。,将,x(0),+,W,8,0,x(4),送回,A(0),单元,将,x(0),-,W,8,0,x(4),送回,A(1),单元,X,3,(0),X,3,(1),x(0),x(4),1）,原位运算 (亦称同址计算),x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7),X(0),X(1),X(2),X(3),X(4),X(5),X(6),X(7),回顾：,N,点,DITFFT,运算流图(,N=8),如上所述，,N,点,DITFFT,运算流图中，每级都有,N/2,个蝶形。每个蝶形都要乘以因子,W,N,P,，,称其为,旋转因子,，,p,称为旋转因子的指数。,2)旋转因子的变化规律,观察,FFT,运算流图发现，第,L,级共有2,L-1,个不同的旋转因子。,N=2,3,=8,时的各级旋转因子表示如下：,L=1,时，,W,N,p,=W,N/4,J,， N/4 =2,1,=2,L,， J=0,L=2,时，,W,N,p,=W,N/2,J,， N/2 =2,2,=2,L,， J=0，1,L=3,时，,W,N,p,=W,N,J,， N =2,3,=2,L,， J=0，1，2，3,对,N=2,M,的一般情况，第,L,级的旋转因子为：,设序列,x(n),经时域抽选(倒序)后，存入数组,X,中。如果蝶形运算的两个输入数据相距,B,个点(,B=2,L-1,)，,应用原位计算，则蝶形运算可表示成如下形式：,下标,L,表示第,L,级运算，,X,L,(J),则表示第,L,级运算后数组元素,X(J),的值。,3) 编程思想及流程图,开始,送入,x(n),和,N=2,M,调整输入,x(n),的顺序,for(L=1; L=M; L+),B = 2,L-1,for(J=0; J=B-1; J+),p = J2,M-L,for(k = J; k= N-1; k=k+2,L,),输出结果,结束,4）,码位倒序,由,N=8,蝶形图看出：原位计算时，FFT输出的X(k)的次序正好是顺序排列的，即X(0)X(7),但输入x(n)都不能按自然顺序存入到存储单元中，而是按x(0),x(4),x(2),x(6) ,x(,1,),x(,5,),x(,3,),x(,7,)的顺序存入存储单元,，,即为乱序输入，顺序输出。,这种顺序看起来相当杂乱，然而它是,有规律,的。即,码位倒读规则,。,自然顺序,n,二进制码表示,码位倒读,码位倒置顺序,n,以,N=8,为例：,0,1,2,3,4,5,6,7,000,001,010,011,100,101,110,111,000,100,010,110,001,101,011,111,0,4,2,6,1,5,3,7,看出：,码位倒读后,的顺序刚好是数据送入计算机内的顺序。,倒序规律,对于数,N，,在其二进制最高位加1，等于加,N/2。,若已知某个反序号为,J，,为求下一个反序号，可先判,J,的最高位：,1) 若为0，则把该位变成1（即加,N/2）,就得到下,一个反序号，,2) 若为1，则需判断次高位：, 若次高位为0，则把最高位变0（相当减去,N/2）,后，再把次高位变1（即加,N/4）。, 若次高位为1，则需判断次次高位,分析：,倒,序,排,列,算,法,的,流,程,图,正序序列已在数组,A ,中，输 入,N,LH= N/2 , j=LH , N1 = N-2,j=j-k,k=k/2,k=LH,jk,j=j+k,T = A(j),A(i) = A(j),A(j) = T,for(i=1; i=N1; i+),i j,N,Y,Y,N,第四节 按频率抽选的基2-,FFT,算法,在基2快速算法中，频域抽取法,FFT,也是一种常用的快速算法，简称,DIFFFT。,设序列,x(n),长度为,N=2,M,，,首先将,x(n),前后对半分开，得到两个子序列，其,DFT,可表示为如下形式,DIFFFT,一次分解运算流图(,N=8),4点,DFT,4点,DFT,x(0),x(1),x(2),x(3),x(4),x(5),x(6),x(7),X(0),X(2),X(4),X(6),X(1),X(3),X(5),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),DIFFFT,二次分解运算流图(,N=8),DIFFFT,运算流图(,N=8),时间抽取算法与频率抽取算法的比较,1) 频率抽选法和时间抽选法总的计算量是相同的,复乘：,复加：,2) 频率抽取法和时间抽取法一样，都适用于,原位运,算,， 即蝶形的输入和输出占用同一个存储单元。,3) 均存在码位倒序问题。,4) 频率抽选法和时间抽选法一样，基本运算也是蝶形,运算。但两者的蝶形形式略有不同。,第五节,IDFT,的快速算法,IFFT,上述,FFT,算法流图也可以用于离散傅里叶逆变换(,Inverse Discrete Fourier Transform，,简称,IDFT)。,比较,DFT,和,IDFT,的运算公式：,1) 旋转因子：,2) 系数：,DITIFFT,运算流图,DITIFFT,运算流图(防止溢出),如果希望直接调用,FFT,子程序计算,IFFT，,则可用下面的方法：,对上式两边同时取共轭，得：,例1、如果通用计算机的速度为平均每次复乘需要,5,s,，每次复加需要0.5,s，,用它来计算512点,的,DFTx(n)，,问：,1）直接计算需要多少时间？,2）用,FFT,需要多少时间？,解：1）用,DFT,进行运算：,复乘：,T,1,=,N,2,510,-6,=1.31072,秒,复加：,T,2,=,N(N-1),0.510,-6,=0.130816,秒,总共：,T=T,1,+T,2,=1.441536,秒,2）用,FFT,进行运算：,复乘：,T,1,=,(N/2)log,2,N,510,-6,=0.01152,秒,复加：,T,2,=,Nlog,2,N,0.510,-6,=0.002304,秒,总共：,T,=T,1,+T,2,=0.013824,秒,例2、对一个连续时间信号,x,a,(t),采样1秒得到4096个采,样点的序列，求：,1）若采样后没有发生混叠现象，,x,a,(t),的最高频率是,多少？,解：1秒内采样4096个点，说明采样频率是,4096,Hz,。,2）若计算采样信号的4096点,DFT，DFT,系数之间,的频率间隔是多少？,解：(,要求解的是频谱分辨的间隔,F,),例3、长度为240点的序列,x(n),与长度为,N,点的,h(n),卷,积。当,N=,10和240时，直接进行卷积,x(n)*h(n),和用,IFFTX(K)H(K),的方法相比，那种方法求,解,y(n),的效率更高？,x(n),h(n),y(n)=x(n)*h(n),LN,1,+N,2,-1,X(k),补,L,-,N,1,个零,x,(,n,),L,点,DFT,补,L,-,N,2,个零,h,(,n,),L,点,DFT,L,点,IDFT,y,(,n,),=,x,(,n,)*,h,(,n,),H(k),Y(k),直接进行卷积（,N=10）：,乘法次数：,24010,= 2400,次,用,FFT,的方法（,N=10）：,添零到256点，,L=256,乘法次数：,3,(L/2)log,2,L,L = 31288256 = 3328,次,直接进行卷积（,N=240）：,乘法次数：,240240,= 57600,次,用,FFT,的方法（,N=240）：,添零到512点，,L=512,乘法次数：,3,(L/2)log,2,L,L = 32569512 = 7424,次,