第五章 快速傅里叶变换(FFT)(数字信号处理).ppt

第五章 快速傅里叶变换(FFT),5.1 引言,DFT是信号分析与处理中的一种重要变换因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换(简称FFT)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的直到1965年发现了DFT的一种快速算法以后，情况才发生了根本的变化5.2 基2FFT算法,5.2.1 直接计算DFT的特点及减少运算量的基本途径 长度为N的有限长序列x(n)的DFT为 考虑x(n)为复数序列的一般情况，对某一个k值，直接按(5.2.1)式计算X(k)值需要N次复数乘法、(N-1)次复数加法5.2.1),如前所述，N点DFT的复乘次数等于N2显然，把N点DFT分解为几个较短的DFT，可使乘法次数大大减少另外，旋转因子WmN具有明显的周期性和对称性其周期性表现为,(5.2.2),其对称性表现为,或者,,5.2.2 时域抽取法基2FFT基本原理 FFT算法基本上分为两大类：时域抽取法FFT(Decimation In Time FFT,简称DIT-FFT)和频域抽取法FFT(Decimation In Frequency FFT,简称DIF―FFT)。

下面先介绍DIF―FFT算法 设序列x(n)的长度为N，且满足,为自然数,按n的奇偶性将x(n)分解为两个N/2点的子序列,则x(n)的DFT为,由于,所以,其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，即,(5.2.5),(5.2.6),由于X1(k)和X2(k)均以N/2为周期，且 ，所以X(k)又可表示为,(5.2.7),(5.2.8),图5.2.1 蝶形运算符号,图5.2.2 N点DFT的一次时域抽取分解图(N=8),与第一次分解相同，将x1(r)按奇偶分解成两个N/4长的子序列x3(l)和x4(l)，即,那么，X1(k)又可表示为,(5.2.9),式中,同理，由X3(k)和X4(k)的周期性和Wm N/2的对称性 Wk+N/4 N/2=-Wk N/2 最后得到：,(5.2.10),用同样的方法可计算出,(5.2.11),其中,,,图5.2.3 N点DFT的第二次时域抽取分解图(N=8),图5.2.4 N点DIT―FFT运算流图(N=8),5.2.3 DIT―FFT算法与直接计算DFT运算量的比较 每一级运算都需要N/2次复数乘和N次复数加(每个蝶形需要两次复数加法)。

所以，M级运算总共需要的复数乘次数为,复数加次数为,例如，N=210=1024时,图5.2.5 FFT算法与直接计算DFT所需乘法次数的比较曲线,5.2.4 DIT―FFT的运算规律及编程思想 1.原位计算 由图5.2.4可以看出，DIT―FFT的运算过程很有规律N=2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成 2.旋转因子的变化规律 如上所述，N点DIT―FFT运算流图中，每级都有N/2个蝶形每个蝶形都要乘以因子WpN，称其为旋转因子，p称为旋转因子的指数观察图5.2.4不难发现，第L级共有2 L-1个不同的旋转因子N=23=8时的各级旋转因子表示如下： L=1时，WpN=WJ N/4=WJ2L， J=0 L=2时， WpN =WJ N/2=WJ2L， J=0，1 L=3时， WpN =WJN=WJ2L， J=0，1，2，3 对N=2M的一般情况，第L级的旋转因子为,(5.2.12),(5.2.13),3. 蝶形运算规律 设序列x(n)经时域抽选(倒序)后，存入数组X中。

如果蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶形运算可表示成如下形式： X (J) XL-1(J)+X L-1(J+B)WpN XL(J+B) XL-1(J)-X L-1(J+B)WpN 式中 p=J·2 M-L；J=0,1,…，2 L-1-1；L=1,2,…，M,下标L表示第L级运算，XL(J)则表示第L级运算后数组元素X(J)的值如果要用实数运算完成上述蝶形运算，可按下面的算法进行 设 T=X L-1(J+B)WpN=TR+jTI X L-1(J)=X′R(J)+jX′I(J) 式中下标R表示取实部，I表示取虚部，,则,4. 编程思想及程序框图,图5.2.6 DIT―FFT运算和程序框图,5. 序列的倒序 DIT―FFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种倒序是很有规律的由于N=2M，所以顺序数可用M位二进制数(nM-1nM-2…n1n0)表示。

图5.2.7 形成倒序的树状图(N=23),表5.2.1 顺序和倒序二进制数对照表,图5.2.8 倒序规律,图5.2.9 倒序程序框图,5.2.5 频域抽取法FFT(DIF―FFT) 在基2快速算法中，频域抽取法FFT也是一种常用的快速算法，简称DIF―FFT 设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式：,偶数,奇数,将X(k)分解成偶数组与奇数组，当k取偶数(k=2r,r=0,1,…,N/2-1)时,(5.2.14),当k取奇数(k=2r+1,r=0,1,…,N/2-1)时,(5.2.15),将x1(n)和x2(n)分别代入(5.2.14)和(5.2.15)式，可得,(5.2.16),图5.2.10 DIF―FFT蝶形运算流图符号,图5.2.11 DIF―FFT一次分解运算流图(N=8),图5.2.12 DIF―FFT二次分解运算流图(N=8),图5.2.13 DIF―FFT运算流图(N=8),图5.2.14 DIT―FFT的一种变形运算流图,图5.2.15 DIT―FFT的一种变形运算流图,5.2.6 IDFT的高效算法 上述FFT算法流图也可以用于离散傅里叶逆变换(Inverse Discrete Fourier Transform，简称IDFT)。

比较DFT和IDFT的运算公式：,图5.2.16 DIT―IFFT运算流图,图5.2.17 DIT―IFFT运算流图(防止溢出),如果希望直接调用FFT子程序计算IFFT，则可用下面的方法： 由于,对上式两边同时取共轭，得,5.3 进一步减少运算量的措施,5.3.1 多类蝶形单元运算 由DIT―FFT运算流图已得出结论，N=2M点FFT共需要MN/2次复数乘法  由(5.2.12)式，当L=1时，只有一种旋转因子W0N=1，所以，第一级不需要乘法运算综上所述，先除去第一、二两级后，所需复数乘法次数应是 从L=3至L=M共减少复数乘法次数为,(5.3.1),(5.3.2),因此，　DIT―FFT的复乘次数降至,(5.3.3),从实数运算考虑，计算N=2M点DIT―FFT所需实数乘法次数为,(5.3.4),5.3.2 旋转因子的生成 在FFT运算中，旋转因子WmN=cos(2πm/N)-jsin(2πm/N)，求正弦和余弦函数值的计算量是很大的5.3.3 实序列的FFT算法 设x(n)为N点实序列，取x(n)的偶数点和奇数点分别作为新构造序列y(n)的实部和虚部，即,对y(n)进行N/2点FFT，输出Y(k)，则,根据DIT―FFT的思想及式(5.2.7)和(5.2.8)，可得到,由于x(n)为实序列，所以X(k)具有共轭对称性，X(k)的另外N/2点的值为,5.4 分裂基FFT算法,5.5.1 分裂基FFT算法原理 当n=pq，且p=N/4，q=4时，n可表示为,并有,再将上式中的k表示为,可得,对k0=0，1,2,3,并用k表示k1，用n表示n0，可以写出,(5.5.1),(5.5.2),(5.5.3),令,则(5.5.2)式可写成如下更简明的形式：,(5.4.4),图5.5.1 分裂基第一次分解L形流图,例如，N=16，第一次抽选分解时，由式(5.5.3)得 x2(n)=x(n)+x(n+8), 0≤n≤7x14(n)={［x(n)-x(n+8)］-j［x(n+4)-x(n+12)］}Wn16, 0≤n≤3x24(n)={［x(n)-x(n+8)］+j［x(n+4)-x(n+12)］}W3n16, 0≤n≤3 把上式代入式(5.4.4)，可得 X(2k)=DFT［x2(n)］, 0≤k≤7 X(4k+1)=DFT［x14(n)］, 0≤k≤3 X(4k+3)=DFT［x24(n)］, 0≤k≤3,,图5.5.2 分裂基FFT算法L形排列示意图与结构示意图(a)分裂基FFT算法L形排列示意图；(b)分裂基FFT算法运算流图结构示意图,图5.5.3 16点分裂基第一次分解L形流图(图中省去箭头),第二次分解： 先对图5.5.3中N/2点DFT进行分解。

令X1(l)=X(2l)，则有 X1(2l)=DFT［y2(n)］, 0≤l≤3 X1(4l+1)=DFT［y14(n)］, 0≤l≤1 X1(4l+3)=DFT［y24(n)］, 0≤l≤1,,其中y2(n)=x2(n)+x2(n+4), 0≤n≤3y14(n)={［x2(n)-x2(n+4)］-［x2(n+2)x(n+6)］}Wn8,n=0,1y24(n)={［x2(n)-x2(n+4)］+j［x2(n+2)x2(n+6)］}W3n8,n=0,1,图5.4.4 图5.4.4中N/2点DFT的分解L形流图,图5.5.5 4点分裂基L形运算流图,图5.4.6 16点分裂基FFT运算流图,5.5.2 分裂基FFT算法的运算量 设第j级有lj个L形，j=1,2,…,M-1,M=log2N，则有l1=N/4。

由图5.5.2(b)可见，第j-1列中的L形包含了第j列中的一部分结点的计算，即空白部分，所占结点数刚好等于第j-1列中所有L形对应结点的一半，所以第j列L形个数就减少l j-1/2个，即,。