第五章 快速傅里叶变换(蝶形运算)

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,第五章 快速傅里叶变换,本章目录,直接计算,DFT,的问题及改进的途径,按,时间抽取,的,基,2-FFT,算法,按,频率抽取,的基,2-FFT,算,法,快速,傅里叶逆变换,(IFFT),算法,Matlab,实现,2,5.1,引言,DFT,在实际应用中很重要,:,可以计算信号的频谱、功率谱和线性卷积等。,直接按,DFT,变换进行计算，当序列长度,N,很大时，计算量非常大，所需时间会很长。,FFT,并不是一种与,DFT,不同的变换，而是,DFT,的一种快速计算的算法。,3,5.2,直接计算,DFT,的问题及改进的途径,DFT,的运算量,设复序列,x,(,n,),长度为,N,点，其,DFT,为,k,=0,，,，,N,-1,（,1,）计算一个,X,(,k,),值的运算量,复数乘法次数：,N,复数加法次数：,N,1,4,5.2.1 DFT,的运算量,（,2,）计算全部,N,个,X,(,k,),值的运算量,复数乘法次数：,N,2,复数加法次数：,N,(,N,1),（,3,）对应的实数运算量,5,一次复数乘法,：,4,次实数乘法,2,次实数加法,一个,X,(,k,),：,4,N,次实数乘法,2,N+,2(,N,-1)= 2(2,N,-1),次实数加法,所以,整个,N,点,DFT,运算共需要：,N,2(2,N,-1)= 2,N,(2,N,-1),实数乘法次数：,4,N,2,实数加法次数：,6,DFT,运算量的结论,N,点,DFT,的复数乘法次数举例,N,N2,N,N2,2,4,64,4049,4,16,128,16384,8,64,256,65 536,16,256,512,262 144,32,1028,1024,1 048 576,结论,：当,N,很大时，其运算量很大，对实时性很强的信号处理来说，要求计算速度快，因此需要改进,DFT,的计算方法，以大大减少运算次数。,7,5.2.2,减少运算工作量的途径,主要原理是利用系数 的以下特性对,DFT,进行分解：,（,1,）对称性,（,2,）周期性,（,3,）可约性,另外，,8,5.3,按时间抽取的基,2-FFT,算法,算法原理,按时间抽取基,-2FFT,算法与直接计算,DFT,运算量的比较,按时间抽取的,FFT,算法的特点,按时间抽取,FFT,算法的其它形式流程图,9,5.3.1,算法原理,设,N,2,L,，将,x,(,n,),按,n,的奇偶分为两组：,r,=0,，,1,，,，,则,10,式中，,X,1,(,k,),和,X,2,(,k,),分别是,x,1,(,n,),和,x,2,(,n,),的,N,/2,的,DFT,。,另外，式中,k,的取值范围是：,0,，,1,，,，,N/2,1,。,11,因此， 只能计算出,X,(k,),的前一半值。,后一半,X,(k,),值，,N/2,，,N/2,1,，,，,N,？,利用,可得到,同理可得,12,考虑到,因此可得后半部分,X,(,k,),及前半部分,X,(,k,),k=,0,，,1,，,，,N/2,1,k,=0,，,1,，,，,N/2,1,13,蝶形运算,蝶形运算式,蝶形运算信号流图符号,因此，只要求出,2,个,N/2,点的,DFT,，即,X,1,(k),和,X,2,(k),，再经过蝶形运算就可求出全部,X,(k,),的值，运算量大大减少。,14,以,8,点为例第一次按奇偶分解,以,N=8,为例，,分解为,2,个,4,点的,DFT,，然后做,8/2=4,次蝶形运算即可求出所有,8,点,X,(k,),的值。,15,蝶形运算量比较,复数乘法次数：,N,2,复数加法次数：,N,(,N,1),复数乘法次数：,2*(,N/,2,),2,+,N,/2=,N,2,/2+,N,/2,复数加法次数：,2*(,N/2),(,N/2,1)+2*,N,/2=,N,2,/2,N,点,DFT,的运算量,分解一次后所需的运算量,2,个,N/2,的,DFT,N/2,蝶形：,因此通过一次分解后，运算工作量减少了差不多一半。,16,进一步按奇偶分解,由于,N,2,L,，因而,N,/2,仍是偶数 ，可以进一步把每个,N,/2,点,子序列再按其奇偶部分分解为两个,N,/4,点的子序列。,以,N,/2,点序列,x,1,(,r,),为例,则有,k,=0,1,17,且,k,=0,1,由此可见，一个,N,/2,点,DFT,可分解成两个,N,/4,点,DFT,。,同理，也可对,x,2,(n),进行同样的分解，求出,X,2,(k),。,18,以,8,点为例第二次按奇偶分解,19,算法原理,对此例,N,=8,，最后剩下的是,4,个,N,/4= 2,点的,DFT,，,2,点,DFT,也可以由蝶形运算来完成。以,X,3,(,k,),为例。,k,=0, 1,即,这说明，,N=2,M,的,DFT,可全部由蝶形运算来完成。,20,以,8,点为例第三次按奇偶分解,N=8,按时间抽取法,FFT,信号流图,21,5.3.2,按时间抽取基,2-FFT,算法与直接计算,DFT,运算量的比较,由按时间抽取法,FFT,的信号流图可知，当,N,=2,L,时，共有,级,蝶形运算；每级都由,个蝶形运算组成，而每个蝶形有,次复乘、,次复加，因此每级运算都需,次复乘和,次复加。,L,N,/2,N,/2,1,2,N,22,这样,级运算总共需要：,L,复数乘法：,复数加法：,直接,DFT,算法运算量,复数乘法：,复数加法：,N,2,N,(,N,1),直接计算,DFT,与,FFT,算法的计算量之比为,M,23,FFT,算法与直接,DFT,算法运算量的比较,N,N,2,计算量之比,M,N,N,2,计算量之比,M,2,4,1,4.0,128,16 384,448,36.6,4,16,4,4.0,256,65 536,1 024,64.0,8,64,12,5.4,512,262 144,2 304,113.8,16,256,32,8.0,1024,1 048 576,5 120,204.8,32,1028,80,12.8,2048,4 194 304,11 264,372.4,64,4049,192,21.4,24,5.3.3,按时间抽取的,FFT,算法的特点,序列的逆序排列,同址运算,（原位运算）,蝶形运算两节点间的距离,的确定,25,序列的逆序排列,由于,x,(,n,),被反复地按奇、偶分组，所以流图输,入端的,排列不再是顺序的，但仍有规律可循：,因为,N,=2,M,，,对于任意,n,（,0,n,N,-1),，,可以用,M,个二进制码表示为：,n,反复按奇、偶分解时，即按二进制码的“,0” “1”,分解。,序列的逆序排列,26,倒位序的树状图（,N=8,）,27,码位的倒位序,(N=8),自然顺序,n,二进制数,倒位序二进制数,倒位序顺序数,0,000,000,0,1,001,100,4,2,010,010,2,3,011,110,6,4,100,001,1,5,101,101,5,6,110,011,3,7,111,111,7,28,倒位序的变址处理（,N=8,）,29,同址运算,（原位运算）,某一列任何两个节点,k,和,j,的节点变量进行蝶形运算后，得到结果为下一列,k,、,j,两节点的节点变量，而和其他节点变量无关。这种原位运算结构可以节省存储单元，降低设备成本。,运算前,运算后,例,同址运算,（原位运算）,30,观察原位运算规律,31,蝶形运算两节点间的距离,以,N=8,为例：,第一级蝶形，距离为：,第二级蝶形，距离为：,第三级蝶形，距离为：,规律,：对于共,L,级的蝶形而言，其,m,级蝶形运算的节,点间的距离为,1,2,4,蝶形运算两节点间的距离,32,的确定,以,N=8,为例：,的确定,33,5.4,按频率抽取的基,2-FFT,算法,算法原理,再把输出,X,(,k,),按,k,的奇偶分组,先把输入按,n,的顺序分成前后两半,设序列长度为,N,=2,L,，,L,为整数,前半子序列,x,(,n,),后半子序列,0,n,0,n,34,5.4.1,算法原理,由,DFT,定义得,k,=0,，,1,，,，,N,35,由于,所以,则,k,=0,，,1,，,，,N,36,然后按,k,的奇偶可将,X,(,k,),分为两部分,r,=0,，,1,，,，,则式,可转化为,37,令,n,=0,，,1,，,，,代入,r,=0,，,1,，,，,可得,为,2,个,N,/2,点的,DFT,，合起来正好是,N,点,X,(,k,),的值。,38,蝶形运算,将,称为蝶形运算,与时间抽选基,2FFT,算法中的蝶形运算符号略有不同。,39,例 按频率抽取,(N=8),例 按频率抽取，将,N,点,DFT,分解为两个,N/2,点,DFT,的组合,(N=8),40,与时间抽取法的推导过程一样，由于,N,=2,L,，,N,/2,仍然是,一个偶数，因而可以将每个,N,/2,点,DFT,的输出再分解为偶数组,与奇数组，这就将,N,/2,点,DFT,进一步分解为两个,N,/4,点,DFT,。,N=8,41,5.4.2,频率抽取法与时间抽取法的异同,频率抽取法输入是自然顺序，输出是倒位序的；时间抽取法正好相反。,频率抽取法的基本蝶形与时间抽取法的基本蝶形有所不同。,频率抽取法运算量与时间抽取法相同。,频率抽取法与时间抽取法的基本蝶形是互为转置的。,42,5.5,快速傅里叶逆变换,(IFFT),算法,IDFT,公式,DFT,公式,比较可以看出，,IDFT,多出,M,个,1/2,可分解到,M,级蝶形运算中。,43,例 频率抽取,IFFT,流图,(N=8),44,快速傅里叶逆变换另一种算法,45,5.8,Matlab,实现,用,FFT,进行谱分析的,Matlab,实现,用,CZT,进行谱分析的,Matlab,实现,在,Matlab,中使用的线性调频,z,变换函数为,czt,，其调用格式为,X=,czt(x, M, W, A),其中，,x,是待变换的时域信号,x(n,),，其长度为,N,，,M,是变换的长度，,W,确定变换的步长，,A,确定变换的起点。若,M= N,，,A= 1,，则,CZT,变成,DFT,。,46,5.8.1,用,FFT,进行谱分析的,Matlab,实现,例,5.1,设模拟信号 ，以,t,= 0.01,n,(,n,=0:,N,-1),进行取样，试用,fft,函数对其做频谱分析。,N,分别为：,(1),N,=45,；,(2),N,=50,；,(3),N,=55,；,(2),N,=60,。,程序清单如下,%,计算,N,=45,的,FFT,并绘出其幅频曲线,N=45;n=0:N-1;t=0.01*n;,q=n*2*pi/N;,x=2*sin(4*pi*t)+5*cos(8*pi*t);,y=,fft(x,N,);,figure(1),subplot(2,2,1),plot(q,abs(y,),title(FFT,N=45),47,例,5.1,程序清单,%,计算,N,=50,的,FFT,并绘出其幅频曲线,N=50;n=0:N-1;t=0.01*n;,q=n*2*pi/N;,x=2*sin(4*pi*t)+5*cos(8*pi*t);,y=,fft(x,N,);,figure(1),subplot(2,2,2),plot(q,abs(y,),title(FFT,N=50),48,%,计算,N,=55,的,FFT,并绘出其幅频曲线,N=55;n=0:N-1;t=0.01*n;,q=n*2*pi/N;,x=2*sin(4*pi*t)+5*cos(8*pi*t);,y=,fft(x,N,);,figure(1),subplot(2,2,3),plot(q,abs(y,),title(FFT,N=55),49,%,计算,N,=60,的,FFT,并绘出其幅频曲线,N=60;n=0:N-1;t=0.01*n;,q=n*2*pi/N;,x=2*sin(4*pi*t)+5*cos(8*pi*t);,y=,fft(x,N,);,figure(1),subplot(2,2,4),plot(q,abs(y,),title(FFT,N=60),50,例,5.1,程序运行结果,从图中可以看出，这几种情况下均有较好的精度。,51,例,5.1,程序运行结果分析,分析：由,t,=0.01,n,进行取样可得，采样频率,fs,=100Hz,。而连续信号的最高模拟角频率为,8,，由,2,f,可得，最高频率为,8,/2,=4Hz,。因此，满足采样定理的要求。,采样序列为,即,为周期序列，周期,N,=50,。,将程序中,plot,改为,stem,函数，则可以更清楚地看出频谱。,52,例,5.1,修改程序运行结果,53,