快速傅里叶变换

第一节引言一、快速付里叶变换FFT•有限长序列通过离散傅里叶变换(DFT)将其频域离散化成有限长序列.但其计算量太大(与N的平方成正比）,很难实时地处理问题,因此引出了快速傅里叶变换(FFT).•FFT并不是一种新的变换形式,它只是DFT的一种快速算法.并且根据对序列分解与选取方法的不同而产生了FFT的多种算法.•FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。二、FFT产生故事当时加文(Garwin)在自已的研究中极需要一个计算付里叶变换的快速方法。他注意到图基(J.W.Turkey)正在写有关付里叶变换的文章，因此详细询问了图基关于计算付里叶变换的技术知识。图基概括地对加文介绍了一种方法，它实质上就是后来的著名的库利(CooleyJ.W)图基算法。在加文的迫切要求下，库利很快设计出一个计算机程序。1965年库利--图基在计算数学、MathematicofComputation杂志上发表了著名的“机器计算付里级数的一种算法”文章，提出一种快速计算DFT的方法和计算机程序--揭开了FFT发展史上的第一页，促使FFT算法产生原因还有1967年至1968年间FFT的数字硬件制成，电子数字计算机的发明，使DFT的运算大大简化了。三、本章主要内容•1.直接计算DFT算法存在的问题及改进途径。•2.多种DFT算法(时间抽取算法DIT算法，频率抽取算法DIF算法，线性调频Z变换即CZT法)•3.FFT的应用第二节直接计算DFT算法存在的问题及改进途径一、直接计算DFT计算量•问题提出：设有限长序列x(n),非零值长度为N,计算对x(n)进行一次DFT运算，共需多大的运算工作量?1.比较DFT与IDFT之间的运算量10)()()(NnknNDFTWnxkXnx1,1,0Nk101()()()NIDFTknNkXkxnXkWN1,1,0Nn其中x(n)为复数，也为复数所以DFT与IDFT二者计算量相同。knNjknNeW22.以DFT为例，计算DFT复数运算量•计算一个X(k)(一个频率成分)值，运算量为•例k=1则要进行N次复数乘法+(N-1)次复数加法所以，要完成整个DFT运算[X(k)有N个点]，其计算量为：N*N次复数相乘和N*(N-1)次复数加法10)()1(NnnNWnxX3.一次复数乘法换算成实数运算量•复数运算要比加法运算复杂，需要的运算时间长。•一个复数乘法包括4次实数乘法和2次实数加法。(a+jb)(c+jd)=(ac-bd)+j(bc+ad)4次实数乘法2次实数加法4.计算DFT需要的实数运算量•每运算一个X(k)的值，需要进行4N次实数相乘和2N+2(N-1)=2(2N-1)次实数相加.整个DFT运算量为：4N2次实数相乘和2N(2N-1)次实数相加.由此看出：直接计算DFT时，乘法次数与加法次数都是和N2成比例的。当N很大时，所需工作量非常可观。])}Re[)](Im[]Im[)((Re[)]Im[)](Im[]Re[)]({(Re[)(10knNknNNnknNknNWnxWnxjWnxWnxkX例子•例1:当N=1024点时，直接计算DFT需要：N2=220=1048576次，即一百多万次的复乘运算这对实时性很强的信号处理(如雷达信号处理)来讲，它对计算速度有十分苛刻的要求--迫切需要改进DFT的计算方法，以减少总的运算次数。•例2：石油勘探，24道记录，每道波形记录长度5秒，若每秒抽样500点/秒，每道总抽样点数=500*5=2500点24道总抽样点数=24*2500=6万点DFT运算时间=N2=(60000)2=36*108次作业•P200页，第1题二、改善DFT运算效率的基本途径利用DFT运算的系数的固有对称性和周期性，改善DFT的运算效率。1.合并法：合并DFT运算中的某些项。2.分解法：将长序列DFT利用对称性和周期性，分解为短序列DFT。knNW利用DFT运算的系数的固有对称性和周期性，改善DFT的运算效率knNWknNW的对称性：nkNknNWW*)(knNW的周期性：)()(kNnNkNnNknN1)()(22kjkNNjkNNeeW因为：1,1,0Nk1)()(22njNnNjNnNeeW1,1,0Nn由此可得出：nkNknNNkNnN)()()(1sincos)(222/jeWNNjNNkNNkNWW)(2例子•例:1454)54(4941898178258利用以上特性,得到改进DFT直接算法的方法.(1)合并法：步骤1-分解成虚实部•合并DFT运算中的有些项•对虚实部而言•所以代入DFT中：*()()knkNnNNWWknNnNkNWWReRe)(knNnNkNWWImIm)((1)合并法：步骤2-代入DFT中1010]}Re)(ImIm)([Re{]}Im)(ImRe)({[Re)(NnnkNnkNNnnkNnkNWnxWnxjWnxWnxkX展开：1010]Im)][Re(Im)([Re)()(NnnkNnkNNnnkNWjWnxjnxWnxkX(1)合并法：步骤3-合并有些项nkNnNkNnkNWnNxnxWnNxWnxRe)](Re)([ReRe)(ReRe)(Re)(nkNnNkNnkNWnNxnxWnNxWnxIm)](Im)([ImIm)(ImIm)(Im)(knNnNkNWWReRe)(knNnNkNWWImIm)(根据：有：(1)合并法：步骤4-结论由此找出其它各项的类似归并方法：乘法次数可以减少一半。例：1545)15(515Re)]4(Re)1([ReRe)]4(Re)1([Re)15(ReRe)1(ReWxxWxxWxWx2、将长序列DFT利用对称性和周期性分解为短序列DFT--思路•因为DFT的运算量与N2成正比的•如果一个大点数N的DFT能分解为若干小点数DFT的组合，则显然可以达到减少运算工作量的效果。2、将长序列DFT利用对称性和周期性分解为短序列DFT--方法22)()()(NNNkXkXkX把N点数据分成二半：其运算量为：2)2(N2)2(N2N再分二半：44)()(NNkXkX44)()(NNkXkX+=22N2)4(N2)4(N2)4(N2)4(N+++=24N这样一直分下去，剩下两点的变换。2、将长序列DFT利用对称性和周期性分解为短序列DFT--结论•快速付里时变换(FFT)就是在此特性基础上发展起来的，并产生了多种FFT算法，其基本上可分成两大类：•按抽取方法分：时间抽取法(DIT);频率抽取法(DIF)•按“基数”分：基-2FFT算法；基-4FFT算法；混合基FFT算法；分裂基FFT算法•其它方法：线性调频Z变换(CZT法)第三节基--2按时间抽取的FFT算法Decimation-in-Time(DIT)(Coolkey-Tukey)一、算法原理•设输入序列长度为N=2M(M为正整数，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的FFT算法。也称为Coolkey-Tukey算法。•其中基数2----N=2M，M为整数.若不满足这个条件，可以人为地加上若干零值（加零补长）使其达到N=2M例子•设一序列x(n)的长度为L=9,应加零补长为•N=24=16应补7个零值。0123456789101112131415nx(n)二、算法步骤1.分组，变量置换DFT变换：10)()(NnknNWnxkX1,,0Nk先将x(n)按n的奇偶分为两组，作变量置换:当n=偶数时，令n=2r;当n=奇数时，令n=2r+1;得到：x(2r)=x1(r);x(2r+1)=x2(r);r=0…N/2-1;则可得其DFT为两部分：前半部分：后半部分：)()()(21kXWkXkXkN12/,,0Nk)()()2/(21kXWkXNkXkN2.代入DFT中)()()12()2()()21rxDFTrxDFTrxDFTrxDFTnxDFTkX（生成两个子序列x(0),x(2)…x(2r)偶数点x(1),x(3)…x(2r+1)奇数点12/,,0Nr代入DFT变换式：/21/212(21)00/21/212212002222/2/2()(2)(21)()()()()NNrkrkNNrrNNrkrkkNNNrrjjNNNNXkxrWxrWxrWxr3.求出子序列的DFT上式得：12/012/02/2/2212/012/02/2/11212/12/0212/02/1)12()()()2()()(12/1,0)()()()()NrNrrkNrkNNrNrrkNrkNkNkNrkNNrNrrkNWrxWrxkXWrxWrxkXNkkXWkXWWrxWrxkX其中：（12/,,0Nk4.结论1•一个N点的DFT被分解为两个N/2点DFT。X1(k),X2(k)这两个N/2点的DFT按照：12/1,0)()()21NkDFTNkXWkXkXkN中的前半部分点又合成（再应用W系数的周期性，求出用X1(k),X2(k)表达的后半部的X(k+N/2)的值。5.求出后半部的表示式12/012/022/2)2/(2/2212/012/012/1)2/(2/11)2/(2/2/)()()()2()()()()2(NrNrrkNkNrNNrNrrkNkNrNNkrNrkNkXWrxWrxkNXkXWrxWrxkNXWW看出：后半部的k值所对应的X1(k),X2(k)则完全重复了前半部分的k值所对应的X1(k),X2(k)的值。)()()(212/)2/kXWkXkX后半部分：又（6.结论2频域中的N个点频率成分为：)()()2/()()()(2121kXWkXNkXkXWkXkXkNkN后半部分：前半部分：结论：只要求出（0~N/2-1）区间内的各个整数k值所对应的X1(k),X2(k)值，即可以求出(0~N-1)整个区间内全部X(k)值，这就是FFT能大量节省计算的关键。7.结论3•由于N=2M，因此N/2仍为偶数，可以依照上面方法进一步把每个N/2点子序列，再按输入n的奇偶分解为两个N/4点的子序列，按这种方法不断划分下去，直到最后剩下的是2点DFT，两点DFT实际上只是加减运算。)1()0()1()1()0()0(2,)()(10xxXxxXNWnxkXNnknN时三、蝶形结即蝶式计算结构也即为蝶式信号流图上面频域中前/后半部分表示式可以用蝶形信号流图表示。作图要素：(1)左边两路为输入(2)右边两路为输出(3)中间以一个小圆表示加、减运算（右上路为相加输出、右下路为相减输出)(4)如果在某一支路上信号需要进行相乘运算，则在该支路上标以箭头，将相乘的系数标在箭头旁。(5)当支路上没有箭头及系数时，则该支路的传输比为1。X1(k)X2(k)kNW)()(21kXWkXkN)()(21kXWkXkN蝶形结描述的另一种方法X1(k)X2(k)kNW)()(21kXWkXkN)()(21kXWkXkN1例子：求N=23=8点FFT变换（1）先按N=8--N/2=4，做4点的DFT：)()()2/()()()(2121kXWkXNkXkXWkXkXkNkN12/,,0Nk先将N=8DFT分解成2个4点DFT：可知：时域上：x(0),x(2),x(4),x(6)