fft算法原理

第四章快速傅立叶变换FastFourierTransform第一节直接计算DFT的问题及改进途径1、问题的提出设有限长序列x(n)，非零值长度为N，若对x(n)进行一次DFT运算，共需多大的运算工作量？计算成本?计算速度?2.DFT的运算量回忆DFT和IDFT的变换式：10)(1)]([)(10NnWkXNkXIDFTnxNknk10)()]([)(10NkWnxnxDFTkXNnnkN1）x(n)为复数，也为复数。2）DFT与IDFT的计算量相当。nkNjnkNeW2注意：计算机运算时（编程实现）：0k0)1(0100)1()1()0()0(NNNNWNxWxWxX1k0111(1)1(1)(0)(1)(1)NNNNXxWxWxNW2k0212(1)2(2)(0)(1)(1)NNNNXxWxWxNW1Nk0111(1)1(1)(0)(1)(1)NNNNNNNXNxWxWxNWN次复乘，N-1次复加N个点10)()]([)(10NkWnxnxDFTkXNnnkN以DFT为例：复数乘法复数加法一个X(k)NN–1N个X(k)(N点DFT)N2N(N–1)实数乘法实数加法一次复乘42一次复加2一个X(k)4N2N+2(N–1)=2(2N–1)N个X(k)(N点DFT)4N22N(2N–1)10()NnkNnxnW运算量(a+jb)(c+jd)=(ac-bd)+j(bc+ad)例：计算一个N点DFT，共需N2次复乘。以做一次复乘1μs计，若N=4096，所需时间为ss1716777216)4096(2例：石油勘探，有24个通道的记录，每通道波形记录长度为5秒，若每秒抽样500点/秒，1）每道总抽样点数：500*5=2500点2）24道总抽样点数：24*2500=6万点3）DFT复乘运算时间：N2=(60000)2=36*108次ss360010*36)60000(82由于计算量大，且要求相当大的内存，难以实现实时处理，限制了DFT的应用。长期以来，人们一直在寻求一种能提高DFT运算速度的方法。FFT便是Cooley&Tukey在1965年提出的的快速算法，它可以使运算速度提高几百倍，从而使数字信号处理学科成为一个新兴的应用学科。第二节改善DFT运算效率的基本途径knNW1、利用DFT运算的系数的固有对称性和周期性，改善DFT的运算效率。1）对称性2）周期性3）可约性()()nkNnknNkNNN周期性nkmnkNmNWW可约性//nknkmNNmWW2jmnkmNe221NjjNee0/2(/2)11NkNkNNNN特殊点：nkNW的特性2jnknkNNWe*()()()nknkNnknNkNNNN对称性NknkNNWWnNnkNNWW2、将长序列DFT利用对称性和周期性分解为短序列DFT的思路因为DFT的运算量与N2成正比的，如果一个大点数N的DFT能分解为若干小点数DFT的组合，则显然可以达到减少运算工作量的效果。N点DFTN/2点DFTN/2点DFTN/4点DFTN/4点DFTN/4点DFTN/4点DFT…….复乘：2N2222NN22N22224444NNNN42NFFT算法的基本思想：利用DFT系数的特性，合并DFT运算中的某些项把长序列DFT→短序列DFT，从而减少运算量。FFT算法分类:时间抽选法DIT:Decimation-In-Time频率抽选法DIF:Decimation-In-Frequency第三节按时间抽选的基2-FFT算法1、算法原理设输入序列长度为N=2M(M为正整数，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的FFT算法。也称为Coolkey-Tukey算法。其中基2表示：N=2M，M为整数.若不满足这个条件，可以人为地加上若干零值（加零补长）使其达到N=2M。先将x(n)按n的奇偶分为两组，作变量置换:当n=偶数时，令n=2r;当n=奇数时，令n=2r+1;分组，变量置换2、算法步骤10)()]([)(10NkWnxnxDFTkXNnnkN得到：1,...,0)()12()()2(221Nrrxrxrxrx带入DFT中10)()]([)(NnnkNWnxnxDFTkX12/0)12(12/02)12()2(NrkrNNrrkNWrxWrx1010)()(NnnkNNnnkNnnWnxWnx为奇数为偶数12/02212/021)()(NrrkNkNNrrkNWrxWWrx所以12/02212/021)()()(NrrkNkNNrrkNWrxWWrxkX由于nNnNjnNjnNWeeW2/2/2222－－12/02/212/02/1)()(NrrkNkNNrrkNWrxWWrx)()(21kXWkXkN1,,1,02Nk1,,1,0Nk?1,,1,02NkX1(k)、X2(k)只有N/2个点，以N/2为周期；而X(k)却有N个点，以N为周期。要用X1(k)、X2(k)表达全部的X(k)值，还必须利用WN系数的周期特性。rkNkNrNWW2/)2/(2/12/02/112/0)2/(2/11)()()2/(NrrkNNrkNrNWrxWrxkNX)()2/()()2/(2211kXkNXkXkNX后半部分前半部分又考虑到的对称性：kNWkNkNNNkNN2/)2/()2/()2/()2/(2)2/(1kNXWkNXkNXkNN有：1,,1,0)()()(221NkNkkXWkXkX1,,1,0)()(221NkNkkXWkX后半部分前半部分)()()2/(21kXWkXkNXkN)()()(21kXWkXkXkN1,,1,02Nk)(1kX)(2kXkNW)()(21kXWkXkN)()(21kXWkXkN蝶形运算流图符号说明：(1)左边两路为输入(2)右边两路为输出(3)中间以一个小圆表示加、减运算（右上路为相加输出、右下路为相减输出)1个蝶形运算需要1次复乘，2次复加复数乘法复数加法一个N点DFTN2N(N–1)一个N/2点DFT(N/2)2N/2(N/2–1)两个N/2点DFTN2/2N(N/2–1)一个蝶形12N/2个蝶形N/2N总计N2/2+N/2≈N2/2N(N/2-1)+N≈N2/2运算量减少了近一半分解后的运算量：)()()2/()()()(2121kXWkXNkXkXWkXkXkNkN12/,,0Nk先将N=8点的DFT分解成2个4点DFT：可知：时域上：x(0),x(2),x(4),x(6)为偶子序列x(1),x(3),x(5),x(7)为奇子序列频域上：X(0)~X(3)，由X(k)给出X(4)~X(7)，由X(k+N/2)给出例子：求N=23=8点FFT变换按N=8→N/2=4，做4点的DFT：N=8点的直接DFT的计算量为：复乘：N2次=64次复加：N(N-1)次=8×7=56次)()()2/()()()(2121kXWkXNkXkXWkXkXkNkN12/,,0Nk此外，还有4个蝶形结，每个蝶形结需要1次复乘，2次复加。一共是：复乘4次，复加8次。得到X1(k)和X2(k)需要：复乘：(N/2)2+(N/2)2次=32次复加：N/2(N/2-1)+N/2(N/2-1)=12+12=24次用分解的方法得到X(k)需要：复乘：32+4=36次复加：24+8=32次N点DFT的一次时域抽取分解图(N=8)4点DFT4点DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)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)38W28W18W08W奇序列、偶序列、)6()2()4()0(:)(1xxxxrx奇序列、偶序列、同理：)7()3()5()1(:)(2xxxxrx因为4点DFT还是比较麻烦，所以再继续分解。若将N/2(4点)子序列按奇/偶分解成两个N/4点(2点)子序列。即对将x1(r)和x2(r)分解成奇、偶两个N/4点(2点)点的子序列。1,0)1...0()()12()()2(44131lllxlxlxlxN此处，奇序列偶序列设：1,0)1...0()()12()()2(46252lllxlxlxlxN此处，奇序列偶序列设：那么，X1(k)又可表示为14/0)12(2/114/022/11)12()2()(NlklNNllkNWlxWlxkX14/04/42/14/04/3)()(NllkNkNNllkNWlxWWlx)()(42/3kXWkXKN1,......1,0)()()()()()(442/34142/31NkNNkNkkXWkXkXkXWkXkX14/0)12(2/214/022/22)21()2()(NlklNNllkNWlxWlxkX14/04/62/14/04/5)()(NllkNkNNllkNWlxWWlx1,......1,0)()()()()()(462/54262/52NkNNkNkkXWkXkXkXWkXkXX2(k)也可以进行相同的分解：注意：通常我们会把写成。kNW2/kNW2)()(62/5kXWkXKNN点DFT的第二次时域抽取分解图(N=8)2点DFT2点DFT2点DFT2点DFTx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)08W28W08W28WX1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)38W28W18W08WX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)4点DFT4点DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)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)38W28W18W08W)1()0()()]([)(32314/04/333xWxWlxlxDFTkXkNlklN)1()0()1()0()1()1()0()0(30233123330233xWxxWxXxWxX8808WX3(0)X3(1)x(0)=x3(0)x(4)=x3(1)N点DIT―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)0NW0NW0NW0NW0NW2NW0NW2NW0NW2NW1NW3NW3、DIT―FFT算法与直接计算DFT运算量的比较22log2NNN1)、N=2M的DFT运算可分成M级，每一级有N/2个蝶形，每个蝶形有一次复乘两次复加。NN2log2NN2log2)、所以M级共有次复乘和次复加。3)、若直接计算DFT，需N2次复乘和N(N-1)次复加。显然，当N较大时，有：例如，N=210=1024时221048576204.8(/2)log5120NNNFFT算法与直接计算