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

1、第四章第四章快速傅里叶变换快速傅里叶变换（FFT）主要内容主要内容qDIT-FFT算法算法 qDIF-FFT算法算法qIFFT算法算法qChirp-FFT算法算法q线性卷积的线性卷积的FFT算法算法4.1 引言引言q FFT: Fast Fourier Transformq 1965年，年，Cooley-Turky 发表文章发表文章机器计算傅机器计算傅里叶级数的一种算法里叶级数的一种算法，提，提出出FFT算法，解决算法，解决DFT运算量太大，在实际使用中受限制的问题。运算量太大，在实际使用中受限制的问题。q FFT的应用。频谱分析、滤波器实现、实时信的应用。频谱分析、滤波器实现、实时信号处理等

2、。号处理等。q DSP芯片实现。芯片实现。TI公司的公司的TMS 320c30，10MHz时钟，基时钟，基2-FFT1024点点FFT时间时间15ms。q 典型应用：信号频谱计算、系统分析等典型应用：信号频谱计算、系统分析等)()(kXnxDFT )()()(nynhnxFFTnhnyIFFTFFTnx)()()( 系统分析系统分析 频谱分析与功率谱计算频谱分析与功率谱计算4.2 直接计算直接计算DFT的问题及改进途径的问题及改进途径10)()(NnknNWnxkX10)(1)(NkknNWkXNnx1、 DFT与与IDFT( )Nx n点有限长序列2、DFT与与IDFT运算特点运算特点复数乘

3、法复数乘法复数加法复数加法一个一个X(k)NN 1N个个X(k)(N点点DFT)N 2N (N 1)10( )NnkNnx n Wajbcjdacbdj adcb同理：同理：IDFT运算量与运算量与DFT相同。相同。实数乘法实数乘法实数加法实数加法一次复乘一次复乘42一次复加一次复加2一个一个X (k) 4N2N+2 (N 1)=2 (2N 1)N个个X (k)(N点点DFT)4N 22N (2N 1)3、降低、降低DFT运算量的考虑运算量的考虑nkNW 的特性*()() ()nknkN n kn N kNNNNWWWW对称性()() nkN n kn N kNNNWWW周期性 nkmnkNm

4、NWW可约性/nknk mNN mWW0/2(/2) 11NkNkNNNNWWWW 特殊点：2jnknkNNWeNknkNNWWnNnkNNWW2jmnkmNe221NjjNee FFT算法分类算法分类:q 时间抽选法时间抽选法DIT: Decimation-In-Timeq 频率抽选法频率抽选法DIF: Decimation-In-FrequencyFFTDFTDFTDFTDFT算法的基本思想： 利用系数的特性，合并运算中的某些项， 把长序列短序列，从而减少其运算量。4.3 按时间抽取（按时间抽取（DIT）的）的FFT算法算法12/.210) 12()()2()(21Nrrxrxrxrx，（

5、Decimation In Time）1、算法原理、算法原理设序列点数设序列点数 N = 2L，L 为整数。为整数。 若不满足，则补零若不满足，则补零将序列将序列x(n)按按n的奇偶分成两组：的奇偶分成两组：N为为2的整数幂的的整数幂的FFT算法称算法称基基-2FFT算法算法。将将N点点DFT定义式分解为两个长度为定义式分解为两个长度为N/2的的DFT10)()()(NnknNWnxnxDFTkXkrNnNrrkNnNrWrxWrx)12(12/0212/0) 12()2( 为奇为偶 )(12/02/2)(2/12/0121)()(kXNrrkNkNkXrkNNrWrxWWrx)()()(21

6、kXWkXkXkN记：记： （1 1）rkNrkNWW2/2（这一步利用：（这一步利用： ）,0,1,./2 1r kN再利用周期性求再利用周期性求X(k)的后半部分的后半部分/22NkNkkNNNNWWWW 又)(2)()()(222112/02/112/0)2/(2/11kXkNXkXWrxWrxkNXNrrkNNrkNrNrkNkNrNWW2/)2/(2/)2()2()2()2(12/,.2 , 1 , 0)()()(2)2/(121kNXWkNXkNXNkkXWkXkXkNNkN，12/,.2 , 1 , 0)()(21NkkXWkXkN，将上式表达的运算用一个专用将上式表达的运算用一

7、个专用“蝶形蝶形”信流图表示。信流图表示。)(1kX)(2kX)()(21kXWkXkN)()(21kXWkXkNkNW1212( )( )( )()( )( )2kNkNX kX kW XkNX kX kW Xk0,1,.,/2 1kN注：注：a. 上支路为加法，下支路为减法；上支路为加法，下支路为减法； b. 乘法运算的支路标箭头和系数。乘法运算的支路标箭头和系数。用用“蝶形结蝶形结”表示上面运算的分解：表示上面运算的分解： 328N)0(x)1 (x)2(x)3(x)4(x)5(x)6(x)7(x)0(X) 1 (X)2(X)3(X)4(X)5(X)6(X)7(X1NW0NW2NW3NW

8、)0(1X)1 (1X)2(1X)3(1X)0(2X)1 (2X(3)2X)2(2XDFTN点2DFTN点2分解后的运算量：分解后的运算量：复数乘法复数乘法复数加法复数加法一个一个N/2点点DFT(N/2)2N/2 (N/2 1)两个两个N/2点点DFTN2/2N (N/2 1)一个蝶形一个蝶形12N/2个蝶形个蝶形N/2N总计总计22/2/2/2NNN2/2 1/2N NNN运算量减少了近一半运算量减少了近一半进一步分解进一步分解MN2122MN2N4N由于由于 ， 仍为偶数，因此，两个仍为偶数，因此，两个 点点DFTDFT又可同样进一步分解为又可同样进一步分解为4 4个个 点的点的DFTD

9、FT。1314(2 )( )(21)( )xlx lxlx l0,1,.,/4 1lN13/2413/24( )( )( )()( )( )4kNkNX kXkWXkNX kXkWXk0,1,.,14Nk 02/NW12/NW)(3lx)(4lx)2(x)4(x)6(x)0(x)0(1X) 1 (1X)2(1X) 3(1X) 0(3X) 1 (3X)0(4X) 1 (4XDFTN点4DFTN点4“蝶形蝶形”信流图表示信流图表示 N点点DFT分解为四个分解为四个N/4点的点的DFTDFTN点4DFTN点4DFTN点4DFTN点4)2(x)4( x)6( x)0( x) 1 ( x) 3 ( x)

10、5(x)7( x0NW2NW0NW2NW1NW0NW2NW3NW)0(X) 1 (X) 2(X) 3(X) 4(X) 5(X)6(X)7(X)(.kX)(.nxq 类似的分解一直继续下去，直到分解为最类似的分解一直继续下去，直到分解为最后的两类蝶形运算为止后的两类蝶形运算为止(2点点DFT).q 如上述如上述N=8=23，N/4=2点中：点中： 类似进一步分解类似进一步分解1点点DFTx(0)1点点DFTx(4)X3(0)X3(1)02W进一步简化为蝶形流图：进一步简化为蝶形流图：0NWX3(0)X3(1)x(0)x(4)4()0()4()0()0(004/3xWxxWxXNN)4()0()4

11、()0() 1 (014/3xWxxWxXNN因此因此8 8点点FFTFFT时间抽取方法的信流图如下时间抽取方法的信流图如下)2(x)4(x)6(x)0(x) 1 ( x) 3 ( x)5(x)7(x0NW0NW0NW0NW第一级.0NW2NW0NW2NW 第二级.)(0kX1m)(1kX)(2kX)(3kX2m3m1NW0NW2NW3NW)0(X) 1 (X)2(X)3(X)4(X)5(X)6(X)7(X 第三级.FFT运算量与运算特点运算量与运算特点 1 N=2L时，共有时，共有L=log2N级运算；每一级有级运算；每一级有N/2个蝶形结。个蝶形结。2每一级有每一级有N个数据中间数据），且

12、每级只用到本个数据中间数据），且每级只用到本级的转入中间数据，适合于迭代运算。级的转入中间数据，适合于迭代运算。3计算量：计算量： 每级每级N/2次复乘法，次复乘法，N次复加。（每蝶形只乘一次，次复加。（每蝶形只乘一次，加减各一次）。共有加减各一次）。共有L*N/2=N/2log2N 次复乘法；次复乘法；复加法复加法L*N=Nlog2N 次。与直接次。与直接DFT定义式运算定义式运算量相比量相比(倍数倍数) N2/(Nlog2N) 。当。当 N大时，此倍数大时，此倍数很大。很大。222()2()loglog2FFmDFTNNNmFFTNN比较比较DFT 参考参考P150 表表4-1 图图4-6

13、可以直观看出，当点数可以直观看出，当点数N越大时，越大时，FFT的优点更突出。的优点更突出。按时间抽取按时间抽取FFT蝶形运算特点蝶形运算特点 1、关于、关于FFT运算的混序与顺序处理（位倒序处理）运算的混序与顺序处理（位倒序处理） 由于输入序列按时间序位的奇偶抽取，故输入序由于输入序列按时间序位的奇偶抽取，故输入序列是混序的，为此需要先进行混序处理。列是混序的，为此需要先进行混序处理。混序规律：混序规律： x(n)按按n位置进行码位（二进制）倒置位置进行码位（二进制）倒置规律输入，而非自然排序，即得到混序排列。所规律输入，而非自然排序，即得到混序排列。所以称为位倒序处理。以称为位倒序处理。位

14、倒序实现：位倒序实现：（1）DSP实现采用位倒序寻址实现采用位倒序寻址（2）通用计算机实现可以有两个方法：一是严格按）通用计算机实现可以有两个方法：一是严格按照位倒序含义进行；二是倒进位的加照位倒序含义进行；二是倒进位的加N/2。倒位序倒位序自然序自然序00000000100410010102201011063011001141001015510101136110111771112 102( )()x nnn n n倒位序倒位序例例 计算计算 ， 。计算。计算 点点FFTFFT。用时间抽取输入倒序算法，。用时间抽取输入倒序算法，问倒序前寄存器的数问倒序前寄存器的数 和倒序后和倒序后 的的数据值

15、？数据值？2)(nnx31.2 , 1 , 0n32N)13( x)13( x16913)13(2x5232 N210)01101()13(102)22()10110(48422)13(2x解：倒序前解：倒序前 倒序倒序 倒序为倒序为 倒序后倒序后 DIT FFT中最主要的蝶形运算实现中最主要的蝶形运算实现（1）参与蝶形运算的两类结点（信号）间）参与蝶形运算的两类结点（信号）间“距离距离”（码地址）与其所处的第几级蝶形有关；第（码地址）与其所处的第几级蝶形有关；第m级的级的“结距离结距离”为为 (即原位计算迭代）即原位计算迭代）（2）每级迭形结构为）每级迭形结构为，12m)2(,.2 , 1L

16、NLmrNmmmmmrNmmmmWkXkXkXWkXkXkX)2()()2()2()()(1111111q 蝶形运算两节点的第一个节点为蝶形运算两节点的第一个节点为k值，表值，表示成示成L位二进制数，左移位二进制数，左移L m位，把右位，把右边空出的位置补零，结果为边空出的位置补零，结果为r的二进制数。的二进制数。2( )2L mrk（3） 的确定：的确定： 第第m级的级的r取值：取值：rNWkNWDIT算法的其他形式流图算法的其他形式流图q 输入倒位序输出自然序输入倒位序输出自然序q 输入自然序输出倒位序输入自然序输出倒位序q 输入输出均自然序输入输出均自然序q 相同几何形状相同几何形状q

17、输入倒位序输出自然序输入倒位序输出自然序q 输入自然序输出倒位序输入自然序输出倒位序参考参考P154-155时间抽取、时间抽取、 输入自然顺序、输入自然顺序、 输出倒位序的输出倒位序的FFTFFT流图流图 0NW0NW0NW2NW1110NW10NW2NW1111X(0)x(0)X(4)x(1)X(2)x(2)X(6)x(3)X(1)x(4)X(5)x(5)X(3)x(6)X(7)x(7)110NW0NW2NW1NW3NW11 例例 用用FFT算法处理一幅算法处理一幅NN点的二维图像，如用每秒可点的二维图像，如用每秒可做做10万次复数乘法的计算机，当万次复数乘法的计算机，当N=1024时，问需

18、要多少时间时，问需要多少时间（不考虑加法运算时间）？（不考虑加法运算时间）？ 解解 当当N=1024点时，点时，FFT算法处理一幅二维图像所需复数乘算法处理一幅二维图像所需复数乘法约为法约为 次，仅为直接计算次，仅为直接计算DFT所需时间的所需时间的10万万分之一。分之一。 即原需要即原需要3000小时，现在只需要小时，现在只需要2 分钟。分钟。 722210log2NN4.4 按频率抽取（按频率抽取（DIF）的）的FFT算法算法q 与与DIT-FFT算法类似分解，但是抽取的是算法类似分解，但是抽取的是X(k)。即分解即分解X(k)成奇数与偶数序号的两个序列。成奇数与偶数序号的两个序列。q 设

19、：设： N = 2L，L 为整数。将为整数。将X(k)按按k的奇偶分的奇偶分组前，先将输入组前，先将输入x(n)按按n的顺序分成前后两半：的顺序分成前后两半：（Decimation In Frequency)一、算法原理一、算法原理12/12/0)()()(NNnnkNNnnkNWnxWnxkX12/0)(212/02)()(NnknNNNnnkNNWnxWnx12/022/)()(NnnkNNkNNWnxWnxkNkNW) 1(2/2 10( )( 1)2NknkNnNx nx nW 0,1,.,1kN下面讨论下面讨论：的)(12,2kXrrk12/02/212/022) 1 ()()()(

20、)()2(NnnrNNNnrnNNWnxnxWnxnxrX12/02/212/0)12(2)2()()()()() 12(NnnrNnNNNnnrNNWWnxnxWnxnxrX按按k k的奇偶将的奇偶将X(k)X(k)分成两部分：分成两部分：显然：显然：。点的对应两个DFTNrXrX2/) 12(),2(nNNNWnxnxnxnxnxnx)()()()()()(2221)()(2NnxnxnNNNWnxnxnxnx)()()()(22nNW令：令：用蝶型结构图表示为：用蝶型结构图表示为：x1(0)x1(1)-1x1(2)x1(3)-1x2(0)x2(1)-1x2(2)x2(3)-1N/2点DF

21、TN/2点DFTx(0)x(7)x(1)x(2)x(3)x(4)x(5)x(6)X1(0)=X(0)X2(0)=X(1)X1(1)=X(2)X1(2)=X(4)X1(3)=X(6)X2(1)=X(3)X2(2)=X(5)X2(3)=X(7)1NW0NW2NW3NW311411/2( )( )(/4)( ) ( )(/4)nNx nx nx nNx nx nx nNW0,1,.,14Nn 313414( )(2 ) ( )( )(21)( )X kXkDFT x nX kXkDFT x n0,1,.,14Nk N/2仍为偶数，进一步分解：仍为偶数，进一步分解：N/2 N/4x3(0)x3(1)-

22、1-1x4(0)x4(1)N/4点DFTN/4点DFTx1(0)x1(1)x1(2)x1(3)X3(0)=X1(0)=X(0)X4(0)=X1(1)=X(2)X3(1)=X1(2)=X(4)X4(1)=X1(3)=X(6)0/2NW1/2NWq 按照以上思路继续分解，即一个按照以上思路继续分解，即一个N/2的的DFT分解成两个分解成两个N/4点点DFT，直到只计算，直到只计算2点的点的DFT，这就是，这就是DIF-FFT算法。算法。2个个1点的点的DFT蝶形流图蝶形流图 进一步简化为蝶形流图：进一步简化为蝶形流图：1点点DFTx3(0)1点点DFTx3(1)X(0)X(4)02W02WX(0)

23、X(4)x3(0)x3(1)2(x)4(x)6(x)0(x)1(x)3(x)5(x)7(x)0(X)1(X)2(X)3(X)4(X)5(X)6(X)7(X0NW0NW1NW2NW3NW0NW0NW0NW2NW0NW2NW0NW1m第一级： 2m第二级：3m第三级：二、按频率抽取二、按频率抽取FFT蝶形运算特点蝶形运算特点1）原位计算）原位计算1111( )( )( )( )( )( )mmmrmmmNXkXkXjXjXkXj WrNW1( )mXk1( )mXj( )mXk( )mXj-1L级蝶形运算，每级级蝶形运算，每级N/2个蝶形，每个蝶形结构：个蝶形，每个蝶形结构：m表示第表示第m级迭代

24、，级迭代，k，j表示数据所在的行数表示数据所在的行数2）蝶形运算）蝶形运算对对N=2L点点FFT，输入自然序，输出倒位序，输入自然序，输出倒位序，两节点距离：两节点距离：2L-m=N / 2m1111( )( )2( )22mmmmrmmmNmmNXkXkXkNNXkXkXkW第第m级运算：级运算：q 蝶形运算两节点的第一个节点为蝶形运算两节点的第一个节点为k值，表示值，表示成成L位二进制数，左移位二进制数，左移m-1位，把右边空出的位，把右边空出的位置补零，结果为位置补零，结果为r的二进制数。的二进制数。rNW 的确定12( )2mrk存储单元存储单元输入序列输入序列x(n) : N个存储单

25、元个存储单元rNW系数系数 ：N / 2个存储单元个存储单元三、三、DIT与与DIF的异同的异同q 基本蝶形不同基本蝶形不同2log2FNmN2logFaNNq DIT: 先复乘后加减先复乘后加减q DIF: 先减后复乘先减后复乘q 运算量相同运算量相同q 都可原位运算都可原位运算q DIT和和DIF的基本蝶形互为转置的基本蝶形互为转置4.5 IDFT的的FFT算法算法（FFT应用一）应用一） 一、从定义比较分析一、从定义比较分析knNNkWkXNkXIDFTnx10)(1)()(10)()()(NnknNWnxnxDFTkX与与DFT的比较：的比较： 1）、旋转因子）、旋转因子WN-kn 的

26、不同；的不同； 2）、结果还要乘）、结果还要乘 1/N。 )(10*10*)(1)(1)()(kXDFTknNNkknNNkWkXNWkXNkXIDFTnx二、实现算法二、实现算法直接使用直接使用FFT程序的算法程序的算法*)(1)(kXDFTNnx共轭共轭FFT共轭共轭乘乘1/ N( )X k*( )Xk( )x n直接调用直接调用FFT子程序计算子程序计算IFFT的方法：的方法：4.6 线性调频线性调频Z变换（变换（Chirp-Z变换）算变换）算法法 （FFT应用二）应用二）单位圆与非单位圆采样单位圆与非单位圆采样（a） 沿单位圆采样沿单位圆采样; (b) 沿沿AB弧采样弧采样 ooooA

27、BX(ej)RezRezjImzjImzAB(a)(b)X(ej) 螺线采样螺线采样 0jImzRezo0A01.0zM1A0W01z0(M I)0zk=AW-k k=0, 1, , M-1 0000jjeWWeAAChirp-Z变换的线性系统表示变换的线性系统表示 x(n)2/ 2 )(nWnh2/ 2nnWAg(n)h(n)1 / h(n)X(zn)由于系统的单位脉冲响应由于系统的单位脉冲响应 可以想象为频率可以想象为频率随时间随时间(n)(n)呈线性增长的复指数序列。在雷达系统中呈线性增长的复指数序列。在雷达系统中, ,这这种信号称为种信号称为线性调频信号（线性调频信号（Chirp Si

28、gnalChirp Signal），因此，这，因此，这里的变换称为里的变换称为线性调频线性调频Z Z变换变换。22)(nWnh一、基本算法思路一、基本算法思路10)()()()()(MmmnxmhnhnxnyLMMd)1()(nMhnh2/LMmd4.7 线性卷积的线性卷积的FFT算法算法（FFT应用三）应用三）若若L点点x(n)，M点点h(n)，则直接计算其线性卷积则直接计算其线性卷积y(n)需运算量：需运算量：若系统满足线性相位，即：若系统满足线性相位，即：则需运算量：则需运算量：FFT法：以圆周卷积代替线性卷积法：以圆周卷积代替线性卷积21mNML令 ( )01( )01x nnLx n

29、LnN( )01( )01h nnMh nMnNN( )( )* ( )( ) ( )y nx nh nx nh n则 NN2log2NN2log2NN2log2)()() 1nhFFTkH)()()2nxFFTkX)()()()3kXkHkY)()()4kYIFFTnyN总运算量：总运算量： 次乘法次乘法NNNmF2log23比较直接计算和比较直接计算和FFT法计算的运算量法计算的运算量22(13/2*log)dmFmMLKmNN222413/2*(1log)106logmMMKMMM223logmMKL讨论：讨论：ML12NMLM 则1）当）当1NMLL 则2）当）当mLK 需采用分段卷积

30、重叠相加法重叠保留法ML x(n)长度很长时，将长度很长时，将x(n)分为分为L长的若干长的若干小的片段，小的片段，L与与M可比拟。可比拟。nLiniLnxnxi，其它，01) 1()()(iinxnx)()()()()(nhnxnyiinhnx)()(1 1、重叠相加法、重叠相加法iiny)( 则：则： 输出：输出：)()()(nhnxnyii1MLN其中：其中：可以用圆周卷积计算：可以用圆周卷积计算：MN2 选选 ，上面圆周卷积可用，上面圆周卷积可用FFTFFT计算。计算。 )()()(nhnxnyiiN 由于由于yi(n)长度为长度为N，而，而xi(n)长度长度L ，必有，必有M-1点重

31、叠，点重叠， yi(n)应相加才能构成最后应相加才能构成最后y(n)的。的。iinyny)()(h(n)0N 1M 1x(n)0L2L3LnnnnnL 10 x0(n)N 10 x1(n)L2L 1LN 13L 10 x2(n)2L2LN 1重叠相加法图形重叠相加法图形nnnN 10y0(n)x0(n) h (n)N2L2L N 100L N 1Ly1(n)x1(n) h (n)Ny2(n)x2(n) h (n)N和上面的讨论一样，和上面的讨论一样， 用用FFT法实现重叠相加法的步骤如下法实现重叠相加法的步骤如下: 计算计算N点点FFT， H(k)=DFTh(n)； 计算计算N点点FFT，Xi

32、(k)=DFTxi(n)； 相乘，相乘，Yi(k)=Xi(k)H(k)； 计算计算N点点IFFT，yi(n)=IDFTYi(k)； 将各段将各段yi(n)（包括重叠部分）相加，（包括重叠部分）相加， 。重叠相加的名称是由于各输出段的重叠部分相加而得名的。重叠相加的名称是由于各输出段的重叠部分相加而得名的。 )()(0nynyii例例 已知序列已知序列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解解: 重叠相加法重叠相加法x1n=2, 3, 4, 5, 6x2n=7, 8, 9, 10, 11x3n=12,13,