快速傅里叶变换(FFT)

第4章快速傅里叶变换(FFT)4.1引言4.2基2FFT算法4.3进一步减少运算量的措施4.4分裂基FFT算法4.5离散哈特莱变换(DHT)4.1引言DFT是信号分析与处理中的一种重要变换。因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换(简称FFT)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。直到1965年Cooley和Tukey发现了DFT的一种快速算法以后，情况才发生了根本的变化。4.2基2FFT算法4.2.1直接计算DFT的特点及减少运算量的基本途径长度为N的有限长序列x(n)的DFT为考虑x(n)为复数序列的一般情况，对某一个k值，直接按(4.2.1)式计算X(k)值需要N次复数乘法、(N-1)次复数加法。因此，N点DFT的复乘次数等于N2,加法次数N(N-1).当N1时，，即N点DFT的乘法和加法运算次数均与N2成正比，当N较大时，运算量相等可观。10()(),0,1,,1NknNnXkxnWkN(4.2.1)2(1)NNN注意：通常将算术乘法和算术加法的次数作为计算复杂性的度量，因为这种方法使用起来很简单。如果在计算机上用软件实现这些算法，则乘法和加法的次数就直接与计算速度有关。但是，在常用的VLSI实现时，芯片的面积和功率要求往往是最重要的考虑因素，而它们有可能与算法的运算次数没有直接的关系。显然，把N点DFT分解为几个较短的DFT，可使乘法次数大大减少。另外，旋转因子WmN具有明显的周期性、对称性和可约性。其周期性表现为22()jmlNjmmlNmNNNNWeeW(4.2.2)其对称性表现为22[](1)mNmNmmNNNNNNmkmkNNN或者可约性表现在：nkmnkNmNnkmnkNNm()()2()21nNkNnknkNNNNNNkkNN4.2.2时域抽取法基2FFT基本原理FFT算法基本上分为两大类：时域抽取法FFT(DecimationInTimeFFT,简称DIT-FFT)和频域抽取法FFT(DecimationInFrequencyFFT,简称DIF-FFT)。下面介绍DIT-FFT算法。设序列x(n)的长度为N，且满足2,MNM为自然数按n的奇偶把x(n)分解为两个N/2点的子序列12()(2),0,1,12()(21),0,1,12NxrxrrNxrxrr则x(n)的DFT为/21/212(21)00/21/21221200()()()(2)(21)()()knknNNnnNNkrkrNNrrNNkrkkrNNNrrXkxnWxnWxrWxrWxrWWxrW偶数奇数由于22222/2jkrNjkrkrkrNNNWeeW所以/21/211/22/21200()()()()()NNkrkkrkNNNNrrXkxrWWxrWXkWXk其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，即/2111/210/2122/220()()[()]()()[()]NkrNrNkrNrXkxrWDFTxrXkxrWDFTxr(4.2.5)(4.2.6)由于X1(k)和X2(k)均以N/2为周期，且，所以X(k)又可表示为2NkkNNWW1212()()()0,1,12()()()0,1,122kNkNNXkXkWXkkNNXkXkWXkk(4.2.7)(4.2.8)2NkkNNWW图4.2.1蝶形运算符号CABA＋BCA－BCX1(k)X2(k)WNKX1(k)+WNKX2(k)X1(k)-WNKX2(k)经过一次分解后，计算复数乘和复数加的次数：复数乘：复数加：一次分解后，运算量减少近一半，故可以对N/2点DFT再作进一步分解。22(1)222NNNNNNNNN2（）222[(-1)]+22222图4.2.2N点DFT的一次时域抽取分解图(N=8)N/2点DFTWN0N/2点DFTWN1WN2WN3x(0)X1(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(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)与第一次分解相同，将x1(r)按奇偶分解成两个N/4长的子序列x3(l)和x4(l)，即3141()(2),0,1,,1()(21)4xlxlNlxlxl那么，X1(k)又可表示为/41/412(21)11/21/200/41/413/4/24/4003/24()(2)(21)()()()(),0,1,/21NNklklNNllNNklkklNNNllkNXkxlWxlWxlWWxlWXkWXkkN(4.2.9)式中/4133/430/4144/440()()[()]()()[()]NklNlNklNlXkxlWDFTxlXkxlWDFTxl同理，由X3(k)和X4(k)的周期性和的对称性，最后得到：13/2413/24()()(),0,1,,/41(/4)()()kNkNXkXkWXkkNXkNXkWXk(4.2.10)2mNW422kNkNNWW用同样的方法可计算出25/2625/26()()(),0,1,/41(/4)()()kNkNXkXkWXkkNXkNXkWXk(4.2.11)其中/4155/450/4166/4605262()()[()]()()[()]()(2),0,1,/41()(21)NklNlNklNlXkxlWDFTxlXkxlWDFTxlxlxllNxlxl图4.2.3N点DFT的第二次时域抽取分解图(N=8)N/4点DFTWN12WN12WN0WN1WN2WN3X1(0)X1(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)x(0)X3(0)X3(1)X4(0)X4(1)x(4)x(2)x(6)x(1)x(5)x(3)x(7)N/4点DFTN/4点DFTN/4点DFTWN02WN02图4.2.4N点DIT―FFT运算流图(N=8)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)4.2.3DIT-FFT算法与直接计算DFT运算量的比较运算流图有M级蝶型，每一级都有N/2个蝶型运算。每一级运算都需要N/2次复数乘和N次复数加(每个蝶形需要两次复数加法)。所以，M级运算总共需要的复数乘次数为2(2)log22MNNCMN复数加次数为例如，N=210=1024时221048576204.8(/2)log5120NNN2(2)logACNMNN图4.2.5FFT算法与直接计算DFT所需乘法次数的比较曲线MATLAB提供了一个fft的函数用于计算一个向量x的DFT。调用X=fft(x,N)就计算出N点的DFT。如果向量x的长度小于N，那么就将x补0。如果略去N，则DFT的长度就是x的长度。如果x是一个矩阵，那么fft(x,N)计算x中每一列的N点的DFT。fft由机器语言写成的，执行速度快。当N为2的幂次方，则使用基2FFT算法，如果不是，那么将N分解为若干素因子并用一个较慢的混合基FFT算法。如果N为某个素数，则fft算法就蜕化为原始的DFT算法。4.2.4DIT-FFT的运算规律及编程思想1.原位计算1)由图4.2.4可以看出，DIT―FFT的运算过程很有规律。N=2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成。2)同一级，每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据节点又同在一水平线上，即计算完一个蝶形后，所得的数据可立即存入原输入数据所占用的存储单元。3)经过M级运算后，原来存放输入序列数据的N个存储单元中依次存放X(k)的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位计算，可以大大节省内存。2.旋转因子的变化规律如上所述，N点DIT-FFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子WpN，称其为旋转因子，p称为旋转因子的指数.观察图4.2.4不难发现，第L级共有2L-1个不同的旋转因子。N=23=8时的各级旋转因子表示如下：L=1时，L=2时，L=3时，对N=2M的一般情况，第L级的旋转因子为12212,0,1,2,,212222,0,1,2,,212LMLLMpJLNLMLMLMPJJLNNNMLWWJN(4.2.12)(4.2.13)42222,0,0,1,0,1,2,3LLLpJJNNpJJNNpJJNN3.序列的倒序DIT-FFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种倒序是很有规律的。由于N=2M，所以顺序数可用M位二进制数(nM-1nM-2…n1n0)表示。图4.2.7形成倒序的树状图(N=23)01010101010101(n2n1n0)200004261537100010110001101011111表4.2.1顺序和倒序二进制数对照表4.2.5频域抽取法FFT(DIF-FFT)在基2快速算法中，频域抽取法FFT也是一种常用的快速算法，简称DIF-FFT。设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式：10/2110/2/21/21(/2)00/21/20()[()]()()()()()2[()()]2NknNnNNknknNNnnNNNknknNNNnnNkNknNNnXkDFTxnxnWxnWxnWNxnWxnWNxnWxnW/21,(1)1,kNkNkWk偶数奇数将X(k)分解成偶数组与奇数组，当k取偶数(k=2r,r=0,1,…,N/2-1)时/2120/212/20(2)[()()]2[()()]2NrnNnNrnNnNXrxnxnWNxnxnW(4.2.14)x1(n)当k取奇数(k=2r+1,r=0,1,…,N/2-1)时/21(21)0/21/20(21)[()()]2[()()]2NnrNnNnnrNNnNXrxnxnWNxnxnWW(4.2.15)将x1(n)和x2(n)分别代入(4.2.14)和(4.2.15)式，可得/211/20/212/20(2)()(21)()NrnNnNrnNnXrxnWXrxnW(4.2.16)x2(n)图4.2.10DIF-FFT蝶形运算流图符号4.2.6IDFT的高效算法上述FFT算法流图也可以用于离散傅里叶逆变换(InverseDiscreteFourierTransform，简称IDFT)。比较DFT和IDFT的运算公式：只要将DFT运算式中的系数改变为，最后乘以，就是IDFT的运算公式。故只要将上