第四章快速傅里叶变换 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(m)的DFT为 (k)=∑x(n)W,k=01…,N 考虑x(n)为复数序列的一般情况,对某一个k值,直接按(1)式计算X(k)值 需要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 次复数
乘法及NN1次复数加法运算。当N很大时,运算量相当可。 2、减少运算量的基本途径 因为(a+b)≥a2+b2,所以,把N点DFT分解为几个较短的DFT,可 是乘法次数大大减少。 旋转因子Wm具有明显的周期性和对称性,即 m+IN Wm+=eM =W Wm=WXN-m或者「WA-m=W 因此,可纵通过不断地把长序列的DFT分解成几个短序列的DFT,并利 W如的周期性和对称性来减少运算次数。 42.2时域抽取法基2FFT基本原理 FFT算法大致可以分为的域抽取法FFT(DIFT, Decimation-In- Time fft) 和颜域抽取法( DIF-FFT, Decimation- n-Frequency FFT)两大类。首先介绍 DIT-FFT算法。 1、序列分解 设序列x(n)的长度为N,且满足 N=2M,M为自然数 按n的奇偶把x(n)分解为两个点的子序列 x()=x(2).r=01-…N-1 x()2=x(2x+1,r=0.1…,1 则x(m)的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 为
(k)=∑x(m)+∑x(m) ∑x(2r)W+∑x(2r+1)x =∑x(r)形2+W∑x1()W 由于 所以 X(k)=X1(k)+X2(k),k=01,2…,N-1 其中x(k)和x(k)分别是x()和x:()的点DF,即 X()=∑x()w=DT[x(r) x2(k)=23()W*=DFT[x2(r)I ·由于x()和x()均以为周期,且=-形,所以X(k)又可表示为 X(k)=X1(k)+WX2(k),k=0, Xk+2=X{(4)-X2(k),k=0.1 这样,就将N点DFT分解为两个点的DFT和(7式以及(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
图421蝶形运算符号 采用该符号上述分解算法可表示为 x(1) M2点x(1) (3) 图422N点DFT的一次时域抽取分解图(N=8) 3、分解运算对提高运算效率的影响 由图42.1可见,要完成一次蝶形运算,需要一次复数乘和两次复数加运算。 由图422可见,经过一次分解后,计算一个N点DFT共需要计算两个一点DFT 和个蝶形运算。而计算一个点DFT需要 2丿次复数乘和(N 次复 数加运算。所以,按上述方法,N点DFT共需2 ≈(N>1时)次 复数乘法和N(2-1)+2=2次复数加法运算。由此可见,仅经过次分 就使运算量减少近一半 4、进一步处理 进一步,如果N=2m,还可以再分解。将x()按奇偶分解成两个长的子 序列x()和x(),即 -X N l=0,1, x x 那么,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(2)形2+x(2+)形 ∑x:()形4+W∑x() 2l= (k)+WX(k),0,1 N 式中 X()=∑x()=DF[x()] X()=∑x(H=DFT[x:(O] 同理,由X(k)和x4(k)的周期性和W的对称性WN4=-W,最后得到 X1(k)=X2(k)+WX4(k) N 0,1, 4 x(4)-2x() 用同样的方法可得 0.1 N Xs (k)-WNX(k) 其中 X(k)=∑x(0)=DF X()=∑x()w=DFT[x()] x5()=x2(2) x()=x1(21+1) 这样,经过第二次分解,又将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 式。依此类推
x(4)· DFT X() DFT XD) 4点 DFT 图4.2.3N点DFT的第二次时域抽取分解图(N=8) X(0) A(3) 图4.2.4N点DIT一FFT运算流图(N=8) 42.3 DIT-FFT算法与直接计算DFT运算量的比较 l、DI-FFT算法的运算量 由DI-FFT算法的分解过程可见,N=2"时,其运算流图应有M级蝶形, 每一级都由N2次复数城和N次复数加(每个蝶形需要两次复数加),所以,M 级运算总共需要的复数乘次数为 N 复数加次数为 C4(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 时
N2>log,从而,DFT算法比直接计算DFT的运算次数大大减少。 3、举例说明 例如:N=20=1024时, N210048567204.8 5120 直接计算 FFT算法 64128256512 N(取样点数 图4.2.5FFT算法与直接计算DFT所需乘法次数的比较曲线 4.2.4DII-FFT的运算规律及编程思想 1、原位计算 N=2M点的FFT共进行M级运算,每级由N2个蝶形运算组成。同一级中 每个蝶形的两个输入数据只对计算本蝶形有用,而且每个蝶形的输入、输出数据 节点又同在一条水平线上,这就意味着计算完一个蝶形后,所得的输出数据可立 即存入原输入数据所占的存储单元。这样,经过M级运算后,原来存放输入序 列数据的N个存储单元中便依次存放X(k)的N个值。这种利用同一存储单元存 储蝶形计算输入、输出数据的方法称为原位计算 原位计算可节省大量内存。 2、旋转因子的变化规律 N点DIFT运算流图中,每级都有N/2个蝶形。每个蝶形都要乘以因子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称为旋转因子的指嫩数。但各级的旋转因子和循环方法都有所 不同 旋转因子WP和运算级数的关系 用L表示从左到右的运算级数(L=1,2,,M)。仔细观察图4.24,可以发 现,对N=2的一般情况,第L级的旋转因子为 =Wy,J=0, 由于 所以 J=0.1 (12) P=J.24-L (13) 这样,就可按(12)式和(13)式确定第L级运算的旋转因子。 3、蝶形运算规律 用复数运算完成上述蝶形运算 设序列x(m)经时域抽选(倒序)后,存入数组ⅹ中。若蝶形运算的两个输 入数据相距B个点,应用原位计算,则蝶形运算可表示为: X()=X-()+X1(J+B)H X2(J+B)<X1()-X-(J+B)W 式中 p=J·2-;J=0,1,…,2--1,L=1,2,…,M 下标L表示第L级运算,x2()则表示第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(+B)WP=TR+jT X-(J)=X()+j() 式中下标R表示取实部,I表示取虚部 TR=XR(+B)cosp+X(J+B)sinp T=X,(J+B)cos p-XR(+B)sinp X2(J)=X(J)+(J) 则 XR()=XR()+TR x(J)=X(/)+T (+B)=X()-7R 4、编程思想及程序框图 第L级中,每个蝶形的两个输入数据相距B=2-个点; 同一旋转因子对应着间隔为2点的2-个蝶形 因此,先从输入端开始,逐级进行,共进行M级运算。在进行第L级运 算时,依次求出2-个不同的旋转因子,每求出一个旋转因子,就计算完 它对应的所有2M-个蝶形。 为了适应原位计算,其输入序列不是按x(m)的自然顺序排序,计算前应 先对序列进行倒序
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=2MJ k=J,N-1,2 x(k)∈X(k)+X(k+BW k(k+B)∈X(k)-X(k+B)W 结束 5、序列倒序 DIT-FFT算法的输出X(k)为自然顺序,但为了适应原位计算,其输入序列 不是按x(m)的自然顺序排序,这种经过M1次偶奇抽选后的排序称为序列x(m) 的倒序。因此,在运算之前应先对序列x(n)进行倒序。 0000 (n21h2 图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)
