第四章 快速傅里叶变换(FFT)

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,第四章快速傅里叶变换（,FFT,）,主要内容,DIT-FFT,算法,DIF-FFT,算法,IFFT,算法,Chirp-FFT,算法,线性卷积的,FFT,算法,4.1,引言,FFT:,Fast Fourier Transform,1965,年，,Cooley-,Turky,发表文章,机器计算傅里叶级数的一种算法,，提,出,FFT,算法，解决,DFT,运算量太大，在实际使用中受限制的问题。,FFT,的应用。频谱分析、滤波器实现、实时信号处理等。,DSP,芯片实现。,TI,公司的,TMS 320c30,，,10MHz,时钟，基,2-FFT1024,点,FFT,时间,15ms,。,典型应用：信号频谱计算、系统分析等,系统分析,频谱分析与功率谱计算,4.2,直接计算,DFT,的问题及改进途径,1,、,DFT,与,IDFT,2,、,DFT,与,IDFT,运算特点,复数乘法,复数加法,一个,X,(,k,),N,N ,1,N,个,X,(,k,),(,N,点DFT),N,2,N,(,N ,1),同理：,IDFT,运算量与,DFT,相同。,实数乘法,实数加法,一次复乘,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),3,、降低,DFT,运算量的考虑,FFT,算法分类,:,时间抽选法,DIT: Decimation-In-Time,频率抽选法,DIF: Decimation-In-Frequency,4.3,按时间抽取（,DIT,）的,FFT,算法,（,D,ecimation,I,n,T,ime）,1,、算法原理,设序列点数,N = 2,L,，,L,为整数。,若不满足，则补零,将序列,x(n),按,n,的奇偶分成两组：,N,为,2,的整数幂的,FFT,算法称,基,-2FFT,算法,。,将,N,点,DFT,定义式分解为两个长度为,N/2,的,DFT,记：,（,1,）,（这一步利用： ）,再利用周期性求,X(k),的后半部分,将上式表达的运算用一个专用“蝶形”信流图表示。,注：,a.,上支路为加法，下支路为减法；,b.,乘法运算的支路标箭头和系数。,用“蝶形结”表示上面运算的分解：,分解后的运算量：,复数乘法,复数加法,一个,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,总计,运算量减少了近一半,进一步分解,由于 ， 仍为偶数，因此，两个 点,DFT,又可同样进一步分解为,4,个 点的,DFT,。,“蝶形”信流图表示,N,点,DFT,分解为四个,N/4,点的,DFT,类似的分解一直继续下去，直到分解为最后的两类蝶形运算为止,(,2,点,DFT).,如上述,N=8=2,3,，,N/4=2,点中：,类似进一步分解,1,点,DFT,x(0),1,点,DFT,x(4),X,3,(0),X,3,(1),进一步简化为蝶形流图：,X,3,(0),X,3,(1),x(0),x(4),因此,8,点,FFT,时间抽取方法的信流图如下,FFT,运算量与运算特点,1,N=2,L,时，共有,L=log,2,N,级运算；每一级有,N/2,个蝶形结。,2,每一级有,N,个数据（中间数据），且每级只用到本级的转入中间数据，适合于迭代运算。,3,计算量：,每级,N/2,次复乘法，,N,次复加。（每蝶形只乘一次，加减各一次）。共有,L*N/2=N/2log,2,N,次复乘法；复加法,L*N=Nlog,2,N,次。与直接,DFT,定义式运算量相比,(,倍数,) N,2,/(Nlog,2,N),。,当,N,大时，此倍数很大。,比较,DFT,参考,P150,表,4-1,图,4-6,可以直观看出，当点数,N,越大时，,FFT,的优点更突出。,按时间抽取,FFT,蝶形运算特点,1,、关于,FFT,运算的混序与顺序处理（位倒序处理）,由于输入序列按时间序位的奇偶抽取，故输入序列是混序的，为此需要先进行混序处理。,混序规律：,x(n),按,n,位置进行码位（二进制）倒置规律输入，而非自然排序，即得到混序排列。所以称为位倒序处理。,位倒序实现：,（,1,）,DSP,实现采用位倒序寻址,（,2,）通用计算机实现可以有两个方法：一是严格按照位倒序含义进行；二是倒进位的加,N/2,。,倒位序,自然序,000,0,0,000,100,4,1,001,010,2,2,010,110,6,3,011,001,1,4,100,101,5,5,101,011,3,6,110,111,7,7,111,倒位序,例计算 ， 。计算 点,FFT,。,用时间抽取输入倒序算法，问倒序前寄存器的数 和倒序后 的数据值？,解：倒序前,倒序,倒序为,倒序后,DIT FFT,中最主要的蝶形运算实现,（,1,）参与蝶形运算的两类结点（信号）间“距离”（码地址）与其所处的第几级蝶形有关；第,m,级的“结距离”为,(,即原位计算迭代）,（,2,）每级迭形结构为,蝶形运算两节点的第一个节点为,k,值，表示成,L,位二进制数，左移,L m,位，把右边空出的位置补零，结果为,r,的二进制数。,（,3,） 的确定：,第,m,级的,r,取值：,DIT,算法的其他形式流图,输入倒位序输出自然序,输入自然序输出倒位序,输入输出均自然序,相同几何形状,输入倒位序输出自然序,输入自然序输出倒位序,参考,P154-155,时间抽取、 输入自然顺序、 输出倒位序的,FFT,流图,例 用,FFT,算法处理一幅,N,N,点的二维图像，如用每秒可做,10,万次复数乘法的计算机，当,N,=1024,时，问需要多少时间（不考虑加法运算时间）？,解 当,N,=1024,点时，,FFT,算法处理一幅二维图像所需复数乘法约为 次，仅为直接计算,DFT,所需时间的,10,万分之一。 即原需要,3000,小时，现在只需要,2,分钟。,4.4,按频率抽取（,DIF,）的,FFT,算法,与,DIT-FFT,算法类似分解，但是抽取的是,X(k),。,即分解,X(k),成奇数与偶数序号的两个序列。,设：,N = 2,L,，,L,为整数。将,X(k),按,k,的奇偶分组前，先将输入,x(n),按,n,的顺序分成前后两半：,（Decimation In Frequency),一、算法原理,下面讨论,按,k,的奇偶将,X(k),分成两部分：,显然：,令：,用蝶型结构图表示为：,x,1,(0),x,1,(1),-1,x,1,(2),x,1,(3),-1,x,2,(0),x,2,(1),-1,x,2,(2),x,2,(3),-1,N,/2,点,DFT,N,/2,点,DFT,x,(0),x,(7),x,(1),x,(2),x,(3),x,(4),x,(5),x,(6),X,1,(0)=,X,(0),X,2,(0)=,X,(1),X,1,(1)=,X,(2),X,1,(2)=,X,(4),X,1,(3)=,X,(6),X,2,(1)=,X,(3),X,2,(2)=,X,(5),X,2,(3)=,X,(7),N/2,仍为偶数，进一步分解：,N/2,N/4,x,3,(0),x,3,(1),-1,-1,x,4,(0),x,4,(1),N,/4,点,DFT,N,/4,点,DFT,x,1,(0),x,1,(1),x,1,(2),x,1,(3),X,3,(0)=,X,1,(0)=,X,(0),X,4,(0)=,X,1,(1)=,X,(2),X,3,(1)=,X,1,(2)=,X,(4),X,4,(1)=,X,1,(3)=,X,(6),按照以上思路继续分解，即一个,N/2,的,DFT,分解成两个,N/4,点,DFT,，,直到只计算,2,点的,DFT,，,这就是,DIF-FFT,算法。,2,个,1,点的,DFT,蝶形流图,进一步简化为蝶形流图：,1,点,DFT,x,3,(0),1,点,DFT,x,3,(1),X(0),X(4),X(0),X(4),x,3,(0),x,3,(1),二、按频率抽取,FFT,蝶形运算特点,1,）原位计算,-1,L,级蝶形运算，每级,N/2,个蝶形，每个蝶形结构：,m,表示第,m,级迭代，,k,，,j,表示数据所在的行数,2,）蝶形运算,对,N=2,L,点,FFT,，,输入自然序，输出倒位序，,两节点距离：,2,L-m,=N / 2,m,第,m,级运算：,蝶形运算两节点的第一个节点为,k,值，表示成,L,位二进制数，左移,m,-1,位，把右边空出的位置补零，结果为,r,的二进制数。,存储单元,输入序列,x(n) : N,个存储单元,系数 ：,N / 2,个存储单元,三、,DIT,与,DIF,的异同,基本蝶形不同,DIT,:,先复乘后加减,DIF,:,先减后复乘,运算量相同,都可原位运算,DIT,和,DIF,的基本蝶形互为转置,4.5 IDFT,的,FFT,算法（,FFT,应用一）,一、从定义比较分析,与,DFT,的比较：,1,）、旋转因子,W,N,-,kn,的不同；,2,）、结果还要乘,1/N,。,二、实现算法,直接使用,FFT,程序的算法,共轭,FFT,共轭,乘1/,N,直接调用,FFT,子程序计算,IFFT,的方法：,4.6,线性调频,Z,变换（,Chirp-Z,变换）算法,（,FFT,应用二）,单位圆与非单位圆采样,（,a,）,沿单位圆采样,; (b),沿,AB,弧采样,螺线采样,z,k,=,AW,-k,k,=0, 1, ,M,-1,Chirp-Z,变换的线性系统表示,由于系统的单位脉冲响应 可以想象为频率随时间,(n),呈线性增长的复指数序列。在雷达系统中,这种信号称为,线性调频信号,（,Chirp Signal,）,，,因此，这里的变换称为,线性调频,Z,变换,。,一、基本算法思路,4.7,线性卷积的,FFT,算法（,FFT,应用三,）,若,L,点,x(n),，,M,点,h(n),，,则直接计算其线性卷积,y(n),需运算量：,若系统满足线性相位，即：,则需运算量：,FFT,法：以圆周卷积代替线性卷积,N,总运算量： 次乘法,比较直接计算和,FFT,法计算的运算量,讨论：,1,）当,2,）当,见,P183,表,x(n),长度很长时，将,x(n),分为,L,长的若干小的片段，,L,与,M,可比拟。,1,、重叠相加法,则：,输出：,其中：,可以,用圆周卷积计算：,选 ，上面圆周卷积可用,FFT,计算。,N,由于,y,i,(n,),长度为,N,，而,x,i,(n,),长度,L,，,必有,M-1,点重叠，,y,i,(n,),应相加才能构成最后,y(n),的。,重叠相加法图形,和上面的讨论一样， 用,FFT,法实现重叠相加法的步骤如下,:,计算,N,点,FFT,，,H,(,k,)=DFT,h,(,n,),；,计算,N,点,FFT,，,X,i,(,k,)=DFT,x,i,(,n,),；,相乘，,Y,i,(,k,)=,X,i,(,k,),H,(,k,),；,计算,N,点,IFFT,，,y,i,(,n,)=IDFT,Y,i,(,k,),；,将各段,y,i,(,n,),（,包括重叠部分）相加， ,。,重叠相加的名称是由于各输出段的重叠部分相加而得名的。,例 已知序列,xn=n+2,0,n,12, hn=1,2,1,试,利用重叠相加法计算线性卷积,取,L=5,。,yn=2, 7, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 41, 14,解,:,重叠相加法,x,1,n=2, 3, 4, 5, 6,x,2,n=7, 8, 9, 10, 11,x,3,n=12,13, 14,y,1,n=2, 7, 12, 16, 20, 17, 6,y,2,n= 7, 22, 32, 36, 40, 32, 11,y,3,n=12, 37, 52, 41, 14,2,、重叠保留法,此方法与上述方法稍有不同。,先将,x,(,n,),分段，每段,L,=,N,-,M,+1,个点，这是相同的。,不同之处是，序列中补零处不补零，而在每一段的前边补上前一段保留下来的（,M,-1,）,个输入序列值， 组成,L,+,M,-1,点序列,x,i,(,n,),。,如果,L,+,M,-12,m,，,则可在每段序列末端补零值点，补到长度为,2,m,，,这时如果用,DFT,实现,h,(,n,),和,x,i,(,n,),圆周卷积，则其每段圆周卷积结果的前（,M,-1,）,个点的值不等于线性卷积值，必须舍去。,重叠保留法示意图,重叠保留法示意图,y,k,=2, 7, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 41, 14,x,1,k,=,0, 0,2, 3, 4,x,2,k,=,3, 4,5, 6 ,7,x,3,k,=,6 ,7 ,8, 9, 10,y,1,k,=,x,1,k,h,k,= ,11, 4,2, 7, 12,x,4,k,=,9, 10 ,11, 12,13,y,2,k,=,x,2,k,h,k,= ,23, 17, 16, 20, 24,y,3,k,=,x,3,k,h,k,= ,35, 29, 28, 32, 36,y,4,k,=,x,4,k,h,k,= ,47, 41, 40, 44, 48,x,5,k,=,12,13,14, 0, 0,y,5,k,=,x,5,k,h,k,= ,12, 37, 52, 41, 14,解,:,重叠保留法,例 已知序列,xn=n+2,0,n,12, hn=1,2,1,试,利用重叠,保留法,计算线性卷积,取,L=5,。,语音信号消噪过程,信号淹没在啸叫噪声中；,(,b,),信号与噪声的功率谱；,(,c,),去噪后的功率谱；,(,d,),重构原语音信号,FFT,应用举例,