1、Fast Fourier Transform FFTFFT产生故事产生故事 虽然频谱分析和虽然频谱分析和DFT运算很重要，但在很长一段时间里运算很重要，但在很长一段时间里，并没有得到真正的运用。而频谱分析仍大多采用模拟信，并没有得到真正的运用。而频谱分析仍大多采用模拟信号滤波的方法解决。号滤波的方法解决。Cooley J W,Tukey J W.An algorithm for the machine computation of complex Fourier series.Mathematics of Computation,1965,pp297301 提出一种快速计算提出一种快速计算DF
2、T的方法，揭开了的方法，揭开了FFT发展史上发展史上的第一页。的第一页。1984年，法国的杜梅尔和霍尔曼将基年，法国的杜梅尔和霍尔曼将基2分解和基分解和基4分解分解糅合在一起，提出了分裂基糅合在一起，提出了分裂基FFT算法算法,使运算效率进上步提使运算效率进上步提高。高。一、直接计算一、直接计算DFTDFT的问题的问题1、问题的提出、问题的提出 设有限长序列设有限长序列x(n)，非零值长度为非零值长度为N，若对若对x(n)进行一次进行一次DFT运算，共需运算，共需多大的运算工作量多大的运算工作量？2.DFT的运算量的运算量 10)(1)()(10 NnWkXNkXIDFTnxNknk10)()
3、()(10 NkWnxnxDFTkXNnnkN1）x(n)为为复数复数，也为也为复数复数。2）DFT与与IDFT的的计算量相当。计算量相当。nkNjnkNeW 2 注意：注意：计算机运算时（编程实现）：计算机运算时（编程实现）：0k0)1(0100)1()1()0()0(NNNNWNxWxWxX1k 0111(1)1(1)(0)(1)(1)NNNNXxWxWx NW 2k 0 212(1)2(2)(0)(1)(1)NNNNXxWxWx NW 1Nk0111(1)1(1)(0)(1)(1)NNNNNNNX NxWxWx NW 10)()()(10 NkWnxnxDFTkXNnnkN以以DFT为例
4、：为例：复数乘法复数乘法复数加法复数加法一个一个X(k)NN 1N个个X(k)(N点点DFT)N 2N(N 1)实数乘法实数乘法实数加法实数加法一次复乘一次复乘42一次复加一次复加2一个一个X(k)4N2N+2(N 1)=2(2N 1)N个个X(k)(N点点DFT)4N 22N(2N 1)10()()NnkNnX kx n W运算量运算量(a+jb)(c+jd)=(ac-bd)+j(bc+ad)例：例：设做一次复乘用时设做一次复乘用时1s。石油勘探，有石油勘探，有24个通道的个通道的记录记录 ，每通道波形记录长度为，每通道波形记录长度为5秒，若每秒抽样秒，若每秒抽样500点点/秒，秒，1）每道
5、总抽样点数：）每道总抽样点数：500*5=2500点点 2）24道总抽样点数：道总抽样点数：24*2500=6万点万点 3）DFT复乘运算时间：复乘运算时间：N2=(60000)2=36*108次次ss360010*36)60000(82 所需时间为所需时间为 由于计算量大，且要求由于计算量大，且要求相当大的内存相当大的内存，难以实现实难以实现实时处理时处理，限制了，限制了DFT的应用。长期以来，人们一直在寻的应用。长期以来，人们一直在寻求一种能求一种能提高提高DFT运算速度运算速度的方法。的方法。1、利用、利用DFT运算的系数的固有对称性和周期运算的系数的固有对称性和周期 性，合并性，合并D
6、FT中中的某些项。的某些项。1）对称性）对称性 2）周期性周期性 3）可约性）可约性二、改善二、改善DFTDFT运算效率的基本途径运算效率的基本途径)(rNnNnNWW*2(),NmmmmN mmNNNNNNWWWWWW rmmrNNWW2、将长序列、将长序列DFT利用对称性和周期性分解为短利用对称性和周期性分解为短 序列序列DFT来来减少计算量。减少计算量。因为因为DFT的运算量与的运算量与N2成正比的，如果一个成正比的，如果一个大大点数点数N的的DFT能分解为能分解为若干小点数若干小点数DFT的组合的组合，则，则显然可以达到显然可以达到减少运算工作量减少运算工作量的效果。的效果。N点点DF
7、TN/2点点DFTN/2点点DFTN/4点点DFTN/4点点DFTN/4点点DFTN/4点点DFT.复乘：复乘：2N2222 NN22N22224444 NNNN42NFFT算法分类算法分类:时域抽取法时域抽取法 DIT:Decimation-In-Time频域抽取法频域抽取法 DIF:Decimation-In-Frequency三、时域抽取法的基三、时域抽取法的基2-2-FFTFFT算法算法1、算法原理、算法原理 设输入序列长度为设输入序列长度为N=2M（M为正整数为正整数），将该序列，将该序列按时按时间顺序的奇偶分解间顺序的奇偶分解为越来越短的子序列，称为基为越来越短的子序列，称为基2按
8、时间抽取按时间抽取的的FFT算法。也称为算法。也称为Coolkey-Tukey算法。算法。其中其中基基2表示：表示：N=2M，M为整数为整数.若不满足这个条件，可若不满足这个条件，可以人为地加上若干零值（以人为地加上若干零值（加零补长加零补长）使其达到）使其达到 N=2M。先将先将x(n)按按n的奇偶分为两组，作变量置换的奇偶分为两组，作变量置换:当当n=偶数时，令偶数时，令n=2r;当当n=奇数时，令奇数时，令n=2r+1;分组，变量置换分组，变量置换2、算法步骤、算法步骤得到：得到：1,.,0)()12()()2(221 Nrrxrxrxrx 带入带入DFT中中10()()()NnkNnX
9、 kDFT x nx n W1100()()nnNNnknkNNnnx n Wx n W为偶数为奇数 12/0)12(12/02)12()2(NrkrNNrrkNWrxWrx 12/02212/021)()(NrrkNkNNrrkNWrxWWrx所以所以 12/02212/021)()()(NrrkNkNNrrkNWrxWWrxkX由于由于nNnNjnNjnNWeeW2/2/2222 12/02/212/02/1)()(NrrkNkNNrrkNWrxWWrx)()(21kXWkXkN 1,1,02 Nk1,1,0 Nk1,1,02 Nk?X1(k)、X2(k)只有只有N/2个点，以个点，以N/
10、2为周期；为周期；1122(/2)()(/2)()X NkX kX NkX k即后半部分后半部分前半部分前半部分又考虑到又考虑到 的对称性：的对称性：kNWkNkNNNkNNWWWW 2/)2/()2/()2/()2/(2)2/(1kNXWkNXkNXkNN 有：有：1,1,0)()()(221 NkNkkXWkXkX1,1,0)()(221 NkNkkXWkX后半部分后半部分前半部分前半部分)()()2/(21kXWkXkNXkN )()()(21kXWkXkXkN 1,1,02 Nk蝶形运算流图符号蝶形运算流图符号-1-1)()()(21kXWkXkXkN)()()2(21kXWkXkNX
11、kN)(1kX)(2kXkNW复数乘法复数乘法复数加法复数加法一个一个N 点点DFTN 2N(N1)一个一个N/2点点DFT(N/2)2N/2(N/2 1)两个两个N/2点点DFTN 2/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(
12、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=8N/2=4，做做4点的点的DFT：N点点DFT的一次时域抽取分解图的一次时域抽取分解图(N=8)x1 1(0)=(0)=x(0)(0)x1(1)=(1)=x(2)(2)N/2N/2点点 x1(2)=(2)=x(4)DFT (4)DFT x1(3)=(3)=x(6)(6)x2(0)=(0)=x(1)(1)x2(1)=(1)=x(3)(3)N/2N/
13、2点点 x2(2)=(2)=x(5)(5)DFT DFT x2(3)=(3)=x(7)(7)X X1(0)(0)X X1(1)(1)X X1(2)(2)X X1(3)(3)X X2(0)(0)X X2(1)(1)X X2(2)(2)X X2(3)(3)WN2WN1WN0WN3-1-1-1-1X(0)X(0)X(1)X(1)X(2)X(2)X(3)X(3)X(4)X(4)X(5)X(5)X(6)X(6)X(7)X(7)奇奇序序列列、偶偶序序列列、)6()2()4()0(:)(1xxxxrx 奇奇序序列列、偶偶序序列列、同同理理：)7()3()5()1(:)(2xxxxrx因为因为4点点DFT还是
14、比较麻烦，所以再继续分解。还是比较麻烦，所以再继续分解。若将若将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)又可表示为又可表示为(x1(n)分成偶序列和奇序列分成偶序列和奇序列)14/0)
15、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)也可以进行
16、相同的分解：也可以进行相同的分解：注意：通常我们会把注意：通常我们会把 写成写成 。kNW2/kNW2)()(62/5kXWkXKN 10)()()(10 NkWnxnxDFTkXNnnkN两点两点DFT及其蝶形图及其蝶形图120()()()01nknX kDFT x nx n Wk100010022220(0)()(0)(1)(0)(1)nnXx n WxWW xxW x110 11 1022220(1)()(0)(1)(0)(1)nnXx n WxWWxxW x-1-102(0)(0)(1)XxW x02(1)(0)(1)XxW x(0)x(1)x02WN点点DFT的第二次时域抽取分解图的
17、第二次时域抽取分解图(N=8)(0)(4)(2)(6)(1)(5)(3)(7)WN0WN0WN0W0N-1-1-1-1X(0)X(1)X(0)X(1)X(0)X(1)X(0)X(1)33445566WN0WN2WN0WN2-1-1-1-1X(0)X(1)X(2)X(3)X(0)X(1)X(2)X(3)11121222WWWWN0N1N2N3-1-1-1-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)xxxxxxxxN点点DITFFT运算流图运算流图(N=8)-1-1-1-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)xxxxxx(0)x(2)x(6)x(
18、1)x(5)x(3)x(7)x(4)000，000001，100010，010011，110100，001101，101110，011111，111X3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN0WN2WN0WN2X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3-1-1-1-1WN0WN0WN0WN0-1-1-1-13、DITFFT算法与直接计算算法与直接计算DFT运算量的比较运算量的比较22log2NNN 1)、N=2M的的DFT运算可分成运算可分成M级，每一级有级，每一级有N/2个蝶形个蝶形 ，每
19、个蝶形有一次复乘两次复加。，每个蝶形有一次复乘两次复加。NN2log2NN2log2)、所以、所以M级共有级共有 次复乘和次复乘和 次复加。次复加。3)、若直接计算、若直接计算DFT，需需N2次复乘和次复乘和N(N-1)次复加。次复加。显然，当显然，当N较大时，有：较大时，有：例如，例如，N=210=1024时时221048576204.8(/2)log5120NNNFFT算法与直接计算算法与直接计算DFT所需乘法次数的比较曲线所需乘法次数的比较曲线4、DITFFT的运算规律及编程思想的运算规律及编程思想 FFT的每级（列）计算都是由的每级（列）计算都是由N个复数数据（输入）两个复数数据（输入
20、）两两构成一个蝶型（共两构成一个蝶型（共N/2个蝶形）运算而得到另外个蝶形）运算而得到另外N个复数个复数数据（输出）。数据（输出）。当数据输入到存储器以后，每一组运算的结果，当数据输入到存储器以后，每一组运算的结果，仍然存仍然存放在这同一组存储器中放在这同一组存储器中直到最后输出。直到最后输出。例：将例：将x(0)放在单元放在单元A(0)中，将中，将x(4)放在单元放在单元A(1)中，中，W80 放在一个暂存器中。放在一个暂存器中。将将x(0)+W80 x(4)送回送回A(0)单元单元将将x(0)-W80 x(4)送回送回A(1)单元单元1）原位运算原位运算 (亦称同址计算亦称同址计算)-1-
21、13(0)X3(1)X(0)x(4)x08W回顾：回顾：N点点DITFFT运算流图运算流图(N=8)如上所述，如上所述，N点点DITFFT运算流图中，每级都有运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子个蝶形。每个蝶形都要乘以因子WNP，称其为称其为旋旋转因子转因子，p称为旋转因子的指数。称为旋转因子的指数。2)旋转因子的变化规律旋转因子的变化规律 观察观察FFT运算流图发现，第运算流图发现，第L级共有级共有2L-1个不同的个不同的旋转因子。旋转因子。N=23=8时的各级旋转因子表示如下：时的各级旋转因子表示如下：L=1时，时，WNp=WN/4J，N/4=21=2L，J=0L=2时，
22、时，WNp=WN/2J，N/2=22=2L，J=0，1L=3时，时，WNp=WNJ，N =23=2L，J=0，1，2，3对对N=2M的一般情况，第的一般情况，第L级的旋转因子为：级的旋转因子为：12,.,1,012 LJPNJWWLMLMLMLN 222212,.,1,0122 LJNJNPNJWWWLMMLLMJp 23)编程思想及流程图编程思想及流程图(自学自学)开始开始送入送入x(n)和和N=2M调整输入调整输入x(n)的顺序的顺序for(L=1;L=M;L+)B=2L-1for(J=0;J=B-1;J+)p=J2M-Lfor(k=J;k=N-1;k=k+2L)pNpNWBkXkXBkX
23、WBkXkXkX)()()()()()(输出结果输出结果结束结束4）码位倒序码位倒序 由由N=8蝶形图看出：原位计算时，蝶形图看出：原位计算时，FFT输出的输出的X(k)的次序正好是顺序排列的，即的次序正好是顺序排列的，即X(0)X(7),但输但输入入x(n)都不能按自然顺序存入到存储单元中，而是都不能按自然顺序存入到存储单元中，而是按按x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)的顺序存入的顺序存入存储单元存储单元，即为乱序输入，顺序输出。即为乱序输入，顺序输出。这种顺序看起来相当杂乱，然而它是这种顺序看起来相当杂乱，然而它是有规律有规律的。的。即即码位倒读
24、规则码位倒读规则。自然顺序自然顺序n二进制码表示二进制码表示码位倒读码位倒读码位倒置顺序码位倒置顺序n以以N=8为例：为例：0123456700000101001110010111011100010001011000110101111104261537看出：看出：码位倒读后码位倒读后的顺序刚好是数据送入计算机内的顺序。的顺序刚好是数据送入计算机内的顺序。四、按频域抽取的基四、按频域抽取的基2-2-FFTFFT算法算法 设序列设序列x(n)长度为长度为N=2M，首先将首先将x(n)前后对半前后对半分开，得到两个子序列，其分开，得到两个子序列，其DFT可表示为如下形式可表示为如下形式 1-NN/2
25、nknN1-N/20nknN1-N0nknNWnxWnxWnxX(k)()()(-1N/20n)2/(-1N/20n)2/()(NnkNknNWNnxWnx -1N/20n-1N/20n)2/()1()(knNkknNWNnxWnx1,.1,0 Nk kjkkNNjkNNeeW1222/式中，式中，12/01,.1,0,)2/(1)()(NnknNkNkWNnxnxkX:)(k分为两部分分为两部分的奇偶可把的奇偶可把按按kX12/.2,1,0,12,2 Nrrkrk及及令令为偶数时，为偶数时，K 12/02/12/02)2/()()2/()(2NnrnNNnrnNWNnxnxWNnxnxrX为
26、奇数时，为奇数时，K 12/02/12/012)2/()()2/()(12NnrnNnNNnnrNWWNnxnxWNnxnxrX）（1,.,1,0)2/()()()2/()()(221 NnNnWNnxnxnxNnxnxnx令令 nrNNnnrNNnWnxrXkXWnxrXkX2/12/0222/12/011)()12()()()2()(-1-1)2()(Nnxnx1,1,02NnnNWNnxnx)2()(nNW)2(Nnx)(nxDIFFFT一次分解运算流图一次分解运算流图(N=8)-1-1-1-1-1-1-1-1W WN N0 0W WN N1 1W WN N2 2W WN N3 3)0(
27、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两点两点DFT及其蝶形图及其蝶形图100010022220(0)()(0)(1)(0)(1)nnXx n WxWW xxW x110 11 1022220(1)()(0)(1)(0)(1)nnXx n WxWWxxW x-1-102(0)(0)(1)XxW x02(1)(0)(1)XxW x(0)x(1)x02W考虑到考虑到 上图可变形为上图可变形为:0021NWW-1-1(0)(0)(1)Xxx02(1)(0)(1)XxW x(0)x(1)x02W
28、DIFFFT运算流图运算流图(N=8)0NW1NW2NW3NW-1-1-1-1-1-1-1-1x(0)x(3)x(1)x(2)x(4)x(5)x(6)x(7)0NW2NW2NW0NWX(0)X(6)X(4)X(2)X(1)X(5)X(3)X(7)0NW0NW0NW0NW-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1 时间抽取算法与频率抽取算法的比较时间抽取算法与频率抽取算法的比较1)频率抽选法和时间抽选法总的计算量是相同的频率抽选法和时间抽选法总的计算量是相同的NN2log2NN2log复乘：复乘：复加：复加：2)频率抽取法和时间抽取法一样，都适用于频率抽取法和时间抽取法一样
29、，都适用于原位运原位运 算算，即蝶形的输入和输出占用同一个存储单元。即蝶形的输入和输出占用同一个存储单元。3)均存在码位倒序问题。均存在码位倒序问题。4)频率抽选法和时间抽选法一样，基本运算也是蝶形频率抽选法和时间抽选法一样，基本运算也是蝶形 运算。但两者的蝶形形式略有不同。运算。但两者的蝶形形式略有不同。-1-1)()()(21kXWkXkXkN)()()2(21kXWkXkNXkN)(1kX)(2kXkNW-1-1)2()(Nnxnx1,1,02NnnNWNnxnx)2()(nNW)2(Nnx)(nx五、五、IDFTIDFT的快速算法的快速算法IFFTIFFT上述上述FFT算法流图也可以用
30、于离散傅里叶逆变换算法流图也可以用于离散傅里叶逆变换(Inverse Discrete Fourier Transform，简称简称IDFT)。比较比较DFT和和IDFT的运算公式：的运算公式：10)(1)()(10 NnWkXNkXIDFTnxNknk10)()()(10 NkWnxnxDFTkXNnnkN1)旋转因子：旋转因子：kNkNWW 2)系数：系数：MNN21，如果希望直接调用如果希望直接调用FFT子程序计算子程序计算IFFT，则可用则可用下面的方法：下面的方法：10101()()1()()NknNkNknNkx nX k WNx nXk WN对上式两边同时取共轭，得：对上式两边同时取共轭，得：1011()()()NknNkx nX k WDFT X kNN
