4-快速傅里叶变换资料

►Down2020/2/15◙Main2020/2/15►Down◄Up◙MainReturn§4.1引言§4.2基2FFT算法§4.3IFFT的计算方法§4.4实序列的FFT高效算法2020/2/15►Down◄Up◙MainReturn§4.1引言•离散傅里叶变换在实际应用中是非常重要的，利用它可以计算信号的频谱、功率谱和线性卷积等。•但是，如果使用定义式来直接计算DFT，当N很大时，即使使用高速计算机，所花的时间也太多。因此，如何提高计算DFT的速度，便成了重要的研究课题。•1965年库利(Cooley)和图基(Tukey)在前人的研究成果的基础上提出了快速计算DFT的算法，之后，又出现了各种各样快速计算DFT的方法，这些方法统称为快速傅里叶变换(FastFourierTransform)，简称为FFT。2020/2/15►Down◄Up◙MainReturn快速傅里叶变换(FFT)是离散傅里叶变换(DFT)的快速算法•它是DSP领域中的一项重大突破，它考虑了计算机和数字硬件实现的约束条件、研究了有利于机器操作的运算结构，使DFT的计算时间缩短了1~2个数量级，还有效地减少了计算所需的存储容量。•FFT技术的应用极大地推动了DSP的理论和技术的发展，成为数字信号处理强有力的工具。2020/2/15►Down◄Up◙MainReturn其中N1knNn0X(k)DFTx(n)x(n)W0kN-1N1knNk01x(n)IDFTX(k)X(k)W0nN-1NjN2NWe在导出FFT算法之前，先估计一下直接计算DFT所需计算量。DFT的定义2020/2/15►Down◄Up◙MainReturn将方程组写成矩阵形式00010(N1)NNN10111(N1)NNN(N1)0(N1)1(N1)(N1)NNNX(0)x(0)Wx(1)Wx(N1)WX(1)x(0)Wx(1)Wx(N1)WX(N1)x(0)Wx(1)WX(N1)W00010(N1)NNN10111(N1)NNN(N1)0(N1)1(N1)(N1)NNNX(0)x(0)(1)x(1)(N1)x(N1)将DFT定义式展开成方程组N1knNn0X(k)DFTx(n)x(n)W2020/2/15►Down◄Up◙MainReturnNXWx212N1NNN222(N1)2NNNNN12(N1)(N1)NNN11111同理1N1xWXN212(N1)NNN222(N1)21NNNNN12(N1)(N1)NNN11111用复数表示：NNNNNXWxRe[W]Re[x]Im[W]Im[x]jRe[W]Im[x]Re[x]Im[W]用向量表示：2020/2/15►Down◄Up◙MainReturn从矩阵形式表示可以看出每计算一个X(k)值需要N次复乘法和(N-1)次复数加法，因而计算N个X(k)值，共需N2次复乘法和N(N-1)次复加法。每次复乘法包括4次实数乘法和2次实数加法，每次复加法包括2次实数加法，因此计算N点的DFT共需要4N2次实数乘法和(2N2+2N·(N-1))次实数加法。当N很大时，这是一个非常大的计算量。例如当N=1024，N2=1048576。在实际中，N往往是比较大的，因此有必要对DFT的计算予以改进。2020/2/15►Down◄Up◙MainReturn(1)对称性jnkN2nkNWenknk(n)k(Nn)kNNNN(W)nknkn(k)n(Nk)NNNN(W)或(2)周期性n(kN)nknNnkNNNN(nN)knkkNnkNNNN(3)可约性nkmnkNmNWWnknk/mNN/mWW或N/2NW1(4)其它kN/2kNNWW2kkNN/2WW首先考察WNnk的一些特殊性质：2020/2/15►Down◄Up◙MainReturn§4.2基2FFT算法基2时间抽选法（库里—图基算法）一、基本原理：利用旋转因子WNnk的对称性和周期性，将一个长序列x(n)的DFT计算，在时域上按奇偶抽选分解成短序列的DFT来计算。设N＝2M，M为正整数。为推导方便，取N＝23＝8，即离散时间信号为x(n)={x(0),x(1),x(2),x(3),x(4),x(5),x(6),x(7)}将序列x(n)分为奇偶两组，一组序号为偶数，另一组序号为奇数，即{x(0),x(2),x(4),x(6)|x(1),x(3),x(5),x(7)}分别表示为：g(r)=x(2r)——偶数项h(r)=x(2r+1)——奇数项r=0,1,…,N/2-12020/2/15►Down◄Up◙MainReturn因为WN2=WN/21，所以上式变为上式中的G(k)和H(k)都是N/2点的DFT。N1nknknkNNNn0nnX(k)x(n)Wx(n)Wx(n)W偶数奇数N/21N/212rkk2rkNNNr0r0g(r)(W)Wh(r)(W)N/21N/21rkkrkN/2NN/2r0r0X(k)g(r)WWh(r)W根据DFT的定义N/21N/212rk(2r1)kNNr0r0x(2r)Wx(2r1)WkNG(k)WH(k)k=0,1,…,N/2-1(4.1)2020/2/15►Down◄Up◙MainReturn要用G(k)和H(k)来表达全部的X(k)值还必须考虑后半部分项数（k+N/2）。利用系数周期性，因为所以rkr(kN/2)N/2N/2WWN/2N/2W1kN/2kNNWWN/21N/21r(kN/2)(kN/2)r(kN/2)N/2NN/2n0n0NX(k)g(r)WWh(r)W2N/21N/21rkkrkN/2NN/2n0n0NX(k)g(r)WWh(r)W2X(k)是N点序列，式4.1计算得到的只是前一半项数的结果kNG(k)WH(k)k=0,1,…,N/2-1(4.2)将式4.1中的k用k+N/2代入，可得到2020/2/15►Down◄Up◙MainReturn2点的DFT流程图ab-1a-bWNka+bWNkWNkkNX(k)G(k)WH(k)kNNX(k)G(k)WH(k)2这样就把一个N点的DFT分解成了两个N/2点的DFT，只是最后求和的符号不同。式4.1与4.2的运算可用蝶形信号流图来表示2020/2/15►Down◄Up◙MainReturn图1N点的DFT分解成两个N/2点DFT的的时间抽选流程图N=8按照式4.1和式4.2可画出下图所示的信号流程图。N/2点DFTN/2点DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)WN1G(1)H(1)X(1)=G(1)+WN1H(1)-1X(5)=G(1)-WN1H(1)G(0)H(0)WN0-1X(0)=G(0)+WN0H(0)X(4)=G(0)-WN0H(0)X(6)=G(2)-WN2H(2)X(2)=G(2)+WN2H(2)G(2)H(2)-1WN2X(7)=G(3)-WN3H(3)X(3)=G(3)+WN3H(3)G(3)H(3)-1WN32020/2/15►Down◄Up◙MainReturnN/21N/41N/41rk2mk(2m1)kN/2N/2N/2r0m0m0G(k)g(n)Wg(2m)Wg(2m1)WN/41N/41mk2kmkN/4NN/4m0m0g(2m)WWg(2m1)W2kNM(k)WN(k)式4.1和式4.2把原N点DFT计算分解成两个N/2点DFT计算。•照此可进一步把每个N/2点DFT的计算再各分解成两个N/4点DFT的计算。•具体说来，是把{x(0)，x(2)，x(4)，x(6)}和{x(1)，x(3)，x(5)，x(7)}分为{x(0)，x(4)|x(2)，x(6)}和{x(1)，x(5)|x(3)，x(7)}。•G(k)和H(k)分别计算如下:k=0,1,…,N/4-1(4.3)2020/2/15►Down◄Up◙MainReturnN/41N/412m(kN/4)(2m1)(kN/4)N/2N/2m0m0NG(k)g(2m)Wg(2m1)W4N/41N/41mk2kmkN/4NN/4m0m0g(2m)WWg(2m1)W2kNM(k)WN(k)k=0,1,…,N/4-1(4.4)N/41N/41mk2kmkN/4NN/4m0m0H(k)h(2m)WWh(2m1)W2kNP(k)WQ(k)N/41N/41mk2kmkN/4NN/4m0m0NH(k)h(2m)WWh(2m1)W42kNP(k)WQ(k)k=0,1,…,N/4-1(4.5)k=0,1,…,N/4-1(4.6)2020/2/15►Down◄Up◙MainReturnM(0)M(1)N(0)N(1)G(0)G(1)G(2)G(3)WN0WN2-1-1N/4点DFTx(0)x(4)x(2)x(6)N/4点DFT图2N/2点的DFT分解成两个N/4点DFT的的时间抽选流程图N=8这样，用式(4.3)～(4.6)就可计算图1中的两组N/2点DFT。图2所示的是其中一组G(k)的计算。2020/2/15►Down◄Up◙MainReturnx(0)x(4)x(2)x(6)WN0WN2-1-1N/4点DFTN/4点DFTx(1)x(5)x(3)x(7)N/4点DFTN/4点DFTWN0WN2-1-1将图1与图2合并，便得到图3所示的信号流程图。X(0)X(1)X(2)X(3)WN0WN1WN2WN3-1-1-1-1X(4)X(5)X(6)X(7)G(0)G(1)G(2)G(3)H(0)H(1)H(2)H(3)2020/2/15►Down◄Up◙MainReturn将2点DFT的信号流程图移入图3，得到图4所示的8点时间抽选的完整的FFT流程图。2点的DFT流程图x(0)-1x(0)+x(4)WN0WN0x(4)x(0)-x(4)WN0N=8，图3中N/4点的DFT就是2点的DFT。2020/2/15►Down◄Up◙MainReturn图4基2时间抽选法8点FFT流程图X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)WN0WN1WN2WN3-1-1-1-1WN0WN2-1-1WN0WN2-1-1x(1)x(5)x(3)x(7)x(0)x(4)x(2)x(6)WN0-1WN0-1WN0-1WN0-12020/2/15►Down◄Up◙MainReturnab-1a-bWNka+bWNkWNk从图4中可看出几个特点：(1)流程图的每一级的基本计算单元都是一个蝶形；(2)输入x(n)不按自然顺序排列，称为“混序”排列，而输出X(k)按自然顺序排列，称为“正序”排列，因而要对输入进行“变址”；(3)由于流程图的基本运算单元为蝶形，所以可以进行“同址”或“原位”计算。2020/2/15►Down◄Up◙MainReturnM=log2N由图可看出完成一个蝶形计算需一次复数乘法和两次复数加法二、运算量的减少（蝶形计算）•任何一个N为2的整数幂（即N=2M）的DFT，都可以通过M次分解，最后成为2点的DFT来计算。M次分解构成了从x(n)到X(