范例-FFT算法程序及分析

《数字信号处理》—西南交通大学FFT算法程序及分析摘要：《FFT的算法程序分析》主要分析了按时间抽取（DIT）的快速傅立叶变换的基2FFT算法，通过对基2FFT算法的原理的分析及与DFT算法运算量的比较，进一步推导出了基rFFT算法，重点是基rFFT算法的推导。在具体的实例中，我们重点分析了FFT过程中幅值大小与FFT选用点数N的关系，验证FFT变换的可靠性，考察在FFT中数据样本的长度与DFT的点数对频谱图的影响。关键字：基2FFT算法，基rFFT算法，样本长度，选用点数要求：学习书上第六节的内容，自己编程实现FFT算法。给出典型信号的时域和频域图，并加以分析。可尝试实现分段卷积程序。论文内容含原程序、运行结果，理论分析和典型信号时域图。一．快速傅立叶变换（FFT）简介离散傅立叶变换（DFT）是信号分析与处理中的一种重要的变换。因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大。所以在快速傅立叶变换（FFT）出现以前，直接用DFT算法进行频谱分析和信号的实时处理是不切实际的。1965年，库利（J.W.Cooley）和图基（J.W.Tukey）在《计算数学》杂志上发表了“机器计算傅立叶级数的一种算法”的文章，这是一篇关于计算DFT的一种快速有效的计算方法的文章。它的思路建立在对DFT运算内在规律的认识之上。这篇文章的发表使DFT的计算量大大减少，并导致了许多计算方法的发现。这些算法统称为快速傅立叶变换（FastFourierTransform），简称FFT。1984年，法国的杜哈梅尔（P.Dohamel）和霍尔曼（H.Hollmann）提出的分裂基快速算法，使运算效率进一步提高。快速傅立叶变换（FFT）不是一种新的变换，而是离散傅立叶变换（DFT）的一种快速算法。FFT分成两大类，即按时间抽取（decimation-in-time，缩写为DIT）法和按频率抽取（decimation-in-frequency，缩写为DIF）法。FFT算法主要包括基2FFT算法，基4FFT算法，混合基FFT，基rFFT算法和分裂基FFT算法。二．快速傅立叶变换（FFT）定义设x(n)为N点有限长序列，其DFT为10)()(NnnkNWnxkXk=0,1,…,N-1（1）其反变换（IDFT）为《数字信号处理》—西南交通大学WnkNNokkXNnx1)(1)(n=0,1,…,N-1（2）二者的差别只在于WN的指数符号不同，以及差一个常数因乘子1/N。一般x(n)和WnkN为复数序列。对某一个k值，直接按（1）式计算X(k)值需要N次复数乘法、（N-1）次复数加法。因此，对所有N个k值，共需要N2次复数乘法及N（N-1）次复数加法运算。当N»1时，N（N-1）N2。（1）式可写为]Im[]Re[])(Im[)](Re[)()(10110]Re[)](Im[]Im[)])((Re[]Im[)](Im[]Re[)](Re[NnnkNnkNnkNnkN所以，一次复数乘法需要四次实数乘法和二次实数加法；一次复数加法则需二次实数加法。因而每运算一个X（K）需4N次实数乘法及2N+2（N-1）=2（2N-1）次实数加法。所以整个DFT运算共需要4N2次实数乘法和N*2（2N-1）=2N（2N-1）次实数加法。由上述可见，，N点DFT的乘法和加法运算次数均与N2成正比。当N较大时，运算量相当可观。所以，必须减少其运算量，才能使DFT在各种科学和工程计算中得到运用。而DFT运算时间能否减少，关键在于时间运算是否存在规律以及如何利用这些规律。仔细观察DFT的运算可以看出，利用系数WnkN的一些固有特性，就可以减少运算量：1)WnkN的对称性)(*WnkN=WnkN2)WnkN的周期性WnkN=WkNnN)(=WNknN)(3）WnkN的可约性WnkN=WmnkmNWnkN=WmnkmN//即有：WkNnN)(=WknNN)(=WnkNWNN2/=-1WNkN)2/(=—WkN因此，利用这些性质，可以合并DFT运算中的有些项；利用这些性质可以将长序列的DFT分解为短序列的DFT。从而，减少运算次数，真正做到提高运算效率。三.FFT基本形式1.基2FFT算法基2FFT算法可以分为按时间抽取（DIT）的基-2FFT算法（库利-图基算法）和按频率《数字信号处理》—西南交通大学抽取的基4FFT算法。我们具体分析按时间抽取（DIT）的基-2FFT算法。算法原理：先设序列点数为N=2L，L为整数。将N=2L的序列x(n)（n=0，1，…，N-1）先按n的奇偶分成两组：)()12()()2(21rrxrrxxx，r=0,1,…,12N则可以将一个N点DFT分解成两个N/2点的DFT，而x1(r)和x2(r)以及X1(k)和X2(k)都是N/2点的序列。X（k）表达为前后两部分：前半部分)()()(21kkkXXWXkNk=0,1,…,2N-1后半部分)2()2()2(2)2/(1NkNkNkXXWXNkN)()(21kkXWXkNk=0,1,…,2N-1只要求出0到（N/2-1）区间的所有X1(k)和X2(k)值，即可求出0到（N-1）区间内的所有X（k）的值。此时，我们可以利用蝶形信号流图符号表示DFT的运算。下图表示N=23=8的情况：)0()0(1xx)0(1X)0()0(1xx)1(1X)1(X)0()0(1xx)2(1X)2(X)0()0(1xx)3(1X)3(X)0()0(1xx)0(2XWN0)4(X)0()0(1xx)1(2XWN1-1)5(X)0()0(1xx)2(2XWN2-1)6(X)0()0(1xx)3(2XWN3-1)7(X-1因为N=2L，2N仍能被2整除。将X（k）分解后的两部分按奇偶各分解为两个4N点序列，而这两个序列又可再分解为两个4N点序列，依次类推，可以一直分解到只有两点的序列。由于N=2L,分解共需要L次。DFT与FFT运算量的比较：由按时间抽选法FFT的流图可见，当N=2L时，共有L级蝶形，每级都由N/2个蝶形运算组成，每个蝶形有一次复乘、二次复加，因此每级运算都需要N/2次复乘和N次复加，这样L级运算总共需要N/2点DFTN/2点DFT《数字信号处理》—西南交通大学复乘数NNLNmFlog222复加数NNNLaFlog2以乘法为例，直接DFT复数乘法次数是N2，FFT复数乘法次数是NNlog22。所以二者的计算量之比为NNNNLNNNloglog22222222.基rFFT算法设序列x(n)长度为N=rM,r和M均为任意整数。仿照基2DIT-FFT算法的推导方法，用M位r进制数表示k和nrkkkkiMiirMMrk10)(0211,2,1,0rki10)(021MiiirrnnnnMMn1,,2,1,0rni那么x（n）的DFT可写成如下形式：)()(021kkkMMXkX101010021011)(rrrpNMMnnnxMWnnn（1-1）式中nprkrniMiijMjjk1010为导出基rDIT-FFT算法，将n分解得WpN=WMiMmiinrrkN1011WMiMmiinrrkN1022WMiiinrkN100＝WnrkMMN110WiMmiinrrkN1022WMiiinrkN100将上式代入（1－1）式得)(021kkkMMX101010110rrrnnnM)(021nnnMMx《数字信号处理》—西南交通大学WnrkMMN110WiMmiinrrkN1022WMiiinrkN100（1-2）由上式同样可得到M个递推公式100132010121012110021010131021002102101010021221110rNMMMMMMMMrNMMrNMMMnnrknnrrknnrkxWkkkknxkkkkxkkkkXWknnnxkknnnxWnnnknnnxriiiMMMiiMMM（1-3）下面以N=42为例说明基rDIT-FFT算法。将M=2,n＝4代入（1-3）式得基4DIT-FFT递推公式kkxkkXWknxkkxWnnknxnnkknnkxNN01230400101230401001010001110(1-4)式中;3,2,1,0ki3,2,1,0ni（k1k0）和(n1n0)分别表示k和n的M(M=2)位4进制数。knx001表示4个4点DFT。300120nkkxknx001WnkN00Wnk014也表示4个4点DFT。第一级的4点DFT与第二级的4点DFT通过旋转因子WnkN00连接起来，这与前面讨论的基2DIT-FFT是一样的。（1-4）式也可以写成矩阵形式。x1(n0k0)的矩阵形式：《数字信号处理》—西南交通大学030201001111xxxx＝30201000xxxx＝jjjj11111111111130201000xxxx=Q30201000xxxx式中Q=jjjj111111111111Qxxxx13121101111131211101xxxxQxxxx23222102111132221202xxxxQxxxx33323103111133231303xxxxkkx012的矩阵形式：Qxxxx302010002222302010001111xxxx《数字信号处理》—西南交通大学Qxxxx312111012222WxWxWxWxNNNN3121110131211101Qxxxx322212022222WxWxWxWxNNNN6141210132221202Qxxxx332313032222WxWxWxWxNNNN9161310133231303完成4点DFT的矩阵乘Q[ABCD]T运算的蝶形结符号如图1（a）所示，对应实际运算流图如图1（b）所示。图（b）只需要8次复加，但需要一次倒序。由x1(n0k0)和x2(k1k0)的矩阵表示式容易画出16点基4DIT-FFT运算流图如图2所示。(a)图1（a）ABCDTQ的蝶形结符号(b)图1（b）ABCDTQ的实际运算流图《数字信号处理》—