详解快速傅里叶变换FFT算法

详解快速傅里叶变换FFT算法快速傅里叶变换FFT是离散傅里叶变换DFT的一种快速算法，只有FFT才能在现实中有实际应用的意义。虽然许多学过数字信号处理这门课的同学都知道DFT和FFT，但实际上真正理解其算法原理的屈指可数，绝大部分同学知其然而不知其所以然，况且限于高校课程教学体制，课堂上不可能把这些原理和算法讲得明明白白的。为此，特意以本文讲解FFT算法的原理与实际应用，给欲往电子信息类专业进修和发展的同学一些课外参考。N点有限长序列x(n)的DFT为10()()WNnkNnXkxn其中2jnknkNNWe其逆变换IDFT为101x(n)X(k)WNnkNnN正逆变换的运算量都是相同的。x(n)和X（k）都是复数序列，计算一个X(k)值，需要N次复数乘法和N-1次复数加法。X(k)有N个点，所以总共需要N*N次复数乘法和N(N-1)次复数加法。复数运算实际上是通过实数运算来完成的。上式可以写成：10()Re()Im()ReImNnknkNNnXkxnjxnWjW10Re()ReIm()ImRe()ImIm()ReNnknknknkNNNNnxnWxnWjxnWxnW由此可见，一次复数乘法需要4次实数乘法和2次实数加减法。一次复数加法需要2次实数加法。所以每一个X（k）计算需要4N次实数乘法以及2N+2(N-1)=2(2N-1)次实数加法。整个DFT运算总共需要4N*N次实数乘法和N*2(2N-1)=2N(2N-1)次实数加法。当N足够大，N1时，直接计算DFT的乘法次数和加法次数都是和N的平方成正比。当N=1024时，DFT的运算量为1048576次，即一百多万次复乘运算，一块嵌入式32位处理器的最高速度为105百万指令每秒，那么它要完全计算这个DFT的时间最快也要1秒，期间还是独占CPU所有运算资源且不能有任何其他的中断请求。这样计算量太庞大，计算速递太慢了，谈不上实时性，根本没有实用意义。所以，我们就要利用DFT的系数的固有特性来简化计算，减少运算量。特性如下：1.共轭对称性：()*nknkNNWW2.周期性：()()nknNknkNNNN3.可约性：//nknmknkmNNNm得出：()()nNkNnknkNNN/21NNW(/2)kNkNNWW利用上述这些系数性质就可以合并DFT某些项的计算从而减少计算量。下面仅给出按时间抽选DIT的基-2FFT算法，也就是库利-图基算法的思路，其详细推导过程教科书上讲得更明白。设点数N=2^m，N=8，16，32，64„„1024，2048先将x(n)序列按n的奇偶性分成奇偶两组，每组的点数为N/2-1。偶：x(2r)=x1(r)奇：x(2r+1)=x2(r)则有1/21/212212000()()W(r)(W)W(r)(W)rkrkNNNnkkNNNNnrrXkxnxx利用系数W的可约性：2222/2/2jjNNNNWeeW上式可以表示为/21/211/22/21200()(r)WW(r)W()W()NNrkkrkkNNNNrrXkxxXkXk一个N点DFT分解成两个N/2点的DFT，此时X(k)只计算了前N/2个，而后N/2个的计算需要应用系数的周期性(/2)/2/2rkrkNNNWW,根据欧拉公式推导即得：11()(k)2NXkX22()(k)2NXkX22NNkkkNNNN于是得出前半部分和后半部分分别为：12()(k)W(k)kNXkXX12()(k)W(k)2kNNXkXX只需求0~N/2-1部分的X1和X2，即可求出0~N-1的所有X(k)。根据这两个式子，可以画出蝶形信号流图：X1(k)112(k)W(k)kNXXX2(k)kNW-112(k)W(k)kNXX每个蝶形运算需要一次复数乘法和两次复数加减法。当分解成两个N/2点DFT的蝶形运算后，一个N/2点需要(N/2)^2=N^2/4次复数乘法和N/2(N/2-1)次复数加法。两个N/2点DFT需要2(N/2)^2=N^2/2次复数乘法和N(N/2-1)次复数加法。把两个N/2点DFT合成N点DFT有N/2个蝶形运算，则还需要N/2个复数乘法和N个复数加法。总共需要N^2/2+N/2=N(N+1)/2，约等于N^2/2次复数乘法，N(N/2-1)+N=N^2/2次复数加法。对比前面所说的直接DFT的运算量便知此次减少了多少运算量。再将N/2分解为两个N/4点DFT，那么8点DFT则可分解成四个N/4=2点DFT。利用四个N/4点的DFT及两级蝶形组合运算来计算N点DFT，比只用一次分解的组合方式的计算量又减少了一半。DFT的运算量集中在N点DFT计算部分，尽力把多点直接DFT分解成少点直接DFT，尽量减少直接DFT的点数，就可以压低整个DFT的运算量。在把原序列多次逐级分解为奇偶两组的过程中，其序列标号是如何变化的呢？观察下面便知：(100)(101)(110)(111)原序列：x0,x1,x2,x3,x4,x5,x6,x7,(000)(001)(010)(011)因此可得：偶:x0(000),x4(100)偶原序列奇:x2(010),x6(110)偶:x1(001),x5(101)奇奇:x3(011),x7(111)3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X3(0)x3(1)x4(0)x4(1)x5(0)x5(1)x6(0)x6(1)x上图中无下标的x(n)和X(k)分别代表原序列和转换后的序列。有下标的x(n)和X(k)分别指的是原序列的分组和转换后的分组，下标相同代表同属一组，括号里的序号为点数。A(0)-A(7)是指内存中的储存单元。把原序列x(0),x(1)„„x(7)按奇偶分组后得到x3(n),x4(n),x5(n),x6(n)四组，则每组有N/4=2点，即其中的n=0,1。用前面分解N/2点的方法再分解一次便得：1433/40(k)()NlkNlXxlW1444/40(k)()NlkNlXxlW1455/40(k)()NlkNlXxlW1466/40(k)()NlkNlXxlW把上面四式具体化表示，仅列出X4(k)作为例子：11444/44/400(k)()NlklkNNllXxlWxW只有两点k=0,100044242(0)(0)(1)(2)(6)(2)(6)NXxWxxWxxWx11044242(1)(0)(1)(2)(6)(2)(6)NXxWxxWxxWx上式中2110221jjNWeeW，故上式不需要乘法，只有加减。同样按这个方法可以求出其他分组，而且这些两点DFT都可用一个蝶形结表示。下图是直接DFT与FFT算法的运算量的差距：蝶形运算的一般规律：原位运算细心观察8点DIF-FFT运算流程图可以归纳出蝶形运算的一般规律，其每级每列计算都是由N/2个蝶形运算构成，每一个蝶形运算如下式所示：11()()()rmmmNXkXkXjW11(j)()()rmmmNXXkXjW由一次复数乘法和两次复数加减法组成。某一列的任何两个节点k和j的节点变量进行蝶形运算后，得到结构为下一列k，j两节点的节点变量，与其他节点无关，所以可以采用原位运算。计算完一级后立刻把结果存进原来的内存单元，再进行下一级的运算。那么输入到结果输出其内存单元不变，如N=1024，则整个FFT只需要1024个长整型单元即可。这样就可以大量节省储存单元，方便单片机或微处理器低成本运行。输入序列按标号的奇偶的不断分组造成倒位序。DIT-FFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种倒序是很有规律的。由于N=2M，因此顺序数可用M位二进制数(nM－1nM－2„n1n0)表示。M次偶奇时域抽选过程如下图所示。细心观察8点DIF-FFT运算流程图还可以看出每蝶形的两个节点距离为2^(m-1)。一级，两节点为1二级，两节点为2三级，两节点为4根据上述规律可将蝶形运算的一般等式写成通式：（m是第m级运算）11()()()rmmmNXkXkXjW111()()(k2)mrmmmNXkXkXW11(j)()()rmmmNXXkXjW1111(2)()(k2)mmrmmmNXkXkXW求旋转因子W中r的算法：将上式中k换成L位二进制数（N=2^L），再乘以2^(L-m)，即将L位二进制数左移L-m位，右边补零，即可求出r。所以最终所需要的储存单元有N个+旋转因子W的N/2个=N+N/2个储存单元。要说明的是，按时域抽取法DIT与按频域抽取法DIF的基本蝶形是互为转置的，所以两种方法的本质是一样的，只是表现形式不同而已。转置就是将流图的所有支路方向都反相，并且交换输入与输出，但节点变量值不交换。在此就不详说DIF法了。再说明一下离散傅里叶逆变换IDFT的算法IFFT。只要对IDFT公式取共轭：101()IDFT()(k)WNnkNkxnXkXN1**01()(k)WNnkNkxnXN得到：*1***011()(k)W()NnkNkxnXDFTXkNN此式表明只要先将X(k)取共轭，就可以直接利用FFT子程序，最后再将运算结果取一次共轭，并乘以1/N即得到x(n)值，因此FFT和IFFT可以共用一个程序。倒位序算法：观察上图可知，倒序二进制数是从最高位加一并向次高位进位的，所以程序要实现这一“反向进位加法”的功能。算法思路是，若顺序二进制数等于倒序二进制数则跳过值交换部分，例如这个8点的倒位序，从000开始计算下面的倒位序二进制数，若倒序数最高位为0，则要将此位取反为1，其余位不变，然后顺序二进制数加一指向下一个。如果顺序数加一后小于倒序数，则两者可以交换。若倒序数最高位为1，则要将此位取反为0，再判断次高位若为0，则再将次高位取反为1。总结来说就是倒序数最高位若为0，则从0变为1，停止下一位的判断；最高位若是1，则从1变为0，还需要判断次高位，若次高位是从0变为1，则可以停止下一位的判断，否则，次高位也是从1变为0，那么就还需要判断下一位了。最高位取反为1，实际操作是此倒序数加N/2（序列点数的1/2点），若取反为0则实际操作是此倒序数减N/2；若次高位取反为1（0），则加（减）N/4，若次次高位取反为1（0），则加（减）N/8，并依此类推。设I为顺序二级制数，J为倒序二进制数，I=J=000。设常量N2为点数的一半，N/2，作为序列中点，设常量N1为序列总点数N-1，作为序列长度。首先判断I=J，I=J=000满足则需要进行倒序，跳过变址交换部分。设变量K=N2=100，判断KJ，是则J=J+K=000+100=100，I=I+1=001。再判断I=J，否则进入变址交换部分：将此时的I和J的储存单元内的值交换。再判断此时的KJ，否则J=J-K=100-100=000，K=K/2=010，再判断此时的KJ，是则J=J+K=000+010=010，I=I+1=001+1=010，判断I没到序列总点数N1=N-1，则返回到I=J,以此类推，一直循环判断和计算，直到I达到总点数后退出程序。下图是程序框图：假设已有原序列一维数组A(N)，N为原序列的长度点数。N2=N/2;N1=N-1;J=I=0;I=J?YNT=A(J);A(J)=A(I);A(I)=T;K=N2;KJ?J=J-K;NK=K/2;YJ=J+K;I=I+1;NI=N1?YFFT算法流程图：数据倒位序原序列先要倒位序，存入一维数组A(N)