快速傅里叶和小波变换

1、快速傅里叶变换快速傅里叶变换(FFT)(FFT)陈琤陈琤Email：本讲在分析直接计算本讲在分析直接计算DFT的特点的基础上介绍的特点的基础上介绍DFT的快速算法的快速算法-快速傅里叶变换快速傅里叶变换(FFT)；同时；同时简要介绍了简要介绍了FFT算法的发展历程；此外还要介算法的发展历程；此外还要介绍绍FFT的两种最常用的算法基于时间抽取的两种最常用的算法基于时间抽取的的FFT（DIT：库利图基算法）和基于频率：库利图基算法）和基于频率抽取的抽取的FFT（DIF：桑德图基算法）。：桑德图基算法）。一、直接计算一、直接计算DFTDFT存在的问题存在的问题N点序列点序列x(n)的的DFT变换定义

2、为：变换定义为：1, 2 , 1 , 0 )(1)(1, 2 , 1 , 0 )()(1010NnWkXNnxNkWnxkXNkknNNnknN计算计算X(k)的运算量：需要的运算量：需要N2次复数乘法，次复数乘法，N(N1)次复数加法。在次复数加法。在N较大时计算量很大。较大时计算量很大。例如：例如：N1024时时, 需要需要1,048,576次复数乘法次复数乘法, 即即4,194,304次实数乘法次实数乘法对于象雷达、通信、声纳等需要实时处理的信号，对于象雷达、通信、声纳等需要实时处理的信号，因为其运算量更大，所以无法满足信号处理的实时因为其运算量更大，所以无法满足信号处理的实时性要求。迫

3、切需要有新的算法。性要求。迫切需要有新的算法。 二、二、DFT运算的特点运算的特点实际上，实际上，DFT运算中包含有大量的重复运算。在运算中包含有大量的重复运算。在WN 矩阵中，虽然其中有矩阵中，虽然其中有N2个元素，但由于个元素，但由于WN的周的周期性，其中只有期性，其中只有N个独立的值，即个独立的值，即 ，且这且这N个值也有一些对称关系。总之，个值也有一些对称关系。总之，WN因子具有因子具有如下所述周期性及对称性：如下所述周期性及对称性：110,nNNNWWW1.1.对称性对称性nNkNNnkNknNWWW)()(2.2.周期性周期性)()(knNknNnNkNWWW由上述特性还可得出：由

4、上述特性还可得出：kN)N/(kNNNWWW22/ , 1利用上述对称特性，可使利用上述对称特性，可使DFT运算中有些项可以合运算中有些项可以合并，这样，可使乘法次数减少大约一半；利用并，这样，可使乘法次数减少大约一半；利用WN矩阵的对称性及周期性，可以将长序列的矩阵的对称性及周期性，可以将长序列的DFT分解分解为短序列的为短序列的DFT，N越小，运算量能够减少。越小，运算量能够减少。例如，对于四点的例如，对于四点的DFT，直接计算需要，直接计算需要16次复数乘次复数乘法，根据上述特性可以有以下形式的算法：法，根据上述特性可以有以下形式的算法：, 10W第二列和第三列交换，则：第二列和第三列交

5、换，则：)3()2() 1 ()0(111111111111)3()2() 1 ()0(1111xxxxWWWWXXXX) 3()2() 1 ()0() 3()2() 1 ()0(1010000010100000 xxxxWWWWWWWWWWWWWWWWXXXX则有：则有：由此得出：由此得出：11)3() 1 ()2()0()3()3() 1 ()2()0()2()3() 1 ()2()0() 1 ()3() 1 ()2()0()0(WxxxxXxxxxXWxxxxXxxxxX从上例可知，通过应用对称性和周期性，从上例可知，通过应用对称性和周期性，4点的点的DFT实际上只需要进行一次复数乘法。

6、实际上只需要进行一次复数乘法。)3() 1 ()2()0(111111111111)3()2() 1 ()0(1111xxxxWWWWXXXX三、三、FFT发展简介发展简介FFT的实质：的实质：快速傅里叶变换快速傅里叶变换(FFT)并不是一种新并不是一种新的变换，是为了改进和提高离散傅里叶变换的变换，是为了改进和提高离散傅里叶变换(DFT)运算速度基于运算速度基于DFT运算特点而发展起来的运算特点而发展起来的DFT快快速算法，其实质还是速算法，其实质还是DFT。FFT发展的原因发展的原因：DFT是信号分析与处理中的一是信号分析与处理中的一种重要变换，广泛应用于通信、图像处理、雷达种重要变换，广

7、泛应用于通信、图像处理、雷达及声纳等领域，由于其计算量与变换区间长度及声纳等领域，由于其计算量与变换区间长度N的平方成正比，在的平方成正比，在N较大时，计算量很大，使得较大时，计算量很大，使得直接应用直接应用DFT进行实时处理信号是不现实的。进行实时处理信号是不现实的。FFT的发展历程：的发展历程：1965年，年，J. W. Cooley和和J. W. Tukey巧妙应用巧妙应用DFT中中W因子的周期性及对称性提出了最早的因子的周期性及对称性提出了最早的FFT，这是，这是基于时间抽取的基于时间抽取的FFT。具有里程碑式的贡献。具有里程碑式的贡献(运算量运算量缩短两个数量级缩短两个数量级)196

8、6年，年，G. Sand提出了基于频率抽取的提出了基于频率抽取的FFT算法算法1975年，年，Winogard提出提出WFTA法；法；1977年年Kolha和和Parks提出素因子算法（提出素因子算法（PFA）1984年，年，P. Dohamel和和H. Hollmann提出分裂基快速提出分裂基快速算法，进一步减少了计算量，提高了计算速度（目算法，进一步减少了计算量，提高了计算速度（目前最理想的算法）前最理想的算法）FFT的各种算法的各种算法纵观纵观FFT的发展历程，的发展历程，FFT算法分成算法分成两大类两大类：(1) 针对针对N等于等于2的整数次幂的算法，如基的整数次幂的算法，如基2算法、

9、基算法、基4算法、实因子算法和分裂基算法算法、实因子算法和分裂基算法(2)针对针对N不等于不等于2的整数次幂的算法，其以的整数次幂的算法，其以Winograd为代表的一类算法为代表的一类算法(素因子法素因子法PFA、Winograd算法算法WFTA)简要介绍库利简要介绍库利-图基算法和桑得图基算法和桑得-图基算法图基算法12/02212/0112/012/0)12(212/0212/01)()() 12()()()()(NnkrNkNNnknNNrNrrkNkrNNnknNNnknNWrxWWrxWrxWrxWnxWnxkX设序列设序列x(n)的长度为的长度为N，且满足，且满足N2M，M为自然

10、为自然数，按数，按n的奇偶将的奇偶将x(n)分解为两个分解为两个N/2的子序列：的子序列：x1(r)=x(2r), r=0,1,2,N/2-1x2(r)=x(2r+1) r=0,1,2,N/2-1则则x(n)的的DFT为：为：1.基本原理基本原理四、按时间抽取四、按时间抽取(DIT)的的FFT库利库利-图基算法图基算法因为因为 ，所以：，所以：krNkrNjkrNjkrNWeeW2/2/2222krNNrkrNNrWrxkXWrxkX2/12/0222/12/011)( )()()(上式中上式中X1(k)和和X2(k)分别为分别为x2(r)和和x2(r)的的N/2点点DFT，即即1, 1 ,

11、0 )()()()()(212/12/022/12/01NkkXWkXWrxWWrxkXkNkrNNrkNkrNNr由于由于X1(k)和和X2(k)均以均以N/2为周期，且为周期，且 ，以，以X(k)又可表示为：又可表示为：kNNkNWW2/即将一个即将一个N点的点的DFT分解成为两个分解成为两个N/2点的点的DFT。上述运算可用右下图来表示，称为蝶形运算符号。上述运算可用右下图来表示，称为蝶形运算符号。12/, 1 , 0(2) )()()(1) )()()(21221NkkXWkXkXkXWkXkXkNNkN从右图可知，要完成一从右图可知，要完成一个蝶形运算需要进行一个蝶形运算需要进行一次

12、复数相乘和两次复数次复数相乘和两次复数相加运算。相加运算。下图是下图是N8时的一个分解运算图。时的一个分解运算图。从上图可知，经过一次分解后，计算一个从上图可知，经过一次分解后，计算一个N点的点的DFT共需要计算两个共需要计算两个N/2点点FFT和和N/2个蝶形运算。个蝶形运算。计算一个计算一个N/2点点DFT需要需要(N/2)2复数乘和复数乘和N/2(N/2-1)次复数加法。所以按刚才的方法计算次复数加法。所以按刚才的方法计算N点点DFT总的总的运算量为运算量为2(N/2)2+N/2=N(N+1)/2N2/2(N1时时)复数复数次乘法和次乘法和N(N/2-1)+2N/2=N2/2次复数加法运

13、算。次复数加法运算。由此可见，由此可见，仅仅经过一次分解就能使运算量减少近仅仅经过一次分解就能使运算量减少近一半！一半！因为因为N/2仍然是偶数，可以作进一步的分解：仍然是偶数，可以作进一步的分解：与第一次分解相同，将与第一次分解相同，将x1(r)按奇偶分解成两个按奇偶分解成两个N/4的子序列的子序列x3(l)和和x4(l), 即：即：则，则，X1(k)又可表示为：又可表示为：14/, 1 , 0 )()()()() 12()2()(42/314/04/42/14/04/314/0)12(2/14/04/1NkkXWkXWlxWWlxWlxWlxkXkNNlklNkNNlklNNllkNNlk

14、lN同理，同理，X3(k)和和X4(k)的周期性和的周期性和WN的对称性，到最的对称性，到最后我们能够得到：后我们能够得到：14/, 1 , 0 ,) 12()()2()(1413Nllxlxlxlx同理可得：同理可得：14/, 1 , 0 ,)()()4/()()()(62/5262/52NkkXWkXNkXkXWkXkXkNkN其中：其中：14/, 1 , 0 ,) 12()()2()()()()()(262564/14/06654/14/055NllxlxlxlxlxDFTWxkXlxDFTWxkXklNNlklNNl14/, 1 , 0 ,)()()4/()()()(42/3142/3

15、1NkkXWkXNkXkXWkXkXkNkN这样，又将这样，又将N/2点的点的DFT分解为两个分解为两个N/4点的点的DFT。依次类推，经过依次类推，经过M-1次分解，最后将次分解，最后将N点点DFT分解分解成成N/2个个2点点DFT。一个完整的。一个完整的8点点DFT-FFT运算流运算流图如下图所示。图如下图所示。2.运算量的比较运算量的比较从上述分析过程可知，在从上述分析过程可知，在N2M时，每一级都由时，每一级都由N/2个蝶形运算构成，即每级都需要个蝶形运算构成，即每级都需要N/2次复数乘次复数乘和和N次复数加，所以总的复数乘的次数为：次复数加，所以总的复数乘的次数为：NNMN2logN

16、NMN2log22总的复数加的次数为：总的复数加的次数为：直接计算时复数乘的次数为直接计算时复数乘的次数为N2，加为，加为N(N-1)次。当次。当N1时，时， ，使运算量大大减少。，使运算量大大减少。NNN22log2以以N1024为例，其运算量与直接计算的比例为：为例，其运算量与直接计算的比例为：8 .20451201048576log)2/(22NNN即运算效率提高了即运算效率提高了200多倍。易知多倍。易知N越大，优越性越大，优越性越明显。另外，在越明显。另外，在N2048时，直接运算需要时，直接运算需要3个个小时，而采用小时，而采用FFT则只需不到一分钟就能完成！则只需不到一分钟就能完

17、成！3.DITFFT的运算规律的运算规律(1)原位计算原位计算根据运算流图可知，根据运算流图可知，DIT-FFT的运算很有规律。的运算很有规律。N2M点的点的FFT共进行共进行M级运算，每级运算有级运算，每级运算有N/2个个蝶形运算构成；同一级中，每个蝶形的两个输入数蝶形运算构成；同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，并且每个蝶形的输入、输据只对计算本蝶形有用，并且每个蝶形的输入、输出数据节点又同在一条水平线上，这意味着计算完出数据节点又同在一条水平线上，这意味着计算完一个蝶形后所得数据可立即存入原输入数据所占用一个蝶形后所得数据可立即存入原输入数据所占用的存贮单元，这样，经过的

18、存贮单元，这样，经过M级运算后，原来存放输级运算后，原来存放输入序列数据的入序列数据的N个存贮单元中并依次存放了个存贮单元中并依次存放了X(k)的的N个值。这种利用同一存贮单元存贮计算输入、输出个值。这种利用同一存贮单元存贮计算输入、输出数据的方法称为原位（址）计算。数据的方法称为原位（址）计算。 (0)=(0)=X X0 0(0)(0) X X1 1(0)(0) X X2 2(0) X(0) X3 3(0)(0)=X(0) =X(0) (4)=(4)=X X0 0(1)(1) X X1 1(1) X(1) X2 2(1) X(1) X3 3(1)(1)=X(1)=X(1) (2)=(2)=X

19、 X0 0(2)(2) X X3 3(2)(2)=X(2)=X(2) (6)=(6)=X X0 0(3)(3) X X3 3(3)(3)=X(3)=X(3) (1)=(1)=X X0 0(4)(4) X X1 1(4) X(4) X2 2(4) X(4) X3 3(4)(4)=X(4)=X(4) (5)=(5)=X X0 0(5)(5) X X3 3(5)(5)=X(5)=X(5) (3)=(3)=X X0 0(6)(6) X X3 3(6)(6)=X(6)=X(6) (7)=(7)=X X0 0(7)(7) X X1 1(7) X(7) X2 2(7) X(7) X3 3(7)(7)=X(7

20、)=X(7)WWWWN0N0N0N0-1-1-1-1WWWWN0N2N0N2-1-1-1-1WWWWNNNN0123.xxxxxxxx输入数据、中间运算结果和最后输出均用同一存储器。输入数据、中间运算结果和最后输出均用同一存储器。(2)旋转因子的变化规律旋转因子的变化规律在每个蝶形的运算过程中，都要乘以因子在每个蝶形的运算过程中，都要乘以因子 ，称，称其为旋转因子，其为旋转因子，p称为旋转因子指数。但各级的旋称为旋转因子指数。但各级的旋转因子和循环方式都有所不同。转因子和循环方式都有所不同。pNW旋转因子与运算级数有一定的关系，若用旋转因子与运算级数有一定的关系，若用L表示表示运算级数，对于运

21、算级数，对于N2M的一般情况，第的一般情况，第L级的旋转级的旋转因子为：因子为：1,20,1,J , 2222 12 , 2 , 1 , 0 ,1 -L)(2L12lMJNJNpNMLMLMLJpNWWWNJWWMlLLMJp2 从运算流图可以看出，原位计算时，从运算流图可以看出，原位计算时，FFT的的输出输出X(k)是按正常顺序排列在存储单元中，即按是按正常顺序排列在存储单元中，即按X(0)，X(1),，X(7)的顺序排列，但是这时输入的顺序排列，但是这时输入x(n)都不是按自然顺序存储的，这看起来好象是都不是按自然顺序存储的，这看起来好象是“混乱无序混乱无序”的，实际上是有规律的，我们称之

22、的，实际上是有规律的，我们称之为倒位序。为倒位序。 造成倒位序的原因是输入造成倒位序的原因是输入x(n)按标号按标号n的偶的偶奇的不断分组而造成。奇的不断分组而造成。(3)倒位序规律倒位序规律倒位序实现倒位序实现 输入序列先按自然顺序存入存储单元输入序列先按自然顺序存入存储单元, ,然后经然后经变址运算来实现变址运算来实现倒位序排列倒位序排列，设输入设输入序列的序号序列的序号为为n,n,二进制为二进制为(n(n2 2 n n1 1 n n0 0 ) )2 2 , ,倒位序倒位序顺序用顺序用 表示表示, ,其其倒位序倒位序二进制为二进制为(n(n0 0 n n1 1 n n2 2 ) )2 2

23、。n A(1) A(2) A(3) A(4) A(5) A(6) A(7) A(8)x(0) x(1) x(2) x(3) x(4) x(5) x(6) x(7)x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7)变址处理方法变址处理方法存储单元存储单元自然顺序自然顺序变址变址倒位序倒位序 0 00 0 0 0 0 00 0 0 0 0 0 0 0 1 01 0 0 0 1 11 1 0 0 0 4 0 4 2 02 0 1 1 0 00 0 1 1 0 2 0 2 3 03 0 1 1 1 11 1 1 1 0 6 0 6 4 14 1 0 0 0 00 0 0 0

24、 1 1 1 1 5 15 1 0 0 1 11 1 0 0 1 5 1 5 6 16 1 1 1 0 00 0 1 1 1 3 1 3 7 17 1 1 1 1 11 1 1 1 1 7 1 7 自然顺序自然顺序n n 二进制二进制n n n n n n 倒位序二进制倒位序二进制n n n n n n 倒位顺序倒位顺序n n2 1 0 0 1 2例如例如 ，N=8时如下表：时如下表： 4.4.蝶形运算两节点的距离蝶形运算两节点的距离:2:2m-1m-1 其中其中,m,m表示第表示第mm列列, ,且且m =1, ,Lm =1, ,L 例如例如N=8=2N=8=23 3 , ,第一级第一级( (

25、列列) )距离为距离为2 21-11-1=1,=1, 第二级第二级( (列列) )距离为距离为2 22-12-1=2=2， 第三级第三级( (列列) )距离为距离为2 23-13-1=4=4。五、按频率抽取五、按频率抽取(DIF)的的FFT桑德桑德-图基算法图基算法库利图基法是将输入序列按其顺序是奇数还是偶库利图基法是将输入序列按其顺序是奇数还是偶数来分解为越来越短的序列；桑德图基法是把输数来分解为越来越短的序列；桑德图基法是把输出序列出序列X(k)按其顺序的偶奇来分解为越来越短的按其顺序的偶奇来分解为越来越短的序列。序列。1.原理原理设序列设序列x(n)长度为长度为N2M，首先将，首先将x(

26、n)前后对半分前后对半分开，得到两个子序列，其开，得到两个子序列，其DFT可表示为：可表示为：10)()(NnnkNWnxkX10)(1012/102222)2()()()(NNNNnknNnnkNNNnnkNnnkNWNnxWnxWnxWnx 由于由于 故故 因此因此 X(k)X(k)可表示为可表示为 nkNnkNWWNnxnxNN1022)2()(, 12jNeWNkkNNW) 1(2nkNnkWNnxnxkXN102)2() 1()()(1, 1 , 0Nk 2.N2.N点点DFTDFT按按k k的奇偶分组可分为两个的奇偶分组可分为两个N/2N/2的的DFTDFT 当当k k为偶数为偶数

27、, ,即即 k=2rk=2r时时,(-1),(-1)k k =1 =1； 当当k k为奇数为奇数, ,即即 k=2r+1 k=2r+1 时时 (-1)(-1)k k =-1 =-1 。 这时这时 X(k)X(k) 可分为两部分：可分为两部分： nrnnrNnNNNWNnxnxWNnxnx22210210)2()()2()()2( rX1, 1 , 02Nrk k为偶数时：为偶数时： 可见,上面两式均为N/2的DFT。nrnnNrnNnNNNWWNnxnxWNnxnxrX22210)12(10)2()()2()() 12(k k为奇数时：为奇数时：1, 1 , 02Nr-1-1)2()(Nnxn

28、x1, 1 ,02NnnNWNnxnx)2()(3.3.蝶形运算蝶形运算进行如下碟形运算：和)2()(NnxnxnNW)2(Nnx)(nx-1-1-1-1-1-1-1-1W WW WW WW WN NN NN NN N0 01 12 23 34.N=84.N=8时时, ,按按k k的奇偶分解过程的奇偶分解过程 先蝶形运算，后先蝶形运算，后DFT:DFT:)0( x) 1 ( x)5( x)4( x) 3( x)2( x)7( x)6( x)0(X)2(X)6(X) 1 (X) 3(X)5(X)7(X)4(XDFTN点2DFTN点2 5. 5.仿照仿照DITDIT的方法的方法 再将再将N/2N/

29、2点点DFTDFT按按k k的奇偶分解为两个的奇偶分解为两个N/4N/4点的点的DFT,DFT,如此进行下去如此进行下去, ,直至分解为直至分解为2 2点点DFTDFT。 以下是以下是8 8点的点的DIF-DFTDIF-DFT流程：流程： x(0) X(0) X(0) x(1) X(4) X(4) x(2) X(2) X(2) x(3) X(6) X(6) x(4) X(1) X(1) x(5) X(5) X(5) x(6) X(3) X(3) x(7) X(7) X(7)-1-1-1-1-1-1-1-1W WW WW WW WN NN NN NN N0 01 12 23 3-1-1-1-1-

30、1-1-1-1W WW WW WW WN NN NN NN N0 02 20 02 2-1-1-1-1-1-1-1-1W WW WW WW WN NN NN NN N0 00 00 00 0 B B. .原位运算原位运算 每级每级( (列列) )都是由都是由N/2N/2个蝶形运算构成个蝶形运算构成, ,即即 -1-1W WN Nr rrNmmmmmmWjXkXjXjXkXkX)()()()()()(1111)(1kXm)(1jXm)()()(11jXkXkXmmmrNmmmWjXkXjX)()()(11C.C.蝶形运算两节点的距离蝶形运算两节点的距离 一般公式为一般公式为2 2L-mL-m =

31、N/2 =N/2mm 例如例如 N=2N=23 3 =8 =8 ： (1)m=1 (1)m=1 时的距离为时的距离为 8/2=48/2=4； (2)(2)m=2 m=2 时的距离为时的距离为 8/4=2;8/4=2; (3)(3)m=3 m=3 时的距离为时的距离为 8/8=18/8=1。1.1.相同点相同点 (1)(1)进行原位运算；进行原位运算； (2)(2)运算量相同运算量相同, ,均为（均为（N/2) LogN/2) Log2 2N N次复乘、次复乘、N N LogLog2 2N N次复加。次复加。2.2.不同点不同点 (1)DIT(1)DIT输入为倒位序输入为倒位序, ,输出为自然顺

32、序；输出为自然顺序； DIFDIF正好与此相反。但正好与此相反。但DITDIT也有输入为自然顺序也有输入为自然顺序, ,输出为倒位序的情况。输出为倒位序的情况。D.DIFD.DIF法与法与DITDIT法的异同法的异同rNmmmWjXkXjX)()()(11rNmmmWjXkXkX)()()(11)(1kXm)(1jXmrNW1(2)(2)蝶形运算不同蝶形运算不同A.A.DITDIT用矩阵表示)(kXm)(jXmrNWrNW)(1kXm)(1jXm=11B.B.DIFDIF用矩阵表示)(kXm)(jXmrNWrNW)(1kXm)(1jXm=11rNmmmWjXkXjX)()()(11)()()(

33、11jXkXkXmmm)(1kXm)(1jXmrNW1(3)(3)两种蝶形运算的关系两种蝶形运算的关系互为转置（矩阵）；互为转置（矩阵）； 将流图的所有支路方向都反向将流图的所有支路方向都反向, ,交换输入和输交换输入和输出，即可得到另一种蝶形。出，即可得到另一种蝶形。 A.A. DIT DITB.B.DIFDIFrNWrNW1111rNWrNWFFT算法，同样可以适用于离散博里叶反变换算法，同样可以适用于离散博里叶反变换(IDFT)运算，并简称为运算，并简称为IFFT，即快速博里叶反变，即快速博里叶反变换，从换，从IDFT公式看：公式看：10)()()(NnknNWnxnxDFTkX而而DF

34、T公式公式10)(1)()(NkknNWkXNkXIDFTnx比较以上两式可知，只要把比较以上两式可知，只要把DFT运算中的每一个运算中的每一个系数系数 换成换成 ，并且最后再乘以常数，并且最后再乘以常数1/N, 则则knNWknNW六、六、IDFT的快速算法的快速算法-IFFT以上所有时间抽取或频率抽取的以上所有时间抽取或频率抽取的FFT算法都可以算法都可以拿来运算拿来运算IDFT。此外，有一种直接采用此外，有一种直接采用FFT程序计算程序计算IFFT的方法。的方法。对对 取共轭，则：取共轭，则： 10)(1)()(NkknNWkXNkXIDFTnx*10*)(1)(1)(kXDFTNWkXNnxNkknN两边同时再取共轭，有：两边同时再取共轭，有：10*)(1)(NkknNWkXNnx上述情况表明：在进行上述情况表明：在进行ID