快速傅里叶变换FFT算法源码经典.pdf

快速傅里叶变换FFT算法源码经典.pdf

《快速傅里叶变换FFT算法源码经典.pdf》由会员分享，可在线阅读，更多相关《快速傅里叶变换FFT算法源码经典.pdf（68页珍藏版）》请在冰点文库上搜索。

快速傅里叶变换FFT算法及其应用摘要本文较为系统地阐述了快速傅里叶变换的算法原理及其在数字信号处理等工程技术中的应用。

根据抽取方法的不同，一维基2FFT算法分为两种：

频域抽取的FFT算法和时频域抽取的FFT算法。

第1节阐述了这两种FFT算法的原理。

第2节给出了两种算法的编程思想和步骤。

第3节阐述了一维非基2FFT的两种算法：

Cooley-tukeyFFT算法和素因子算法（PrimeFactorAlgorithm）的思想原理，给出了在把一维非基2DFT的多层分解式转化为二层分解的过程中，如何综合运用这两种算法以达到总运算次数最少的方案；并以20点DFT为例描述了非基2FFT算法实现的一般步骤。

第4节介绍了一维FFT算法在计算连续时间信号的傅里叶变换、离散信号的线性卷积、离散信号压缩和滤波等数字信号处理中的典型应用。

第5节把一维FFT变换推广到二维FFT变换，并在一维FFT算法的基础上，给出了二维FFT算法的原理和实现过程。

最后在附录中给出了一维DFT的基2FFT算法（包括频域抽取的FFT和IFFT算法、时域抽取的FFT和IFFT算法），一维任意非基2FFT算法，二维DFT的基2FFT算法以及二维DFT的任意非基2FFT算法的详细的VisualC+程序。

本文通过各种流程图和表格，较为深入系统地阐述了FFT的算法原理；运用Matlab编程，通过大量生动的实例，图文并茂地列举出了FFT算法的各种应用，并在每个实例中都附上了完整的Matlab程序，可供读者参考。

由于篇幅所限，本文未涉及FFT变换以及其应用的数学理论背景知识。

关键词：

FFT算法的应用，一维基2FFT算法，频域抽取，时域抽取，非基2FFT算法，Cooley-Tukey算法，素因子算法，线形卷积，信号压缩和滤波，二维FFT算法快速傅里叶变换FFT算法及其应用1目录1一维DFT的快速算法FFT.11.1频域抽取的基2算法.11.1.1正变换的计算.11.1.2逆变换的计算.41.2时域抽取的基2算法.52一维基2FFT算法编程.63一维任意非基2FFT算法.103.1COOLEY-TUKEYFFT算法.103.2素因子算法（PRIMEFACTORALGORITHM,PFA）.113.3一维任意非基2FFT算法.134一维FFT算法的应用.164.1利用FFT计算连续时间信号的傅里叶变换.164.2利用FFT计算离散信号的线性卷积.194.3利用FFT进行离散信号压缩.214.4利用FFT对离散信号进行滤波.244.5利用FFT提取离散信号中的最强正弦分量.275二维DFT的快速变换算法及应用简介.325.1二维FFT变换及其算法介绍.325.2二维FFT变换算法的应用.33参考文献.33附录.341一维DFT的基2FFT算法VISUALC+程序.34

（1）频域抽取的FFT和IFFT算法.34

（2）时域抽取的FFT和IFFT算法.392一维任意非基2FFT算法VISUALC+程序.443二维DFT的基2FFT算法VISUALC+程序.494二维DFT的任意非基2FFT算法VISUALC+程序.5711一维DFT的快速算法FFT当序列fn的点数不超过N时，它的N点DFT定义为21001NiknNnFkfnekN=

（1）反变换IDFT定义为210101NiknNkfnFkenNN=

（2）二者形式相似，快速算法的原理一样，这里先就其正变换进行讨论。

令2/iNNWe=，当k依次取为0,1,2,1NL时，可表示为如下的方程组：

0001020

（1）1011121

（1）2021222

（1）

（1）0

（1）1

（1）

（1）0012110121201211011NNNNNNNNNNNNNNNNNNNNNNFfWfWfWfNWFfWfWfWfNWFfWfWfWfNWFNfWfWfNW=+=+=+=+LLLML（3）由上式可见，直接按照定义计算N点序列的N点DFT时，每行含N个复乘和N个加，从而直接按定义计算点的总计算量为2N个复乘和2N个加。

当N较大时，2N很大，计算量过大不仅耗时长，还会因字长有限而产生较大的误差，甚至造成计算结果不收敛。

所谓快速傅里叶变换就是能大大减少计算量而完成全部点计算的算法。

下面介绍两种经典的DFT的快速算法：

频域抽取的FFT算法和时域抽取的FFT算法。

1.1频域抽取的基2算法1.1.1正变换的计算这里仅介绍基2算法，即是2的整次幂的情况。

由定义1001NknNnFkfnWkN=（4）把fn分成两半，即fn和/2fnN+（0/21）nN，代入（4）式得/21/21（/2）00/201NNknknNNNnnFkfnWfnNWkN+=+（5）1一维DFT的快速算法FFT2由于（/2）/2

（1）knNknkNkknNNNNWWWW+=（5）式两项又可合并为/210

（1）/201NkknNnFkfnfnNWkN=+（6）当2kr=为偶数时，注意到

（1）1k=，222/knrnirnNNNWWe=/2rnNW=，（6）式变为/21/20/21/202（/2）（）（）0/21NrnNnNrnNnFrfnfnNWgnWGrrN=+=（7）当21kr=+为奇数时，（21）2（21）/2knrnirnNnrnNNNNWWeWW+=，（6）式变为/21/20/21/2021（/2）（）（）0/21NnrnNNnNrnNnFrfnfnNWWpnWPrrN=+=+=（8）这样就把一个N点序列（fn）的N点DFT（Fk）计算化成了两个/2N点序列（gn和pn）的/2N点DFT（Gr和Pr）计算。

由fn划分成gn和pn的计算量为N个加，即/2fnfnN+和/2,0/21fnfnNnN+和/2N个乘，即（/2）,0/21nNfnfnNWnN+由于gn算出的/2N点DFT，是fn的N点DFT（Fk）中k为偶数的那一半，由pn算出的则是k为偶数的那一半，故需要把偶数k的Fk抽出来放在一起作为gn的DFT（（）Gr）输出，同时把奇数k的Fk抽出来放在一起作为pn的DFT（（）Pr）输出。

由于k是频域变量，故这种算法称为频域抽取的FFT算法。

接着，两个/2N点DFT仍可用上述方法各经/4N个乘/2N个加划分成两个/4N点DFT（同时还要做相应的频域抽取），从而共划分成4个/2N点DFT，总划分计算量仍是N个加和/2N个乘。

这样的划分可一步步做下去，不难看出，快速傅里叶变换FFT算法及其应用3每步的总划分计算量都是N个加和/2N个乘。

经过1M步的划分后就划成了/2N个如下两点DFT的计算问题00012210110222（）AaWbWabBaWbWabW=+=+=+=（9）上式所需计算量是2个加和1个乘，于是完成/2N个两点DFT的总计算量仍是N个加和/2N个乘。

从而完成全部N点DFT的总计算量2logMNNN=个加和2/2（/2）logMNNN=个乘，这比直接按定义计算所需的2N个乘和加要少得多。

例如，1021024N=，10M=，用FFT算法计算所需的乘法个数为/2MN51024=，而直接按定义计算所需的乘法个数为210241024N=，二者相差1024/5200倍。

若直接计算需半小时，而用FFT计算只需9s即可完成，可见其效率之高，而且N越大，FFT的效率提高越明显。

图图1频域抽取的频域抽取的8点点FFT计算流图计算流图一般情况下，由于做了1M次分奇偶的抽取，此算法最后的/2N个两点DFT计算出的Fk不是顺序抽取的。

次序的变化可用二进码来说明：

第一次抽取所分的奇偶是由二进码第1位是1或0来区分的，该位为0时为偶数，该位为1时为奇数，第二次抽奇偶是由二进码第2位是1或0来区分的，每次抽取都是把偶数项放在前（左）边，把奇数项放在后（右）边，从而抽取以后数的二进码是按照二进制位从左向右依次排列的，和普通二进制数从右向左依次排列的的规律正好相反，所以称为倒位序。

在计算出Fk之后要把倒位序变成顺序。

f0F0000F0F0000f1-1W20F4100F2F1001gnf2-1W40F2010F4F2010f3-1W41-1W20F6110F6F3011f4-1W80F1001F1F4100f5-1W81-1W20F5101F3F5101pnf6-1W82-1W40F3011F5F6110f7-1W83-1W41-1W20F7111F7F71111一维DFT的快速算法FFT41.1.2逆变换的计算所谓逆变换是指由Fk求fn的计算，若直接按定义10101NknNkfnFkWnNN=做计算，则除了求和号和正变换相同的计算量外，每算一个fn都还需再多做一个乘1/N的乘法运算。

故按定义完成全部N点DFT的总计算量是2N个加和

（1）NN+个乘。

下面从图导出它的快速算法，先讨论第3列的2点DFT的逆运算如何完成。

由式（7）得，02（）abAabWB+=由上式不难解出02021（）21（）2aABWbABW=+=（10）图图2频域抽取的频域抽取的8点点IFFT计算流图计算流图此计算过程如图2所示，可以看出：

左边各列的划分计算也都是由/2N个碟形运算来完成的，只是各碟形运算所乘的相移因子W不同。

把每个碟形运算都用图的办法变成对应的逆运算，并把它们按输入在左、输出在右重新排列，就F0000F0F00001/8f0F1001F2F41001/8W2-0-1f1F2010F4F20101/8W4-0-1f2F3011F6F61101/8W2-0-1W4-1-1f3F4100F1F10011/8W8-0-1f4F5101F3F51011/8W2-0-1W8-1-1f5F6110F5F30111/8W4-0-1W8-2-1f6F7111F7F71111/8W2-0-1W4-1-1W8-3-1f7快速傅里叶变换FFT算法及其应用5得到了全部N点IDFT的计算流图。

给出了8N=的示例，图中先对顺序输入的Fk做1M次的频域抽取，并把3个乘1/2的运算合成了一个乘1/8的运算放在了最前边，然后就开始做求逆的碟形运算。

1.2时域抽取的基2算法比较正变换DFT和反变换IDFT的定义式1001NknNnFkfnWkN=10101NknNkfnFkWnNN=可见，去掉乘1/N的运算，把1W换成W，交换Fk、fn和k、n，反变换定义式就变成了正变换的定义式。

对图2做这些变换，则得到图3的流程图。

对图1做这些变换，则得到图4的流程图。

这就是时域抽取的算法流图。

进行碟形运算之前，先要对顺序的时域输入序列进行1M次的奇偶抽取，故称为时域抽取的FFT算法。

图图3时域抽取的时域抽取的8点点FFT计算流图计算流图比较图2和图3不难看出，两种算法的计算量是完全一样的。

这里先算出/2N个两点的DFT000122（）（/2）（）（/2）AfnWfnNWfnfnN=+=+f0000f0f0000F0f1001f2f4100W20-1F1f2010f4f2010W40-1F2f3011f6f6110W20-1W41-1F3f4100f1f1001W80-1F4f5101f3f5101W20-1W81-1F5f6110f5f3011W40-1W82-1F6f7111f7f7111W20-1W41-1W83-1F72一维基2FFT算法编程6101122（）（/2）（）（/2）BfnWfnNWfnfnN=+=+（11）图图4时域抽取的时域抽取的8点点IFFT计算流图计算流图然后把/2N个两点的DFT组合成/4N个4点的DFT，再把/4N个4点的DFT组合成/8N个8点的DFT，经过1M次的组合之后，就得到了顺序N点DFT计算结果。

算法原理参考文献【4】。

2一维基2FFT算法编程以频域抽取的基2FFT正变换为例，对FFT的信号流图进行讨论，以找到FFT算法的规律。

1）分级在进行DFT变换的过程中，从N点DFT到两点DFT共分了M级，如图1所示，从左到右依次为1m=级，2m=级，mM=级。

2）倒位序在频域抽取的基2FFT算法中，输出数据不是按照序列的先后顺序排列的，这是由于变换过程中，输出按奇、偶抽取的缘故。

如果将序列xn中标号n用二进制值0212（）MMnnnL表示，那么在FFT信号流图输入端，xn位于f01/8F0000F0F0000f11/8-1W2-0F4100F2F1001f21/8-1W4-0F2010F4F2010f31/8-1W4-1-1W2-0F6110F6F3011f41/8-1W8-0F1001F1F4100f51/8-1W8-1-1W2-0F5101F3F5101f61/8-1W8-2-1W4-0F3011F5F6110f71/8-1W8-3-1W4-1-1W2-0F7111F7F7111快速傅里叶变换FFT算法及其应用71202（）MMnnnL处，称为倒序。

以8点FFT为例，顺序和倒序的关系如表1所示。

表表1顺序和倒序对照表顺序和倒序对照表顺序倒序十进制数二进制数二进制数十进制数0123456700000101001110010111011100010001011000110101111104261537从表1可以看出，一个自然顺序二进制数，是在最低位加1，逢2向左移位；而倒序数的顺序是在最高位加1，逢2向右移位。

用i表示顺序数，j表示倒序数，k表示位权重。

对于一个倒序数j来说，下一个倒序数可以按下面的方法求：

先对最高位加1，相当于十进制运算/2jN+。

如果/2jN，说明二进制最高位为0，则直接由/2jN+得到下一个倒序值；如果/2jN，说明二进制最高位为1，则将j的最高位变为0，通过/2jjN实现，同时令/2kk，接着判断次高位是否为0，直到位为0时，令jjk=+LL。

归结起来算法流程图如图5所示。

2一维基2FFT算法编程8图图5倒位序程序流程图倒位序程序流程图图图6频域抽取频域抽取FFT程序流程图程序流程图3）蝶形运算单元和同址计算频域抽取的信号流图中，基本的运算结构如图7所示，该运算结构的形状像蝴蝶，故称为“蝶形运算单元”。

在这一结构中，DFT和IDFT运算关系分别为1111（）mmmrmmmNFpFpFqFpFpFqW=+=1111（）/2（）/2rmmmNrmmmNfpfpfqWfpfpfqW=+=（12）输入N,信号fNM=log2Nl=1tolbr=（l-1）2m-1n=（l-1）:

la:

N-2m=1toMla=2M+1-mlb=la/2lc=n+lbt=fn+flcflc=（fn-flc）WNrfn=t倒位序重排信号j=N/2i=1toN-2t=fifj=fifi=fjk=N/2j=j-kk=k/2kjj=j+kij是是否否输出快速傅里叶变换FFT算法及其应用9（a）两点DFT的计算（b）两点IDFT的计算图图7频域抽取频域抽取FFT算法流图中第算法流图中第m级碟形单元级碟形单元而在时域抽取的信号流图中，基本的运算结构如图8所示。

在这一结构中，DFT和IDFT运算关系分别为1111rmmmNrmmmNFpFpFqWFpFpFqW=+=1111（）/2（）/2mmmrmmmNfpfpfqfpfpfqW=+=（13）（a）两点DFT的计算（b）两点IDFT的计算图图8时域抽取时域抽取FFT算法流图中第算法流图中第m级碟形单元级碟形单元其中，p、q分别表示该蝶形运算单元的上下节点的序号。

可以看出参与运算的输入序号p，输出序号仍为q，并且该运算不涉及到其它的点，因此我们可以将输出的结果仍然放在数组中，称这样的操作为同址计算。

也就是说，共同占有同一个存储单元。

4）寻址和相移因子rNW的计算时域抽取基2FFT信号流图中，每一级有个蝶形单元。

每一级的个蝶形单元又可以分为若干组，每一组具有相同的结构和因子的分布。

如图1所示，第1级分为1组，第2级分为2组，第m级分为12m组。

fm-1p1/2fmpfm-1q-1WN-r/2fmqFm-1pFmpFm-1q-1WNrFmqfm-1p1/2fmpfm-1qWN-r-11/2fmqFm-1pFmpFm-1qWNr-1Fmq3一维任意非基2FFT算法10在第m级中，相邻组之间的间距（也即每个分组所含节点数）为12Mm+，每个蝶形单元的上下节点之间的距离（也即每个分组所含碟形单元数）为2Mm。

每组的相移因子22cos（）sin（）rNWrirNN=，其中m-1r=（l-1）2,0,1,l=,21MmL综合以上各步骤，得到频域抽取FFT程序流程图如图6所示。

采用类似的步骤可得到频域抽取IFFT流程图、时域抽取FFT与IFFT流程图（略）。

频域抽取FFT算法、时域抽取FFT算法的VisualC+源程序分别见附录1.

（1），1.

（2）。

在Matlab中，傅里叶变换及其逆变换分别用dwt（）和idwt（）函数实现。

3一维任意非基2FFT算法3.1Cooley-TukeyFFT算法FFT的核心是将一层运算映射为二层运算。

作N点FFT时，若N不是素数，则N可分解为12NNN=，那么由fn的DFT1001NnkNnFkfnWkN=（14）通过映射：

2121122112112201,0101,01nNnnnNnNkkNkkNkN=+=+（15）可得到212112211121221122（）（）（）NnnkNkNnkNNnknkNnknkNNNWWW+=而12NNN=，12NNNWW=，21NNNWW=，可化简为11212212nknknknkNNNNWWWW=（16）从而式（14）转化为21222111212111121200,（,）NNnknknkNNNnnFkkWWfnnW=（17）其中112201,01kNkN。

以20点FFT为例，1220,5,4NNN=，映射方式为：

124nnn=+，快速傅里叶变换FFT算法及其应用11125kkk=+，则计算流图如图9所示。

3.2素因子算法（PrimeFactorAlgorithm,PFA）当Cooley-TukeyFFT算法中的1N、2N互素的话，相移因子21nkNW可通过选择12,nn12,kk前的特殊系数而消去，FFT的算法进一步的简化。

12112212112201,0101,01nAnBnnNnNkCkDkkNkN=+=+（18）当ABCD、满足以下条件0mod0modADNBCN（19）图图9Cooley-Tukey20点点FFT算法算法n1k2f00W00F0f41W01F5f82W02F10f123W03F15f164W00F1f10W11F6f51W22F11f92W33F16f133f174W00F2W21F7f20W42F12f61W63F17f102f143W00F3f184W31F8W62F13f30W93F18f71f112W00F4f153W41F9f194W82F14W123F19k1=0n2=0n2=1n2=2n2=3k1=1k1=2k1=3k1=43一维任意非基2FFT算法1221modmodACNNBDNN（19）时，有121211122122211122112212（）（）（）AnBnCkDknkNNACnkADnkBCnkBDnkNNnkNnkNnknkNNWWWWWW+=（20）于是式（14）化简为212211212111121200,（,）NNnknkNNnnFkkWfnnW=（21）从而消去了相移因子21nkNW。

同样以20点FFT为例，修改映射方式为：

20,N=125,4NN=121245,04,03nnnnn=+（22）12121625,04,03kkkkk=+（23）则由式（22）得到的映射关系如表2，由式（23）得到的映射关系如表3，计算流图如图10所示。

表表2由式由式（22）得到的映射关系得到的映射关系n1n2012300510151491419281318331217274161611表表3由式由式（23）得到的映射关系得到的映射关系k1k2012300510151161611212172738131834491419nk快速傅里叶变换FFT算法及其应用133.3一维任意非基2FFT算法对于非素数N点DFT，对N做因式分解12lNNNN=L（24）令22llNNN=L，则12lNNN=。

于是式（24）中多层FFT转化为二层FFT。

如果1N与2lN互素，那么采用PFA算法进行映射，否则采用Cooley-TukeyFFT算法（此时需乘以相移因子）。

22llNNN=L可采用同样的方法分解成2N个3lN点DFT，其中33llNNN=L，依此类推，直到DFT长度为lN。

可以证明，复乘的总次数不大于12（）lNNNNl+L（25）例如，若63337N=，则复乘总次数至多为63（3373）630=。

而用基2的FFT算法计算64点DFT，需要的复乘总次数为192。

这说明，将N分解得越细，运算量越少。

实际中，常常将输入序列补零，使得N成为2的幂次数，这样能够减少复乘运算的次数。

3一维任意非基2FFT算法14图图10素因子素因子（PFA）20点点FFT算法算法再以20点FFT为例，进行如下三层因式分解123522NNNN=即15N=，23224N=，由于5与4互素，所以第一层采用PFA算法，相应的二层映射为12312345,04,03nnnnn=+1231231625,04,03kkkkk=+由于2与2不互素，所以第二层采用Cooley-TukeyFF