第四章快速傅里叶变换 41引言 DFT的计算量与N的平方成正比 DFT是信号分析与处理中的一种重要变换。但因直接计算DFT的计算量与 变换区间长度N的平方成正比,当N较大时,计算量很大,所以在快速傅里叶 变换出现以前,直接用DFT算法进行谱分析和信号实时处理是不切实际的。 直到1965年发现了DFT的一种快速算法,情况才发生了根本变化。 快速傅立叶变换发展 1965年,图基(JW.Tuky)和库力(Tw. Coody)在 Math. Computation 杂志发表著名的“机器计算傅立叶级数的一种算法”一文。桑德图 基快速算法出现,经改进,形成一套高效运算方法,这就是快速傅里 叶变换。 ◇1984年,杜哈梅尔和霍尔曼提出分裂基快速算法,使运算效率进 步提高 42基2FFT算法 -本章主要讨论基2FFT算法及其编程思想 42.1直接计算DFT的特点及减少运算量的基本途径 1、按DFT定义计算DFT的计算量 根据前面知识,长度为N的有限长序列x(n)的DFT为 x(k)=∑x(m)W,k=0,1…,N-1 考虑x(n)为复数序列的一般情况,对某一个k值,直接按(1)式计算X()值 需要N次复数乘法、(N-1)次复数加法。因此刷所有N个k值,共需N2次复数
第四章 快速傅里叶变换 4.1 引言 ⚫ DFT 的计算量与 N 的平方成正比 DFT 是信号分析与处理中的一种重要变换。但因直接计算 DFT 的计算量与 变换区间长度 N 的平方成正比，当 N 较大时，计算量很大，所以在快速傅里叶 变换出现以前，直接用 DFT 算法进行谱分析和信号实时处理是不切实际的。 ⚫ 直到 1965 年发现了 DFT 的一种快速算法，情况才发生了根本变化。 ⚫ 快速傅立叶变换发展 1965 年，图基（J. W. Tuky）和库力（T.W. Coody）在 Math.Computation 杂志发表著名的“机器计算傅立叶级数的一种算法”一文。桑德-图 基快速算法出现，经改进，形成一套高效运算方法，这就是快速傅里 叶变换。 1984 年，杜哈梅尔和霍尔曼提出分裂基快速算法，使运算效率进一 步提高。 4.2 基 2FFT 算法 ——本章主要讨论基 2FFT 算法及其编程思想 4.2.1 直接计算 DFT 的特点及减少运算量的基本途径 1、按 DFT 定义计算 DFT 的计算量 根据前面知识，长度为 N 的有限长序列 x n( ) 的 DFT 为 ( ) ( ) 1 0 , 0,1, , 1 N kn N n X k x n W k N − = = = − (1) 考虑 x n( ) 为复数序列的一般情况，对某一个 k 值，直接按(1)式计算 X k( ) 值 需要 N 次复数乘法、（N-1）次复数加法。因此对所有 N 个 k 值，共需 2 N 次复数
乘法及N(N1执复数加法运算。当N很大时,运算量相当可。 2、减少运算量的基本途径 因为(a+b)≥a2+b2,所以,把N点DFT分解为几个较短的DFT,可 是乘法次数大大减少。 旋转因子W具有明显的周期性和对称性,即 rm+/N 27(m+) 或者[-丁 WN2=-WN 因此,可纵通过不断地把长序列的DFT分解成几个短序列的DFT,并利用 W的期性和对称性水减少运算次数 422时域抽取法基2FFT基本原理 FFT算法大致可以分为时域抽歌法FFT( DIT-FFT, Decimation-In- Time fft) 和频域抽取法( DIF-FFT, Decimation-In-Frequency FFT)两大类。首先介绍 DIT-FFT算法。 1、序列分解 设序列x(n)的长度为N,且满足 N=2M,M为自然数 按n的奇偶把x(n)分解为两个点的子序列 x()=x(2),r=0.1.N x2(r)=x(2r+ 则x(n)的DFT为
乘法及 N(N-1)次复数加法运算。当 N 很大时，运算量相当可观。 2、减少运算量的基本途径 ⚫ 因为 ( ) 2 2 2 a b a b + + ，所以，把 N 点 DFT 分解为几个较短的 DFT，可 是乘法次数大大减少。 ⚫ 旋转因子 m WN 具有明显的周期性和对称性，即 ( ) 2 2 j m lN j m m lN m N N W e e W N N − + − + = = = (2) * m N m N m m W W W W N N N N − − − = = 或者 (3) 2 N m m W W N N + = − 因此，可以通过不断地把长序列的 DFT 分解成几个短序列的 DFT，并利用 kn WN 的周期性和对称性来减少运算次数。 4.2.2 时域抽取法基 2FFT 基本原理 FFT 算法大致可以分为时域抽取法 FFT（DIT-FFT, Decimation-In-Time FFT） 和频域抽取法（DIF-FFT, Decimation-In-Frequency FFT）两大类。首先介绍 DIT-FFT 算法。 1、序列分解 设序列 x n( ) 的长度为 N，且满足 2 , M N M = 为自然数 ⚫ 按 n 的奇偶把 x n( ) 分解为两个 2 N 点的子序列 1 ( ) (2 , 0,1, , 1 ) 2 N x r x r r = = − 2 ( ) (2 1 , 0,1, , 1 ) 2 N x r x r r = + = − 则 x n( ) 的 DFT 为
X(k ∑x(2r)形2+∑x(2r+1)H ∑x1()W2"+W∑x2(r) 由于 所以 X(k)=X(k)+WX2(k),k=0L12,…N-1 其中x1(k)和x2(k)分别是x()和x2()的点DFT,即 X1(k)=∑x()=DF[x( X(k)=∑x(r)形=DFT[x2(T) 由于x()和x2(k)均以2为周期,且W2=-形,所以X(k)又可表示为 X(k)=X1(k)+WX2(k),k=0 2 x1(k)-WX2(k),k=0,1 这样,就将N点DFT分解为两个N点的DFT和(式以及(8)式的运算 2、蝶形运算 为了形象直观地描述(7)和(8)式的运算,引入以下符号
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2 2 2 2 1 0 0 1 1 2 2 2 2 1 2 0 0 2 2 1 kn kn N N n n N N kr kr N N r r N N kr k kr N N N r r X k x n W x n W x r W x r W x r W W x r W = = − − + = = − − = = = + = + + = + 偶数 奇数 由于 2 2 2 2 2 2 j kr N j kr kr kr N W e e W N N − − = = = 所以 ( ) 1 2 ( ) ( ), 0,1,2, , 1 k X k X k W X k k N = + = − N (4) 其中 X k 1 ( ) 和 X k 2 ( ) 分别是 x r 1 ( ) 和 x r 2 ( ) 的 2 N 点 DFT，即 ( ) ( ) ( ) 1 2 1 1 1 0 2 N kr N r X k x r W DFT x r − = = = (5) ( ) ( ) ( ) 1 2 2 2 2 0 2 N kr N r X k x r W DFT x r − = = = (6) ⚫ 由于 X k 1 ( ) 和 X k 2 ( ) 均以 2 N 为周期，且 2 N k k W W N N + = − ，所以 X k( ) 又可表示为 ( ) 1 2 ( ) ( ), 0,1, , 1 2 k N N X k X k W X k k = + = − (7) 1 2 ( ) ( ), 0,1, , 1 2 2 k N N N X k X k W X k k + = − = − (8) 这样，就将 N 点 DFT 分解为两个 2 N 点的 DFT 和(7)式以及(8)式的运算。 2、蝶形运算 为了形象直观地描述(7)和(8)式的运算，引入以下符号 C A B A＋ BC A－BC
图42.1蝶形运算符号 采用该符号上述分解算法可表示为 M点x(1 x1(2) x(1) M2点X DFT 图422N点DFT的一次时域抽取分解图(N=8) 3、分解运算对提高运算效率的影响 由图42.1可见,要完成一次蝶形运算,需要一次复数乘和两次复数加运算。 由图422可见,经过一次分解后,计算一个N点DFT共需要计算两个点DFT 和个蝶形运算。而计算一个,点DT需要[N)次复数和2(2-次复 N 数加运算。所以,按上述方法,N点DFT共需2/N)2NN2(N>1时)次 复数垂注和(N12NN2次复数加法运算。由此可见,仅经过一次分 就做运量减少近一半 进一步处理 进一步,如果N=2",还可以再分解。将x(r)按奇偶分解成两个长的子 序列x:()和x4(),即 x()=x(21+)=0、N入 x()=x1(21) 那么,X(k)又可以表示为
图 4.2.1 蝶形运算符号 采用该符号上述分解算法可表示为 N/2点 DFT W N 0 N/2点 DFT W N 1 W N 2 W N 3 x(0) X1 (0) x(2) x(4) x(6) x(1) x(3) x(5) x(7) X1 (1) X1 (2) X1 (3) X2 (0) X2 (1) X2 (2) X2 (3) X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7) 图 4.2.2 N 点 DFT 的一次时域抽取分解图(N=8) 3、分解运算对提高运算效率的影响 由图 4.2.1 可见，要完成一次蝶形运算，需要一次复数乘和两次复数加运算。 由图 4.2.2 可见，经过一次分解后，计算一个 N 点 DFT 共需要计算两个 2 N 点 DFT 和 2 N 个蝶形运算。而计算一个 2 N 点 DFT 需要 2 2 N 次复数乘和 1 2 2 N N − 次复 数加运算。所以，按上述方法，N 点 DFT 共需 2 2 2 2 2 2 N N N + （N>>1 时）次 复数乘法和 2 1 2 2 2 2 N N N N − + = 次复数加法运算。由此可见，仅经过一次分解， 就使运算量减少近一半。 4、进一步处理 进一步，如果 2 m N = ，还可以再分解。将 x r 1 ( ) 按奇偶分解成两个 4 N 长的子 序列 x l 3 ( ) 和 x l 4 ( ) ，即 ( ) ( ) ( ) ( ) 3 1 4 1 2 , 0,1, , 1 2 1 4 x l x l N l x l x l = = − = + 那么， X k 1 ( ) 又可以表示为
X(k)=∑x1(21)+∑x1(2 x()WM4+W∑x1()WM 2 X3(k)+WX4(k),0, N X3(k)=∑x() x2(k)=x:()=DF[x( 同理,由X(4)和x()的周期性和W的对称丝么Nm,最后得到 2 x1(k)=X3(k)+WX4(k) N N k=0,1, 3 (k)-WNX,(k) 用同样的方法可得 2(k)=X5(k)+WX6(k) N k Y(k-WkX 其中 N xs (k)=E*s(0W*=DFT 3 (7 X(4)=2x(O)=DFr[x(O)] ()=x(2) x(2)=x(21+.1=01" 这样,经过第二次分解,又将N2点DFT分解为两个N4点DFT式。依此类推
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 4 4 2 2 1 1 1 / 2 1 / 2 0 0 1 1 4 4 3 / 4 4 / 4 0 0 2 3 4 2 2 2 1 , 0,1, , 1 2 N N M k l N N l l N N kl k kl N N N l l k N X k x l W x l W x l W W x l W N X k W X k − − + = = − − = = = + + = + = + − (9) 式中 ( ) ( ) ( ) 1 4 3 3 3 0 4 N kl N l X k x l W DFT x l − = = = ( ) ( ) ( ) 1 4 4 4 4 0 4 N kl N l X k x l W DFT x l − = = = 同理，由 X k 3 ( ) 和 X k 4 ( ) 的周期性和 2 m WN 的对称性 4 2 2 N k k W W N N + = − ，最后得到 ( ) ( ) ( ) ( ) ( ) 1 3 4 2 1 3 4 2 , 0,1, , 1 4 4 k N k N X k X k W X k N k N X k X k W X k = + = − + = − (10) 用同样的方法可得 ( ) ( ) ( ) ( ) ( ) 2 5 6 2 2 5 6 2 , 0,1, , 1 4 4 k N k N X k X k W X k N k N X k X k W X k = + = − + = − (11) 其中 ( ) ( ) ( ) 1 4 5 5 5 0 4 N kl N l X k x l W DFT x l − = = = ( ) ( ) ( ) 1 4 6 6 6 0 4 N kl N l X k x l W DFT x l − = = = ( ) ( ) ( ) ( ) 5 2 6 2 2 , 0,1, , 1 2 1 4 x l x l N l x l x l = = − = + 这样，经过第二次分解，又将 N/2 点 DFT 分解为两个 N/4 点 DFT 式。依此类推
DFTLX,(1) X,(1) DFT 图4.2.3N点DFT的第二次时域抽取分解图(N=8) A(0) xX(1) X(2) X(5) 图4.2.4N点DI一FT运算流图(N=8) 423 DIT-FFT算法与直接计算DFT运算量的比较 1、 DIT-FFT算法的运算量 由 DIT-FFT算法的分解过程可见,N=2M时,其运算流图应有M级蝶形, 每一级都由N2次复数城和N次复数加(每个蝶形需要两次复数加),所以,M 级运算总共需要的复数乘次数为 N C(2 复数加次数为 CA(2)=NM=Nlog? 2、根据定义直接计算的运算量 直接计算N点DFT需N2次复数乘,N(N-1)次复数加。当N>>1时
N/4点 DFT W N 1 2 W N 1 2 W N 0 W N 1 W N 2 W N 3 X1 (0) X1 (1) X1 (2) X1 (3) X2 (0) X2 (1) X2 (2) X2 (3) X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7) x(0) X3 (0) X3 (1) X4 (0) X4 (1) x(4) x(2) x(6) x(1) x(5) x(3) x(7) N/4点 DFT N/4点 DFT N/4点 DFT W N 0 2 W N 0 2 图 4.2.3 N 点 DFT 的第二次时域抽取分解图(N=8) W N 0 W N 1 W N 2 W N 3 W N 0 W N 2 W N 0 W N 2 W N 0 W N 0 W N 0 W N 0 x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7) A(0) A(1) A(2) A(3) A(4) A(5) A(6) A(7) A(0) A(1) A(2) A(3) A(4) A(5) A(6) A(7) A(0) A(7) X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7) A(0) A(1) A(2) A(3) A(4) A(5) A(6) A(7) 图 4.2.4 N 点 DIT―FFT 运算流图(N=8) 4.2.3 DIT-FFT 算法与直接计算 DFT 运算量的比较 1、DIT-FFT 算法的运算量 由 DIT-FFT 算法的分解过程可见， 2 M N = 时，其运算流图应有 M 级蝶形， 每一级都由 N/2 次复数城和 N 次复数加（每个蝶形需要两次复数加），所以，M 级运算总共需要的复数乘次数为 ( ) 2 2 log 2 2 N M N N C M = = 复数加次数为 ( ) 2 2 logN C N M N A = = 2、根据定义直接计算的运算量 直接计算 N 点 DFT 需 2 N 次复数乘， N N( −1) 次复数加。当 N>>1 时
log2,从而, DIT-FFT算法比直接计算DFT的运算次数大大减少 3、举例说明 例如:N=20=1024时, 10048567 =204.8 5120 直接计算 512 FHT算法 612825651 N(取样点数) 图4.2.5FFT算法与直接计算DFT所需乘法次数的比较曲线 424 DIT-FFT的运算规律及编程思想 1、原位计算 N=2M点的FFT共进行M级运算,每级由N2个蝶形运算组成。同一级中, 每个蝶形的两个输入数据只对计算本蝶形有用,而且每个蝶形的输入、输出数据 节点又同在一条水平线上,这就意味着计算完一个蝶形后,所得的输出数据可立 即存入原输入数据所占的存储单元。这样,经过M级运算后,原来存放输入序 列数据的N个存储单元中便依次存放X(k)的N个值。这种利用同一存储单元存 储蝶形计算输入、输出数据的方法称为原位计算 原位计算可节省大量内存。 旋转因子的变化规律 N点 DIT-FFT运算流图中,每级都有N2个蝶形。每个蝶形都要乘以因子W
2 2 log 2 n N N ，从而，DIT-FFT 算法比直接计算 DFT 的运算次数大大减少。 3、举例说明 例如： 10 N = = 2 1024 时， 2 2 10048567 204.8 5120 log 2 N N N = = 图 4.2.5 FFT 算法与直接计算 DFT 所需乘法次数的比较曲线 4.2.4 DIT-FFT 的运算规律及编程思想 1、原位计算 2 M N = 点的 FFT 共进行 M 级运算，每级由 N/2 个蝶形运算组成。同一级中， 每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据 节点又同在一条水平线上，这就意味着计算完一个蝶形后，所得的输出数据可立 即存入原输入数据所占的存储单元。这样，经过 M 级运算后，原来存放输入序 列数据的 N 个存储单元中便依次存放 X k( ) 的 N 个值。这种利用同一存储单元存 储蝶形计算输入、输出数据的方法称为原位计算。 原位计算可节省大量内存。 2、旋转因子的变化规律 N 点 DIT-FFT运算流图中，每级都有 N/2 个蝶形。每个蝶形都要乘以因子 p WN
称其为旋因子,p称为旋因孑的指数。但各级的旋转因子和循环方法都有所 不同 旋转因子W和运算级数的关系 用L表示从左到右的运算级数(L=1,2,,M)。仔细观察图424,可以发 现,对N=2M的一般情况,第L级的旋转因子为 0,1,2,…,24 由于 N- 2 所以 WR=H2-=W2,J=0,2,…21-1 (12) J·2 这样,就可按(12)式和(13)式确定第L级运算的旋转因子。 3、蝶形运算规律 ●用复数运算完成上述蝶形运算 设序列x(m)经时域抽选(倒序)后,存入数组ⅹ中。若蝶形运算的两个输 入数据相距B个点,应用原位计算,则蝶形运算可表示为: x2(J)<x-(J)+X-1(J+B)W X1(J+B)cX-l()-Xi+B)W 式中 p=J2-;J=0.1…,24--1,L=1,2…,M 下标L表示第L级运算,X()则表示第L级运算后数组元素X(J)的值 ●用实数运算完成上述蝶形运算
称其为旋转因子，p 称为旋转因子的指数。但各级的旋转因子和循环方法都有所 不同。 ⚫ 旋转因子 p WN 和运算级数的关系 用 L 表示从左到右的运算级数（L=1,2,…，M）。仔细观察图 4.2.4，可以发 现，对 2 M N = 的一般情况，第 L 级的旋转因子为 1 2 L , 0,1,2, ,2 1 p J L W W J N − = = − 由于 2 2 2 2 L M L M L M N − − = = 所以 2 1 2 , 0,1,2, ,2 1 M L L M p J J L W W W J N N N − − − = = = − (12) 2 M L p J − = (13) 这样，就可按(12)式和(13)式确定第 L 级运算的旋转因子。 3、蝶形运算规律 ⚫ 用复数运算完成上述蝶形运算 设序列 x n( ) 经时域抽选（倒序）后，存入数组 X 中。若蝶形运算的两个输 入数据相距 B 个点，应用原位计算，则蝶形运算可表示为： ( ) 1 1 ( ) ( ) p X J X J X J B W L L L N + + − − ( ) 1 1 ( ) ( ) p X J B X J X J B W L L L N + − + − − 式中 1 2 ; 0,1, ,2 1; 1,2, , M L L p J J L M − − = = − = 下标 L 表示第 L 级运算， X J L ( ) 则表示第 L 级运算后数组元素 X J( ) 的值。 ⚫ 用实数运算完成上述蝶形运算 设
T=XL(J+B)W=TR+jT, X-()=X(J)+jX(J) 式中下标R表示取实部,I表示取虚部, TR=XR(J+)cos p+X,(+B)si T,=X'(J+B)cosp-XR(J+B)sin X2()=XA(J)+jx/() 则 XR()=XR()+TR x1()=X/()+T X(+B=XR()-TR (+B)=X(J/)-T 4、编程思想及程序框图 第L级中,每个蝶形的两个输入数据相距B=24-个点; 同一旋转因子对应着间隔为2点的2M-个蝶形 因此,先从输入端开始,逐级进行,共进行M级运算。在进行第L级运 算时,依次求出24-个不同的旋转因子,每求出一个旋转因子,就计算完 它对应的所有2M-个蝶形。 为了适应原位计算,其输入序列不是按x(n)的自然顺序排序,计算前应 先对序列进行倒序
1 ( ) p T X J B W T jT = + = + L N R I − X J X J jX J L R I −1 ( ) = + ( ) ( ) 式中下标 R 表示取实部，I 表示取虚部， ( ) ( ) 2 2 cos sin T X J B p X J B p R R I N N = + + + ( ) ( ) 2 2 cos sin T X J B p X J B p I I R N N = + − + X J X J jX J L R I ( ) = + ( ) ( ) 则 X J X J T R R R ( ) = + ( ) X J X J T I I I ( ) = + ( ) X J B X J T R R R ( + = − ) ( ) X J B X J T I I I ( + = − ) ( ) 4、编程思想及程序框图 ⚫ 第 L 级中，每个蝶形的两个输入数据相距 1 2 L B − = 个点； ⚫ 同一旋转因子对应着间隔为 2 L 点的 2 M L− 个蝶形 ⚫ 因此，先从输入端开始，逐级进行，共进行 M 级运算。在进行第 L 级运 算时，依次求出 1 2 L− 个不同的旋转因子，每求出一个旋转因子，就计算完 它对应的所有 2 M L− 个蝶形。 ⚫ 为了适应原位计算，其输入序列不是按 x n( ) 的自然顺序排序，计算前应 先对序列进行倒序
开始 送入x(n),M L=1,M P=2- k=J,N-1,2 (k)∈X(k)+X(k+B)W x(k+B)∈X(k)-X(k+B)W 结束 5、序列倒序 DIT-FFT算法的输出X(k)为自然顺序,但为了适应原位计算,其输入序列 不是按x(n)的自然顺序排序,这种经过M1次偶奇抽选后的排序称为序列x(n) 的倒序。因此,在运算之前应先对序列x(n)进行倒序 11 图4.2.7形成倒序的树状图(N=23)
开 始 送入x(n)，M N＝ 2 M 倒 序 L＝ 1 , M Ｊ＝ 0 , B－ 1 P＝ 2 M －LJ k＝ J , N－1 , 2L p N p N X k B X k X k B W X k X k X k B W ( ) ( ) ( ) ( ) ( ) ( ) + − + + + 输 出 结 束 B 2 L－1 5、序列倒序 DIT-FFT 算法的输出 X k( ) 为自然顺序，但为了适应原位计算，其输入序列 不是按 x n( ) 的自然顺序排序，这种经过 M-1 次偶奇抽选后的排序称为序列 x n( ) 的倒序。因此，在运算之前应先对序列 x n( ) 进行倒序。 0 1 0 1 0 1 0 1 0 1 0 1 0 1 (n2 n1 n0 ) 2 000 0 4 2 6 1 5 3 7 100 010 110 001 101 011 111 图 4.2.7 形成倒序的树状图(N=23)