第四章快速傅氏变换

第四章快速傅氏变换

引言：离散傅里叶变换在实际应用中是非常重要的，利用它可以计算信号的频谱、功率谱和线性卷积等。但是，如果使用定义式来直接计算DFT，当N很大时，即使使用高速计算机，所花的时间也太多。因此，如何提高计算DFT的速度，便成了重要的研究课题。1965年库利(Cooley)和图基(Tukey)在前人的研究成果的基础上提出了快速计算DFT的算法，之后，又出现了各种各样快速计算DFT的方法，这些方法统称为快速傅里叶变换(Fast

FourierTransform)，简称为FFT。快速傅里叶变换(FFT)是离散傅里叶变换(DFT)的快速算法它是DSP领域中的一项重大突破，它考虑了计算机和数字硬件实现的约束条件、研究了有利于机器操作的运算结构，使DFT的计算时间缩短了1~2个数量级，还有效地减少了计算所需的存储容量。FFT技术的应用极大地推动了DSP的理论和技术的发展，成为数字信号处理强有力的工具。其中在导出FFT算法之前，先估计一下直接计算DFT所需计算量。DFT的定义将方程组写成矩阵形式将DFT定义式展开成方程组同理用复数表示：用向量表示：从矩阵形式表示可以看出每计算一个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的计算予以改进。(1)对称性或(2)周期性(3)可约性或(4)其它首先考察WNnk的一些特殊性质：§4.1

基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-1因为WN2=WN/21，所以上式变为上式中的G(k)和H(k)都是N/2点的DFT。根据DFT的定义k=0,1,…,N/2-1(4.1)要用G(k)和H(k)来表达全部的X(k)值还必须考虑后半部分项数（k+N/2）。利用系数周期性，因为所以X(k)是N点序列，式4.1计算得到的只是前一半项数的结果k=0,1,…,N/2-1(4.2)将式4.1中的k用k+N/2代入，可得到2点的DFT流程图ab-1a-bWNka+bWNkWNk这样就把一个N点的DFT分解成了两个N/2点的DFT，只是最后求和的符号不同。式4.1与4.2的运算可用蝶形信号流图来表示图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)G(1)G(2)G(3)H(1)H(2)H(9)X(1)=G(1)+WN1H(1)X(2)=G(2)+WN2H(2)X(3)=G(3)+WN3H(3)WN1WN2WN3-1-1-1X(5)=G(1)-WN1H(1)X(6)=G(2)-WN2H(2)X(7)=G(3)-WN3H(3)G(0)H(0)WN0-1X(0)=G(0)+WN0H(0)X(4)=G(0)-WN0H(0)式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)k=0,1,…,N/4-1(4.4)k=0,1,…,N/4-1(4.5)k=0,1,…,N/4-1(4.6)M(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)的计算。x(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)将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。图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-1ab-1a-bWNka+bWNkWNk从图4中可看出几个特点：(1)流程图的每一级的基本计算单元都是一个蝶形；

(2)输入x(n)不按自然顺序排列，称为“混序”排列，而输出X(k)按自然顺序排列，称为“正序”排列，因而要对输入进行“变址”；(3)由于流程图的基本运算单元为蝶形，所以可以进行“同址”或“原位”计算。M=log2N由图可看出完成一个蝶形计算需一次复数乘法和两次复数加法二、运算量的减少（蝶形计算）任何一个N为2的整数幂（即N=2M）的DFT，都可以通过M次分解，最后成为2点的DFT来计算。M次分解构成了从x(n)到X(k)的M级迭代计算，每级由N/2个蝶形组成。ab-1a-bWNka+bWNkWNk下图表示了蝶形的一般形式表示。其输入和输出之间满足下列关系：大多数情况下复数乘法所花的时间最多，因此下面仅以复数乘法的计算次数为例来与直接计算进行比较。直接计算DFT需要的乘法次数为aD=N2，于是有即直接计算DFT所需复数乘法次数约为FFT的205倍。显然，N越大，FFT的速度优势越大。复数乘法次数：例：当N=1024时，因此，完成N点的时间抽选FFT计算的总运算量为复数加法次数：N所需复数乘法次数aD/aF直接计算DFT用FFT计算DFT24816326412825651210242048416642561024409616384655362621441048576419430414123280192448102423045120112644.04.05.48.012.821.436.664.0113.8204.8372.4不同N值所对应的两种计算方法的复数乘法次数和它们的比值ab-1a-bWNka+bWNkWNk三、同址(原位)计算图4包含M级迭代运算，每级由N/2个蝶形计算构成。蝶形计算的优点是可进行所谓同址(原位)计算。首先每一个蝶形运算都需要两个输入数据，计算结果也是两个数据，与其它结点的数据无关，其它蝶形运算也与这两结点的数据无关。因此，一个蝶形运算一旦计算完毕，原输入数据便失效了。这就意味着输出数据可以立即使用原输入数据结点所占用的内存。原来的数据也就消失了。这样输出、输入数据利用同一内存单元的这种蝶形计算称为同址计算。-1WN0x(0)x(4)x(0)+x(4)WN0x(0)-x(4)WN0M(0)M(1)M(0)M(1)类似地，M(2)和M(3)中的x(2)和x(6)进入运算器进行蝶形运算后的结果也被送回到M(2)和M(3)保存，等等。现在来考察第一级的计算规律设将输入x(0)，x(4)，x(2)，x(6)，x(1)，x(5)，x(3)，x(7)分别存入计算机的存储单元M(0),M(1),M(2),…，M(6)和M(7)中。首先，存储单元M(0)和M(1)中的数据x(0)和x(4)进入运算器并进行蝶形运算，运算后的结果便可送到M(0)和M(1)存储起来。第二级运算与第一级类似M(0)和M(2)存储单元中的数据进行蝶形运算后的结果被送回M(0)和M(2)保存，M(1)和M(3)中的数据进行蝶形运算后送回M(1)和M(3)保存，等等。这样一直到最后一级的最后一个蝶形运算完成。可以看出，每一级的蝶形的输入与输出在运算前后可以存储在同一地址(原来位置上)的存储单元中，这种同址运算的优点是可以节省存储单元，从而降低对计算机存储量的要求或降低硬件实现的成本。M(0)M(1)M(2)M(3)M(4)M(5)M(6)M(7)M(0)M(1)M(2)M(3)M(4)M(5)M(6)M(7)M(0)M(1)M(2)M(3)M(4)M(5)M(6)M(7)M(0)M(1)M(2)M(3)M(4)M(5)M(6)M(7)如图所示的N点FFT计算，只需N个存储存储单元X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)WN0WN1WN2WN3-1-1-1-1WN0-1WN0-1WN0-1WN0-1WN0-1WN0-1WN2WN2-1-1四、变址计算（倒序重排）从图4所示的流程图看出，同址计算要求输入x(n)是“混序”排列的。所谓输入为“混序”，并不是说输入是杂乱无章的，实际上它是有规律的。如果输入x(n)的序号用二进制码来表示，就可以发现输入的顺序恰好是正序输入的“码位倒置”。下表2列出了这种规律。自然序列二进制表示码位倒置混序x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x(000)x(001)x(010)x(011)x(100)x(101)x(110)x(111)表2x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)x(000)x(100)x(010)x(110)x(001)x(101)x(011)x(111)存储单元M(0)M(1)M(2)M(3)M(4)M(5)M(6)M(7)自然顺序x(n)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)码位倒置顺序x(l)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)变址图5码位倒置的变址处理在实际运算中，按码位倒置顺序输入数据x(n)，特别当N较大时，是很不方便的。因此，数据总是按自然顺序输入存储，然后通过“变址”运算将自然顺序转换成码位倒置顺序存储。实现这种转换的程序可用图5来说明。图中用n表示自然顺序的标号，用l表示码位倒置的标号当l＝n时，x(n)和x(l)不必互相调换。当l≠n时，必须将x(l)和x(n)互相调换，但只能调换一次，为此必须规定每当l>n时，要将x(l)和x(n)相互调换，即把原来存放x(n)的存储单元中的数据调入存储x(l)的存储单元中，而把原来存储x(l)的存储单元中的数据调入到存储x(n)的存储单元中。这样，按自然序输入的数据x(n)经过变址计算后变成了码位倒置的排列顺序，便可进入第一级的蝶形运算。

另外一些形式的时间抽选FFT算法流程图对于任何流程图，只要保持各节点所连支路及其传输系数不变，则不论节点位置怎样排列，所得到的流程图总是等效的，因而都能得到DFT的正确结果，只是数据的提取和存储次序不同而已。把图4中与x(4)水平相邻的所有节点和与x(1)水平相邻的所有节点交换，把与x(6)水平相邻的所有节点和与x(3)水平相邻的所有节点交换，而与x(0)、x(2)、x(5)和x(7)水平相邻各节点位置不变，就可以从图4得到图6。图6与图4的区别只是节点的排列不同，而支路传输比，即WN的各次幂保持不变。显然图6所示流程图的输入是正序（自然顺序）排列的，输出是码位倒置排列的，所以输出要进行变址计算。图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-1x(1)x(5)x(3)x(7)WN0WN2-1-1x(0)x(4)x(2)x(6)WN0-1WN0-1WN0-1WN0-1图6基2时间抽选法8点FFT流程图（输入正序）X(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)WN0WN0WN0WN0-1-1-1-1WN0WN0-1-1x(4)x(5)x(6)x(7)WN2WN2-1-1x(0)x(1)x(2)x(3)WN0-1WN1-1WN2-1WN3-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)WN0WN0WN0WN0-1-1-1-1WN0WN0x(4)x(5)x(6)x(7)WN2WN2x(0)x(1)x(2)x(3)WN0WN1WN2WN3-1-1-1-1-1-1-1-1另一种形式的流程图是将节点排列成输入和输出两者都是正序排列，但这类流程图不能进行同址计算，因而需要两列长度为N的复数存储器。其中r=0,1,…,N/2-1§4.2基2频率抽选法（桑德—图基算法）基2频率抽选法是在频域内将X(k)逐次分解成奇偶序列，再进行DFT计算设N＝2M，M为正整数。为推导方便，取N=23＝8。将频率偶奇分， 即X(k)={X(0),X(2),X(4),X(6),|X(1),X(3),X(5),X(7)}用X(2r)和X(2r+1)分别表示频率的偶数项和奇数项， 则有其中g(n)=x(n)+x(n+N/2)n=0,1,…,N/2-1同理其中h(n)=x(n)-x(n+N/2)n=0,1,…,N/2-1(4.7)(4.8)上面两式所表示的是N/2的DFT。在实际计算中，首先形成序列g(n)和h(n)，然后计算h(n)WNn，最后分别计算g(n)和h(n)WNn的N/2点DFT，便得到偶数输出点和奇数输出点的DFT。计算流程图如图7所示。N/2点DFTN/2点DFTx(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)WN0WN1WN2WN3-1-1-1-1图7N点的DFT分解成两个N/2点DFT的的频率抽选流程图N=8ab-1a-bWNka+bWNkWNkaa+bb-1(a–b)WNkWNk时间抽选蝶形单元频率抽选蝶形单元与时间抽选的FFT算法一样，图7所示的流程图的基本运算也是蝶形运算，但是频率抽选与时间抽选中的蝶形单元有所不同。N是2的整数幂，所以N/2仍然是偶数。这样可以将N/2点DFT的输出再分为偶数组和奇数组，即将N/2点的DFT计算进一步分解为两个N/4点的DFT计算，以此类推直到不能分解。如图8图8基2频率抽选法8点FFT流程图X(0)X(4)X(2)X(6)WN0WN1WN2WN3-1-1-1-1WN0WN2-1-1x(4)x(5)x(6)x(7)WN0WN2-1-1x(0)x(1)x(2)x(3)WN0-1WN0-1WN0-1WN0-1X(1)X(5)X(3)X(7)频率抽选与时间抽选的异同点：比较图4与图8可知，频率抽选FFT算法的计算量与时间抽选FFT算法的计算量相同。与时间抽选算法一样，频率抽选FFT算法也具有同址（原位）计算的优点。但是，与时间抽选不同的是，频率抽选FFT算法的信号输入为正序排列，输出为码位倒置排列，因此输出要进行变址计算。k=0,1,…,N/2-1在时间抽选FFT法中，第一次分解得到G(k)、H(k)分别为输入序列偶数点和奇数点的N/2点DFT(4.1)(4.2)§4.3

IFFT的计算方法（快速傅里叶反变换）一、IFFT算法k=0,1,…,N/2-1根据上式可画出蝶形图G(k)X(k)H(k)-1X(k+N/2)以此类推可求出x(n)的各点，反变换流程如图9根据式4.1与4.2可得到图98点IFFT流程图1/21/21/21/2X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)-1-1-1-1-1-1x(1)x(5)x(3)x(7)-1-1x(0)x(4)x(2)x(6)-1-1-1-11/21/21/21/21/21/21/21/2比较上面两式，可以看出，只要把DFT公式中的系数WNkn改为WN-kn

，并乘以系数1/N，将X(k)作为输入序列，将x(n)作为输出序列，就可用FFT算法来计算IDFT，这就得到了IFFT的算法。（例图10）在IFFT计算中经常把常量1/N分解成M个1/2连乘，即1/N＝(1/2)M，并且在M级的迭代运算中，每级的运算都分别乘上一个1/2因子。二、利用FFT求IFFTFFT算法同样可以应用于IDFT的计算。已知DFT和IDFT公式为x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)WN0WN-1WN-2WN-3-1-1-1-1WN0WN-2-1-1X(4)X(5)X(6)X(7)WN0WN-2-1-1WN0-1WN0-1WN0-1WN0-11/21/21/21/21/21/21/21/21/21/21/21/21/21/21/21/21/21/21/21/21/21/21/21/2图108点IFFT流程图（例：利用频率抽取FFT—图8）IFFT算法分类当把时间抽选FFT算法用于IFFT计算时，由于原来输入的时间序列x(n)现在变为频率序列X(k)，原来是将x(n)偶奇分的，而现在变成对X(k)进行偶奇分了，因此这种算法改称为频率抽选IFFT算法。类似地，当把频率抽选FFT算法用于计算IFFT时，应该称为时间抽选IFFT算法。（图10为时间抽选8点IFFT算法）与DFT比较，只需将X*(k)作为输入序列就可直接利用FFT流程具体步骤如下：1、求X(k)的共轭X*(k)，取共轭只需将虚部乘以(-1)2、以X*(k)作为输入序列，直接调用FFT程序，计算得到Nx*(n)3、对计算结果取共轭并乘以1/N即得到x(n)取共轭法：对DFT的反变换取共轭§4.4

基4FFT算法（略）基4FFT与基2FFT相比，复数乘法减少了约1/4，但复数加法却增加了。无论是软件还是硬件方法，复数乘法计算时间是复数加法的几十倍乃至上百倍，所以基4算法产生的主要目的是为了减少复数乘法。近年来，在一些高速数字信号处理器中，实现一次乘法运算与一次加法运算的时间完全一样，所以应用基4算法时就不一定合算了。§4.5

线性调频Z变换算法（CZT）在某些情况下，我们所需要的频率取样点并不均匀的分布在单位圆上；有时只在单位圆的某一部分；有时要求某一部分取样点密集（例如窄带信号）；有时甚至取样点不在单位圆上。（例如语音信号处理）沿Z平面上的一段螺线作等分角抽样Zk=AW-k

k=0,1,…,M-1且M≤N式中

如图：一、算法基本原理已知x(n)长度为N，则W0、A0为正实数Im(z)z平面Re(z)1A0φ0

θ0

Z0ZM-10Im(z)z平面Re(z)1A0φ0

θ0

Z0ZM-10由图可知：1、A0表示起始抽样点Z0的矢量半径长度，通常A0≤1，否则Z0将在单位圆之外；2、q0表示起始抽样点Z0相角，可正可负;3、f0表示两相邻抽样点间角度差，f0为正时，表示路径逆时针，f0为负时表示路径顺时针。4、W0表示螺线的伸展率，W0＞1表示随k增大螺线内缩，W0＜1表示随k增大螺线外伸。W0=1表示半径为A0的一段圆弧，又如果A0=1,则圆弧是单位圆一部分。与直接计算DFT相似，当N、M很大时，运算量很大；于是可采用布鲁斯坦（Bluestein）等式，将其转换为卷积和形式，从而可采用FFT算法。Bluestein等式令可得：0≤k≤M-1将Zk代入X(z)，有0≤k≤M-1计算过程如图：x(n)h(n)y(n)g(n)令则或 可想象为频率随时间n成线性增长的复指数序列，在雷达系统中，这种信号称为线性调频信号(Chirpsignal)，因此在这里称为线性调频Z变换算法。序列0≤k≤M-10≤k≤M-10≤k≤M-10≤k≤M-1由可知h(n)的长度为n=-(N-1)～M-1点数为L=N+M-1，y(n)的长度为N，