第3章快速付里叶变换FFT

1、第第3 3章章 快速付里叶变换快速付里叶变换(FFT)(FFT) Fast FourierFast Fourier TransformingTransformingu离散傅里叶变换在实际应用中是非常重要的，利用它可以计算信号的频谱、功率谱和线性卷积等.u但是，如果使用定义式来直接计算DFT，当N很大时，DFT的计算量太大，即使使用高速计算机，所花的时间也太多，很难对问题进行实时处理。所以在相当长的时间里，并没有得到真正的运用。因此，如何提高计算DFT的速度，便成了重要的研究课题。快速付里叶变换快速付里叶变换(FFT)(FFT) 1965年库利（J.W.Cooley）和图基（J.W.Tukey）

2、首次提出了DFT运算的一种快速算法，后来又有桑德（G.Sande）和图基（J.W.Tukey）的快速算法相继出现以后，情况才发生了根本的变化。人们开始认识到DFT运算的一些内在规律，从而很快地发展和完善了一套高速有效的运算方法，这就是现在人们普遍称之为快速傅里叶变换（FFT）的算法。 快速付里叶变换快速付里叶变换(FFT)(FFT)FFT的出现，使计算DFT的计算量可缩短一、二个数量级，还有效地减少了计算所需的存储容量，从而成为数字信号处理强有力的工具。FFT技术的应用极大地推动了DSP的理论和技术的发展，使DFT的运算在实际中真正得到了广泛的应用。uFFTFFT:各种各样快速计算DFT的方法

3、，统称为快速傅里叶变换(Fast Fourier Transform)，简称为FFT。快速付里叶变换快速付里叶变换(FFT)(FFT) FFTFFT:快速傅里叶变换（Fast Fourier Transform，简称为FFT）并不是一种新的变换， 而是离散傅里叶变换（DFT）的一种快速算法。由于有限长序列在其频域也可离散化为有限长序列（DFT），因此离散傅里叶变换（DFT）在数字信号处理中是非常有用的。快速付里叶变换快速付里叶变换(FFT)(FFT) 3.1 直接计算直接计算DFTDFT的特点及减少运的特点及减少运 算量的基本途径算量的基本途径 10NknNnX kx n W 101NknNk

4、x nXk WNDFTDFT及及IDFTIDFT的定义的定义0,1,1kN0,1,1nN将将DFTDFT定义式展开成方程组定义式展开成方程组：将方程组写成矩阵形式将方程组写成矩阵形式用向量表示用向量表示： 3.1 直接计算直接计算DFTDFT的特点及减少运的特点及减少运 算量的基本途径算量的基本途径 3.1 直接计算直接计算DFTDFT的特点及减少运的特点及减少运 算量的基本途径算量的基本途径 3.1 直接计算直接计算DFTDFT的特点及减少运的特点及减少运 算量的基本途径算量的基本途径 FFTFFT算法是基于可以将一个长度为算法是基于可以将一个长度为N N的序的序列的离散傅里叶变换逐次分解为

5、较短的列的离散傅里叶变换逐次分解为较短的离散傅里叶变换来计算这一基本原理的。离散傅里叶变换来计算这一基本原理的。这一原理产生了许多不同的算法，但它这一原理产生了许多不同的算法，但它们在计算速度上均取得了大致相当的改们在计算速度上均取得了大致相当的改善。善。3.1 3.1 直接计算直接计算DFTDFT的特点及减少运的特点及减少运 算量的基本途径算量的基本途径 3.1 直接计算直接计算DFTDFT的特点及减少运的特点及减少运 算量的基本途径算量的基本途径 nkNnkNWW*)(（2 2） W WN Nnknk的周期性的周期性 )()(NknNkNnNnkNWWW（3 3） W WN Nnknk的可

6、约性的可约性 mnkmNnkNnmkmNnkNWWWW/,kNNkNWW)2/(（1 1） W WN Nnknk的对称性的对称性 的特性的特性:nkNWnkNknNNkNnNWWW)()(10NW12/NNWjWNN4/ 3.1 直接计算直接计算DFTDFT的特点及减少运的特点及减少运 算量的基本途径算量的基本途径 的特性的特性:nkNWFFT DIT:按时间抽取法 （Decimation-in-Time） DIF:按频率抽取法 （Decimation-in-Frequency）。 按时间抽取（DIT）的基-2 FFT算法的基本出发点是，利用旋转因子WN nk的对称性和周期性，将一个长序列的D

7、FT分解为一些逐次变小的DFT来计算。分解过程遵循两条规则： 对时间进行奇偶分解； 对频率进行前后分解。3 3.2.2 按时间抽取（按时间抽取（DITDIT）的基）的基-2 FFT-2 FFT 算法算法原理原理（库利-图基算法）3 3.2.2 按时间抽取（按时间抽取（DITDIT）的基）的基-2 FFT-2 FFT 算法算法原理原理（库利-图基算法） 设序列x(n)的长度为N，且M为自然数N-point DFTMN22logMN 10,0,1,.,1NknNnX kx n WkN将其一分为二，分成偶数和奇数序列项（the even-indexed and odd-indexed terms）则

8、N/2点的序列为： even: x1(r)=x(2r), r=0, 1, , N/2-1 odd: x2(r)=x(2r+1) , r=0, 1, , N/2-13 3.2.2 按时间抽取（按时间抽取（DITDIT）的基）的基-2 FFT-2 FFT 算法算法原理原理（库利-图基算法）偶数序列偶数序列 the range: 02rN-2 (N is even) 0r(N/2)-1奇数序列奇数序列 the range: 12r+1N-1 (N is even) 02rN-2 0r(N/2)-13 3.2.2 按时间抽取（按时间抽取（DITDIT）的基）的基-2 FFT-2 FFT 算法算法原理原

9、理（库利-图基算法）3 3.2.2 按时间抽取（按时间抽取（DITDIT）的基）的基-2 FFT-2 FFT 算法算法原理原理（库利-图基算法） 则则x(n)的的DFT：rkNNrkNrkNNrkrNNrrkNNrNnnnkNNnNnnnkNnkNWrxWWrxWrxWrxWnxWnxWnxnxDFTkX)()( )() 12()2()()()()()(2120221201) 12(1202120101010为奇数为偶数故故: :2/2/2222NNjNjNWeeWx(n)的的DFT：由于由于:)()()()()()( )()(212/21202/11202120221201kXWkXWrxW

10、WrxWrxWWrxkXkNrkNNrkNrkNNrrkNNrkNrkNNrN/2-point DFTrkNNrrkNNrrkNNrrkNNrWrxWrxkXWrxWrxkX2/1202/212022/1202/11201) 12()()()2()()()()()(21kXWkXkXkNx x( (n n) )的前的前N N/2-point DFT/2-point DFT：k=0, 1, , N/2-1222221NkXWNkXNkXNkNx x( (n n) )的后的后N N/2-point DFT/2-point DFT：rkNNkrNWW2/22/)()()(212/120122/120

11、11kXWrxWrxkNXrkNNrkNrNNr W WN Nnknk的周期性的周期性 )()(NknNkNnNnkNWWWx(n)的的DFT：)(222kXkNX0,1,12Nk kNkNNNkNNWWWW2/2)()(221kXWkXNkXkN222221NkXWNkXNkXNkNx x( (n n) )的后的后N N/2-point DFT/2-point DFT：0,1,12Nk )()(221kXWkXNkXkN)()()(21kXWkXkXkNx x( (n n) )的的N N-point DFT-point DFT：表示上述算法可用蝶形结（表示上述算法可用蝶形结（ butterf

12、ly butterfly） 蝶形运算符号作图要素：作图要素：(1)左边两路为输入(2)右边两路为输出(3)中间以一个小圆表示加、减运算（右上路为相加输出、右下路为相减输出)(4)如果在某一支路上信号需要进行相乘运算，则在该支路上标以箭头，将相乘的系数标在箭头旁。(5)当支路上没有箭头及系数时，则该支路的传输比为1。0,1,12Nk )()(221kXWkXNkXkN)()()(21kXWkXkXkNx x( (n n) )的的N N-point DFT-point DFT：运算量？运算量？22N24N24N24N22N2N运算量？运算量？)()(221kXWkXNkXkN)()()(21kXW

13、kXkXkN运算量？运算量？x(n)x(n)y(n)y(n)数字信数字信号处理号处理？ 例：直接计算直接计算DFT【例例：】 计算有限长序列, N=4的DFT X(k)。【解】 由于x(n)是4点序列，因此X(k)也是4点序列。 ( )2, 4, 3, 4x n 340,03nknX kx n Wk 2j44ejW 所以，若将所以，若将k k的具体值代入，可得的具体值代入，可得 x(n)x(n)y(n)y(n)数字信数字信号处理号处理？ 例：离散傅里叶变换（离散傅里叶变换（DFT）300(0)( )( j)(0)(1)(2)(3)24349nnXx nxxxx 3100123(1)( ) (

14、j)(0) ( j)(1) ( j)(2) ( j)(3) ( j)2 14 ( j)3 ( 1)4j5nnXx nxxxx 3200246(2)( ) ( j)(0) ( j)(1) ( j)(2) ( j)(3) ( j)2 14 ( 1)3 ( 1)4 ( 1)7nnXx nxxxx x(n)x(n)y(n)y(n)数字信数字信号处理号处理？例：例： 离散傅里叶变换（离散傅里叶变换（DFT）因此，得到x(n)的DFT为 3300369(3)( ) ( j)(0) ( j)(1) ( j)(2) ( j)(3) ( j)2 14 (j)3 ( 1)4 ( j)5nnXx nxxxx ( )

15、9,5,7,5X k （1 1）先按）先按N N=8-N/2=4=8-N/2=4，做，做4 4点的点的DFTDFT： 1212/2kNkNX kX kW XkX kNX kW Xk12/, 0Nk 例：求例：求 N N=2=23 3=8=8点点FFTFFT变换变换rkNNrrkNNrrkNNrrkNNrWrxWrxkXWrxWrxkX2/1202/212022/1202/11201) 12()()()2()()( 例：求例：求 N N=2=23 3=8=8点点FFTFFT变换变换 将N=8点DFT分解成2个4点DFT： 时域上时域上：x(0),x(2),x(4),x(6)为偶子序列 x(1),

16、x(3),x(5),x(7)为奇子序列 频域上频域上：X(0)X(3)，由X(k)给出 X(4)X(7),由X(k+N/2)给出将将N=8N=8点分解成点分解成2 2个个4 4点的点的DFTDFT的信号流图的信号流图 4点DFTx1(0)=x(0)x1(1)=x(2)x1(2)=x(4)x1(3)=x(6)4点DFTx2(0)=x(1)x2(1)=x(3)x2(2)=x(5)x2(3)=x(7)08W18W28W38WX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)X1(0)X1(1) 0112812823128128000111 222 333XXXWXXXWXXXWXXXW

17、如：X(4)X(7)=?X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)(2) N/2(4(2) N/2(4点点)-N/4(2)-N/4(2点点)FFT)FFT 若将N/2(4点)子序列按奇/偶分解成两个N/4点(2点)子序列。即对将x1(r)和x2(r)分解成奇、偶两个N/4点(2点)点的子序列。(a)将4点分解成2点的DFT：(2) N/2(4(2) N/2(4点点)-N/4(2)-N/4(2点点)FFT)FFT 104:26xxx rxx、偶序列、 奇序列 13142 0101214xLxLNLLxLXL偶序列若设：，在此，奇序列x1(0)=x(0)；x1(1)=x(2)x1

18、(2)=x(4)；x1(3)=x(6)(2) N/2(4(2) N/2(4点点)-N/4(2)-N/4(2点点)FFT)FFT 215:37xxx rxx、偶序列同理：、奇序列 25252 0101214xLxLNLLxLXL偶序列同理：，在此，奇序列x2(0)=x(1)；x2(1)=x(3)x2(2)=x(5)；x2(3)=x(7)(b) 求求2点的点的DFT (c) 一个2点的DFT蝶形流图2点DFT2点DFTx3(0)=x1(0)=x(0)x3(1)=x1(2)=x(4)X3(0)X3(1)X4(0)X4(1)X1(0)X1(1)X1(2)X1(3)04W14W 011344134401

19、13441344000111200 311XXW XXXW XXXW XXXW X其中2点DFT2点DFTX4(0)X4(1)04W14Wx4(0)=x1(1)=x(2)x4(1)=x1(3)=x(6) (d) 另一个2点的DFT蝶形流图2点DFT2点DFTX5(0)X5(1)X6(0)X6(1)X2(0)X2(1)X2(2)X2(3)04W14W01254625460125462546(0)(0)(0)(1)(1)(1)(2)(0)(0)(3)(1)(1)XXW XXXW XXXW XXXW X同理：同理：其中：x5(0)=x2(0)=x(1)x5(1)=x2(2)=x(5)x6(0)=x2

20、(1)=x(3)x6(1)=x2(3)=x(7)(3) (3) 将将N/4(2N/4(2点点)DFT)DFT再分解成再分解成2 2个个1 1点的点的DFTDFT(a) 求2个一点的DFT 21200000322220100322221220 102200404(1)04040,1;0,1NnknnknknknkNNX kx n WXxWxWxWxWXxWxWxWxWWWWWnkWW 代入上面流图可知：这是一蝶形结这里用到对称性，则，其中 最后剩下两点DFT，它可分解成两个一点DFT，但一点DFT就等于输入信号本身，所以两点DFT可用一个蝶形结表示。取x(0)、x(4)为例。(b) 2个1点的D

21、FT蝶形流图 1点DFTx(0)1点DFTx(4)X3(0)X3(1)02W进一步简化为蝶形流图：02WX3(0)X3(1)x(0)x(4) 032032004 04104 04XxW xxxXxW xxx其中：(4)一个完整一个完整N=8的按时间抽取的按时间抽取FFT的运算流图的运算流图 x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)38W28W18W08W04W02W02W02W02W04W14W14WX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)m=0m=1m=2 N点点DIT-FFT运算流程图运算流程图(N=8)X1(0)X1(1)X1(2)X1(3

22、)X2(0)X2(1)X2(2)X2(3)X3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1) DIT-FFTDIT-FFT算法与直接计算算法与直接计算DFTDFT运算量的比较运算量的比较: : 由前面介绍的DIT-FFT运算流图可见： 每级都由N/2个蝶形单元构成，因此每一级运算都需要N/2次复乘和N次复加（每个结加减各一次）。这样（N=2M)M级运算共需要： 复乘次数： 复加次数： 可以得出如下结论：按时间抽取法所需的复乘数和复加数都是与成正比。而直接计算DFT时所需的复乘数与复加数则都是与N2成正比.(复乘数N2,复加数N(N-1)N2)NNMNCM2log2

23、22 NNMNCA2log2NN2log DIT-FFTDIT-FFT算法与直接计算算法与直接计算DFTDFT运算量的比较运算量的比较: :复数乘法复数乘法：3.3 DIT-FFT DIT-FFT的运算规律及编程思想的运算规律及编程思想1、原位运算：、原位运算：利用同一单元存储蝶形计算的输入、输出数据。每个蝶形的输入和输出均为相同位数。原位运算可节省大量内存，因而硬件成本低。一个完整一个完整N=8的按时间抽取的按时间抽取FFT的运算流图的运算流图 x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)38W28W18W08W04W02W02W02W02W04W14W14WX(0)X(

24、1)X(2)X(3)X(4)X(5)X(6)X(7)m=0m=1m=2 N点点DIT-FFT运算流程图运算流程图(N=8)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)1、原位运算：、原位运算： 蝶形运算的特点蝶形运算的特点: :首先每一个蝶形运算都需要两个输入数据，计算结果也是两个数据，与其它结点的数据无关，其它蝶形运算也与这两结点的数据无关、因此，一个蝶形 运算一旦计算完毕，原输入数据便失效了。这就意味着输出数据可以立即使用原输 人数据结点所占用的内存。原来的数据也就消失了。输

25、出、输人数据利用同一内存单 元的这种蝶形计算称为同位(址)计算。 同址运算的优点同址运算的优点：可以节省存储单元，从而降低对计算机存储量的要求或降低硬件实现的成本。3 3. .3 3 DIT-FFT DIT-FFT的的算法特点算法特点及编程思想及编程思想2、蝶形运算两个节点的距离蝶形运算两个节点的距离：x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)38W28W18W08W04W02W02W02W02W04W14W14WX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)m=0m=1m=2 N点点DIT-FFT运算流程图运算流程图(N=8)X1(0)X1(1)X1(

26、2)X1(3)X2(0)X2(1)X2(2)X2(3)X3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)3 3. .3 3 DIT-FFT DIT-FFT的的算法特点算法特点及编程思想及编程思想2、蝶形运算两个节点的距离蝶形运算两个节点的距离：3 3. .3 3 DIT-FFT DIT-FFT的的算法特点算法特点及编程思想及编程思想3、旋转因子的变化规律：、旋转因子的变化规律：x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)38W28W18W08W04W02W02W02W02W04W14W14WX(0)X(1)X(2)X(3)X(4)X(5)X(6

27、)X(7)m=0m=1m=2X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)2 2(m-1)(m-1)2 2(m-1)(m-1)N=2N=2L L2 2(m-1)(m-1)N=2N=2L L4 4、倒位序倒位序规律规律x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)38W28W18W08W04W02W02W02W02W04W14W14WX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)m=0m=1m=2 N点点DIT-FFT运算流程图运算流程图(N=8)

28、X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)4 4、倒位序倒位序规律规律倒序是有规律的。由于 ，所以顺序数可用M位二进制数（ ）表示。先按n0的0和1进行偶奇分解，再按n1、n2、依次进行分解。 MN20121. nnnnMM4 4、倒位序倒位序规律规律 000 0 100 4 010 2 110 6 001 1 101 5 011 3 111 7 形成倒序的树状图(N=23) 0 00 00 01 11 10 00 00 00 01 11 11 11 11 1(n(n2 2n

29、n1 1n n0 0) )2 2n n0 0为奇数为奇数n n0 0为偶数为偶数二进制数二进制数 十进制数十进制数3 3. .3 3 频率抽取法频率抽取法FFT(DIF-FFT)FFT(DIF-FFT)DITDIT算法算法：把x(n)按奇数偶数分解为越来越短的傅里叶变换；DIFDIF算法算法：把X(k)按奇数偶数分解为越来越短的傅里叶变换；频率抽选基2FFT算法简称为频率抽选，它的推导过程遵循两个规则规则：对时间前后分；对频率偶奇分。3 3. .3 3 频率抽取法频率抽取法FFT(DIF-FFT)FFT(DIF-FFT)设序列x(n)长度为 ，将其前后对半分开，得：MN2 10/2 110/2

30、 NknNnNNknknNNnn NX kDFT x nx n Wx n Wx n W式中 /2 1/2 1/2002NNk n NknNNnnNx n Wx nW/21,11kkNNkWk 偶数，奇数 /2 1/202NkNknNNnNx nWx nW再将X(k)分解成偶数组和奇数组k为偶数时： /2 120/2 1/2022 2NrnNnNrnNnNXrx nx nWNx nx nWk为奇数时： /2 1210/2 1/202122NrnNnNnrnNNnNXrx nx nWNx nx nWW令 122,0,1,.,122nNNx nx nx nNnNxnx nx nW得 /2 11/20

31、/2 12/20221NrnNnNrnNnXrx n WXrxn W12,.,1 , 0Nr DIF-FFT蝶形运算流图nNWx(n)x(n+N/2)x1(n)=x(n)+x(n+N/2)x2(n)=x(n)-x(n+N/2)nNWx(n)x(n)y(n)y(n)数字信数字信号处理号处理？ 例：直接计算直接计算DFT【例例：】 计算有限长序列, N=4的DFT X(k)。【解】 由于x(n)是4点序列，因此X(k)也是4点序列。 ( )2, 4, 3, 4x n 340,03nknX kx n Wk 2j44ejW 所以，若将所以，若将k k的具体值代入，可得的具体值代入，可得 x(n)x(n

32、)y(n)y(n)数字信数字信号处理号处理？ 例：离散傅里叶变换（离散傅里叶变换（DFT）300(0)( )( j)(0)(1)(2)(3)24349nnXx nxxxx 3100123(1)( ) ( j)(0) ( j)(1) ( j)(2) ( j)(3) ( j)2 14 ( j)3 ( 1)4j5nnXx nxxxx 3200246(2)( ) ( j)(0) ( j)(1) ( j)(2) ( j)(3) ( j)2 14 ( 1)3 ( 1)4 ( 1)7nnXx nxxxx x(n)x(n)y(n)y(n)数字信数字信号处理号处理？例：例： 离散傅里叶变换（离散傅里叶变换（DF

33、T）因此，得到x(n)的DFT为 3300369(3)( ) ( j)(0) ( j)(1) ( j)(2) ( j)(3) ( j)2 14 (j)3 ( 1)4 ( j)5nnXx nxxxx ( )9,5,7,5X k DIF-FFT运算流图(N=8)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)0NW1NW2NW2NW2NW3NW0NW0NW0NW0NW0NW0NW38W28W18W08WX(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)04W04W14W14W02W02W02W02Wx(0

34、)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)x3(0)x3(1)x4(0)x4(1)x5(0)x5(1)x6(0)x6(1) N点点DIF-FFT运算流程图运算流程图(N=8)DITDIT与与DIFDIF的异同的异同 运算量相同：L级，每级N/2个蝶形，每个蝶形一次复乘，两次复加； 与时间抽选算法一样，频率抽选FFT算法也具有同址(原位)计算的优点。但是，与时间抽选不同的是，频率抽选FFT算法的信号输入为正序排列，输出为码位倒置排列，因此输出要进行变址计算。 利用转置定理，可以使DIT和DIF的基本蝶

35、形进行相互转换。3 3. .4 4 快速傅里叶反变换快速傅里叶反变换IFFTIFFT的计算的计算 10NknNnX kx n W 101NknNkx nXk WNDFTDFT及及IDFTIDFT的定义的定义0,1,1kN0,1,1nN 比较上面两式，可以看出，只要把DFT公式中 的系数 改为 ，并乘以系数1/N，就可用 FFT算法来计算IDFT，就得到了IFFT的算法。 当把时间抽选FFT算法用于IFFT计算时，由于 原来输入的时间序列x(n)现在变为频率序列 X(k)，原来是将x(n)偶奇分的，而现在变成 对X(k)进行偶奇分了，因此这种算法改称为 频率抽选IFFT算法。类似地，当把频率抽选 FFT算法用于计算IFFT时，称为时间抽选IFFT 算法。 3 3. .4