快速傅里叶变换.

N/2点DFTN/2点DFTx(0)x(1)x(3)x(2)x(4)x(5)x(6)x(7)X1(0)WN1X1(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)WN2WN3WN02第7章快速傅里叶变换7.1引言7.2直接计算DFT的问题及改进的途径7.3.1按时间抽取（DIT）的基-2FFT算法7.3.2按频率抽取（DIF）的基-2FFT算法补充内容：IDFT的高效算法补充内容：实序列的FFT算法37.1引言有限长序列通过离散傅里叶变换（DFT）将其频域离散化成有限长序列，但其计算量太大，很难实时处理，因此引出了快速傅里叶变换（FFT）。FFT并不是一种新的变换形式，它只是DFT的一种快速算法，并且根据对序列分解与选取方法的不同产生了多种算法。FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。447.2直接计算DFT的问题及改进的途径7.2.1直接计算DFT的运算量问题若x(n)是N点序列，实现x(n)的DFT，即求出N个X(k)，需要N2次复数乘法，N(N-1)次复数加法。10)()()(NnknNDFTWnxkXnx1,1,0Nk10)(1)()(NkknNIDFTWkXNnxkX1,1,0Nn55一个复数乘法包括4个实数乘法和2个实数加法。一次复数加法需二次实数加法。例如：x(n)N=1024，N2=1048576])}Re[)](Im[]Im[)]((Re[][Im)](Im[]Re[)]({Re[]}Im[])]}{Re[(Im[)]({Re[)()(101010nkNnkNnkNNnnkNNnnkNnkNNnnkNWnxWnxjWnxWnxWjWnxjnxWnxkX因而每运算一个X(k)需4N次实数乘法和2N+2(N-1)=2(2N-1)次实数加法。所以，整个DFT运算总共需要4N2次实数乘法和2N(2N-1)次实数加法。(a+jb)(c+jd)=(ac-bd)+j(bc+ad)6解决耗时的乘法问题是将数字信号处理理论用于实际的关键问题。特别是40年前，计算机的速度相当慢。因此，很多学者对解决DFT的快速计算问题产生了极大的兴趣。CooleyJW,TukeyJW.AnalgorithmforthemachinecomputationofcomplexFourierseries.MathematicsofComputation,1965,pp297~30177.2.2改善途径FFT的思路：10)()(NnknNWnxkX1,1,0Nk如何充分利用这些关系1.10W1,1.2mNNWWrrNWW.31.42/NWrrNWW2/.5对称性周期性8以四点DFT为例11111111(0)11(1)(2)1111(3)11xWWxxxWW0000012302460369(0)(0)(1)(1)(2)(2)(3)(3)304)()(nknWnxkX3,2,1,0k11111111(0)11(2)(1)1111(3)11xWWxxxWW911(0)[(0)(2)][(1)(3)](1)[(0)(2)][(1)(3)](2)[(0)(2)][(1)(3)](3)[(0)(2)][(1)(3)]XxxxxXxxxxWXxxxxXxxxxW11(0)(2)xx(1)(3)xx(0)(1)XX(2)(3)XX111W(0)(2)xx(0)(2)xx(1)(3)xx(1)(3)xx10利用WNnk的对称性和周期性，将大点数的DFT分解成若干个小点数的DFT，FFT正是基于这个基本思路发展起来的。分类：按时间抽取（DIT）算法和按频率抽取（DIF）算法。问题是如何分最有效？N点DFTN/2点DFTN/4点DFT2点DFT1个2个4个N/2个DecimationinTimeDecimationinFrequency11knNW的周期性knNW的对称性nkNnkNWW*)(nkNNnkNWW)2/()()(NknNkNnNnkN127.3.1按时间抽取（DIT）的基-2FFT算法1.算法原理设序列x(n)长度为N＝2M，M为整数，将x(n)（n=0,1…N-1）按n的奇偶分解成两组N/2点的子序列x1(r)=x(2r)x2(r)=x(2r+1)r=0,1…N/2-1（长度为N/2）则10)()(NnnkNWnxkX为奇数为偶数nnkNnnkNWnxWnx)()(12/0)12(12/02)12()2(NrkrNNrrkNWrxWrx12/02/12/02/)12()2(NrrkNkNNrrkNWrxWWrx1312/02/12/02/)12()2()(NrrkNkNNrrkNWrxWWrxkX)(1kX)(2kX12/,...,1,0)()()2/(21NkkXWkXNkXkN12/,...,1,0)()()(21NkkXWkXkXkNWNk-1X1(k)X2(k)X(k)X(k+N/2)X1(k)X2(k)X(k)X(k+N/2)WNkX1(k)、X2(k)都是N/2点的DFT，它们各自又可分成N/4点的DFT，如此继续分下去，直至两点DFT。两点DFT不需要乘法运算：N/2点DFTN/2点DFTx(0)x(1)x(3)x(2)x(4)x(5)x(6)x(7)X1(0)WN1X1(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)WN2WN3WN0)1()0()1()1()0()0(xxXxxXWN0x(0)x(1)X(0)X(1)2点DFT的运算流图N/4点DFTX1(0)X1(1)X1(2)X1(3)x(0)X3(0)X3(1)X4(0)X4(1)x(4)x(2)x(6)N/4点DFTWN/20WN/21N/4点DFTX2(0)X2(1)X2(2)X2(3)x(1)X5(0)X5(1)X6(0)X6(1)x(5)x(3)x(7)N/4点DFTWN/20WN/2115由于每一步分解都是按输入序列在时间上的奇偶次序，分解成两个半长的子序列，所以称“按时间抽取算法”。WN0WN1WN2WN3WN0WN2WN0WN2WN0WN0WN0WN0x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)第一级第二级第三级蝶形运算单元162.运算量比较对N=2M，共可分M次，即m=0,1…M-1，每一级有N/2个如下的蝶形单元一个蝶形单元只需一次复数乘法，两次复数加法。总共复数乘法：复数加法：xm(p)xm(q)WNrxm+1(p)xm+1(q))()()()()()(11qXWpXqXqXWpXpXmrNmmmrNmmNNMN2log22NNMN2log可以共享17直接计算DFT与FFT算法的计算量之比为NNNMNN222log22183.算法特点1）原位运算两点构成一个蝶形单元，并且这两点不再参与别的蝶形单元的运算。所谓原位运算，就是当数据输入到存储器后，每一级运算的结果仍然贮存在这同一组存储器中，直到最后输出，中间无需其它存储器。Xm-1(k)WNrXm-1(j)Xm(k)Xm(j)存储存储FFT蝶形运算Xm-1(k)Xm-1(j)Xm(k)Xm(j)192）旋转因子的变化规律(4.2.13)如上所述，N点DIT―FFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子，称其为旋转因子，p称为旋转因子的指数。不难发现，第L级共有2L-1个不同的旋转因子。N=23=8时的各级旋转因子表示如下：pNW3,2,1,0,,31,0,,20,,1222/24/J一般情况，第L级的旋转因子指数为L为从左到右的运算级数，编程时L为循环变量。LMJp2MLJL,...,2,1,12,...,2,1,0120如果蝶形运算的两个输入数据相距B个点，应用原位运算，则pNLLLpNLLLWBJXJXBJXWBJXJXJX)()()()()()(1111式中LMJp2MLJL,...,2,1,12,...,2,1,01B——两蝶形节点的“距离”（蝶距）12LB3）蝶形运算214）倒位序输入序列x(n)的排列次序不符合自然顺序。此现象是由于按n的奇偶分组进行DFT运算而造成的，这种排列方式称为“倒位序”。所谓“倒位序”，是指按二进制表示的数字首尾位置颠倒，重新按十进制读数。自然顺序二进制顺序码位倒置码位倒读顺序000000001001100420100102301111064100001151011015611001137111111722倒位序的实现在实际运算时，先按自然顺序将输入序列存入存储单元，再通过变址运算将自然顺序变换成按DIT-FFT算法要求的顺序。变址的过程可以用程序安排加以实现--称为“整序”或“重排”（采用码位倒读）注意:(1)当I=J时，数据不必调换；(2)当I≠J时，必须将原来存放数据x(I)送入暂存器R，再将x(J)送入x(I)，R中内容送x(J)，进行数据对调。(3)为了避免再次调换已调换过的数据，保证调换只进行一次，否则又变回原状。JI（或IJ)时，调换。23倒位序程序框图24DIT―FFT运算和程序框图254.小结☆全部计算分解为M级，或称为M次迭代。☆输入倒序，输出正序。☆每级都包含N/2个蝶形单元。☆每级的若干蝶形单元组成“群”。第1级群数为N/2，第2级群数为N/4，……第i级群数为N/2i，最后一级的群数为1。☆每个蝶形单元都包含乘WNk与-WNk的运算（可简化为乘WNk与加、减法各一次）。☆同一级中，各个群的W分布规律完全相同。WN0x(0)x(1)X(0)X(1)267.3.2按频率抽取（DIF）的基-2FFT算法把输出序列X(k)按k奇偶分解成越来越短的序列，称为按频率抽取FFT算法。1．算法原理设序列x(n)长度为N＝2M，M为整数。在把X(k)按k的奇偶分组之前，先把x(n)按n的顺序分解为前后两部分。12/12/010)()()()(NNnnkNNnnkNNnnkNWnxWnxWnxkX12/0)2/(12/0)2/()(NnkNnNNnnkNWNnxWnx12/0)2/()]2/()([NnnkNkNNWNnxWnxk)1(27现在对频率序列抽取，按k是奇数还是偶数将X(k)分成两部分。频率序列{X(2r)}是时间序列{x(n)+x(n+N/2)}的N/2点DFT，频率序列{X(2r+1)}是时间序列{[x(n)-x(n+N/2)]WNn}的N/2点DFT。这样，将N点DFT分解成两个N/2点DFT运算。12/01-N,1,0,k,)]2/()1()([)(NnnkNkWNnxnxkX12/02)]2/()([)2(NnnrNWNnxnxrX12/02