1、1,第5章 快速傅里叶变换,5.1 直接计算DFT算法存在的问题及改进途径 5.2 按时间抽取（DIT）的基-2 FFT算法 5.3 蝶形、同址、变址运算 5.4 按频率抽取（DIF）的基-2 FFT算法 5.5 快速傅里叶反变换（IFFT） 5.6 快速傅里叶变换FFT的应用 5.7 线性调频Z变换,2,5.1 直接计算DFT算法存在的问题及改进途径,5.1.1 直接计算DFT的运算量5.1.2 改进途径,3,5.1 直接计算DFT算法存在的问题及改进途径,快速傅里叶变换（Fast Fourier Transform，简称为FFT）并不是一种新的变换， 而是离散傅里叶变换（DFT）的一种快速
2、算法。由于有限长序列在其频域也可离散化为有限长序列（DFT），因此离散傅里叶变换（DFT）在数字信号处理中是非常有用的。但是，在相当长的时间里， 由于DFT的计算量太大，即使采用计算机也很难对问题进行实时处理，所以并没有得到真正的运用。 直到1965年库利（J.W.Cooley）和图基（J.W.Tukey）首次提出了DFT运算的一种快速算法，后来又有桑德（G.Sande）和图基（J.W.Tukey）的快速算法相继出现以后，情况才发生了根本的变化。人们开始认识到DFT运算的一些内在规律，从而很快地发展和完善了一套高速有效的运算方法，这就是现在人们普遍称之为快速傅里叶变换（FFT）的算法。 FFT
3、出现后使DFT的运算大大简化，运算时间一般可缩短一、二个数量级之多，从而使DFT的运算在实际中真正得到了广泛的应用。,4,5.1.1 直接计算DFT的运算量, 设x(n)为N点有限长序列，其DFT为,k=0, 1, , N-1,（5-1）,反变换（IDFT）为,n=0, 1, , N-1,（5-2）,二者的差别只在于WN的指数符号不同，以及差一个常数乘因子1/N，所以IDFT与DFT具有相同的运算工作量。 下面我们只讨论DFT的运算量。,5,一般来说，x(n)和WNnk都是复数，X(k)也是复数，因此每计算一个X(k)值，需要N次复数乘法和N-1次复数加法。而X(k)一共有N个点（k从0取到N
4、-1），所以完成整个DFT运算总共需要N2次复数乘法及N(N-1)次复数加法。 在这些运算中乘法运算要比加法运算复杂，需要的运算时间也多一些。因为复数运算实际上是由实数运算来完成的， 这时DFT运算式可写成,5.1.1 直接计算DFT的运算量,（5-3）,6,由此可见，一次复数乘法需用四次实数乘法和二次实数加法；一次复数加法需二次实数加法。 因而每运算一个X(k)需4N次实数乘法和2N+2(N-1)=2(2N-1)次实数加法。所以，整个DFT运算总共需要4N 2次实数乘法和2N(2N-1)次实数加法。上述统计与实际需要的运算次数稍有出入，因为某些WN nk可能是1或j，就不必相乘了，例如WN0
5、=1，WN= -1， WNN/4 =- j等就不需乘法。但是为了便于和其他运算方法作比较， 一般都不考虑这些特殊情况，而是把WN nk都看成复数，当N很大时，这种特例的影响很小。从上面的统计可以看到，直接计算DFT，乘法次数和加法次数都是和N 2成正比的，当N很大时，运算量是很可观的，有时甚至是无法忍受的。例如，当N=8时，DFT需64次复乘，计算量小，而当对一幅NN点的二维图像进行DFT变换，N=1024时直接计算DFT所需复乘次数为（N2)21012次，如用每秒可做10万次复数乘法的计算机，即使不考虑加法运算时间，也需要近3000小时。这对实时性很强的信号处理来说，要么提高计算机的计算速度
6、，而这样，对计算速度的要求太高了。另外，就只能通过改进对DFT的计算方法，以大大减少运算次数。,5.1.1 直接计算DFT的运算量,7,下面讨论减少运算量的途径。仔细观察DFT的运算就可看出， 利用系数（旋转因子） 的以下固有特性，就可减少运算量： （1） WNnk的对称性,（2） WNnk的周期性,5.1.2 改进途径,（3） WNnk的可约性,另外,8,这样，利用这些特性，使DFT运算中有些项可以合并，并能将长序列DFT分解为短序列的DFT进行运算。而前面已经说到，DFT的运算量是与N 2成正比的，所以N越小越有利，因而小点数序列的DFT比大点数序列的DFT的运算量要小。快速傅里叶变换算法
7、正是基于这样的基本思想而发展起来的。它的算法形式有很多种，但基本上可以分成两大类，即按时间抽取（Decimation-in-Time，缩写为DIT）法和按频率抽取（Decimation-in-Frequency，缩写为DIF）法。下面分别进行探讨。,5.1.2 改进途径,9,5.2 按时间抽取（DIT）的基-2 FFT算法,设序列点数为N2M，M为整数。若不满足这个条件，可以人为地加上若干个零值点，达到这一要求。这种N为2的整数幂的FFT称为基-2 FFT。按时间抽取（DIT）的基-2 FFT算法的基本出发点是，利用旋转因子WN nk的对称性和周期性，将一个长序列的DFT分解为一些逐次变小的D
8、FT来计算。分解过程遵循两条规则：对时间进行奇偶分解；对频率进行前后分解。设序列x(n)长度为N，且满足N=2M，M为正整数。按n的奇偶把x(n)分解为两个N/2点的子序列：,（5-4）,10,则可将DFT化为,由于 ， 故上式可表示成,（5-5）,5.2 按时间抽取（DIT）的基-2 FFT算法,11,式中，X1(k)与X2(k)分别是x1(r)及x2(r)的N/2点DFT：,（5-6）,（5-7）,由此，我们可以看到，一个N点DFT已分解成两个N/2点的DFT。 这两个N/2点的DFT再按照式（5-5）组合成一个N点DFT。这里应该看到X1(k)，X2(k)只有N/2个点，即k=0, 1,
9、 , N/2-1。而X(k)却有N个点，即k=0, 1, , N-1， 故用式（5-5）计算得到的只是X(k)的前一半的结果，要用X1(k)，X2(k)来表达全部的X(k)值，还必须应用系数的周期性， 即,5.2 按时间抽取（DIT）的基-2 FFT算法,12,这样可得到,（5-8）,同理可得,（5-9）,式（5-8）、式（5-9）说明了后半部分k值（N/2kN-1）所对应的X1(k)，2(k)分别等于前半部分k值（0kN/2-1）所对应的X1(k)，X2(k)。,5.2 按时间抽取（DIT）的基-2 FFT算法,13,再考虑到WkN 的以下性质:,这样，把式（5-8）、式（5-9）、式（5-
10、10）代入式（5-5），就可将X(k)表达为前后两部分：,（5-10）,（5-11）,（5-12）,5.2 按时间抽取（DIT）的基-2 FFT算法,14,图 5-1 时间抽取法蝶形运算流图符号,式（5-11）、式（5-12）的运算可以用图5-1的蝶形运算流图表示。流图表示法中，当支路没有系数时，该支路的传输系数为1。,5.2 按时间抽取（DIT）的基-2 FFT算法,15,图 5-2 按时间抽取将一个N点DFT分解为两个N/2点DFT（N=8）,5.2 按时间抽取（DIT）的基-2 FFT算法,采用此法，可将上面讨论的分解过程表示于图5-2中。此图为N238的情况，其中输出值X(0)到X(3
11、)由式（5-11）给出，而输出值X(4)到X(7)由式（5-12）给出。,16,可以看出，每个蝶形运算需要一次复数乘法X2(k) WN k及两次复数加（减）法。据此，一个N点DFT分解为两个N/2点DFT，每一个N/2点DFT只需(N/2)2=N2/4次复数乘法，N/2(N/2-1)次复数加法。两个N/2点DFT共需2(N/2)2=N2/2次复数乘法和N(N/2-1)次复数加法。此外，把两个N/2点DFT合成为N点DFT时，有N/2个蝶形运算，还需要N/2次复数乘法及2N/2=N次复数加法。因而通过第一步分解后，总共需要(N2/2)+(N/2)=N(N+1)/2N2/2次复数乘法和N(N/2-
12、1)+N=N2/2次复数加法。由此可见，通过这样分解后运算工作量差不多减少了一半。 ,5.2 按时间抽取（DIT）的基-2 FFT算法,17,（5-13）,5.2 按时间抽取（DIT）的基-2 FFT算法,既然如此，由于N=2M，因而N/2仍是偶数，可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。,18,且,式中：,（5-14）,（5-15）,图5-3给出N=8时，将一个N/2点DFT分解成两个N/4点DFT， 由这两个N/4点DFT组合成一个N/2点DFT的流图。,5.2 按时间抽取（DIT）的基-2 FFT算法,19,图 5-3 N/2点DFT分解为两个N/4点DF
13、T,5.2 按时间抽取（DIT）的基-2 FFT算法,20,X2(k)也可进行同样的分解：,式中：,5.2 按时间抽取（DIT）的基-2 FFT算法,将系数统一为WN/2 kWN 2k，则一个N8点DFT就可以分解为四个N/42点的DFT，如图5-4所示。而根据上面同样的分析可得，这样利用四个N/4点的DFT及两级蝶形组合运算来计算N点DFT，比只用一次分解蝶形组合方式的计算量又减少了大约一半。,21,图 5-4 一个N=8点DFT分解为四个N/4点DFT,5.2 按时间抽取（DIT）的基-2 FFT算法,22,最后，剩下的是2点DFT，对于此例N=8，就是四个N/4=2 点DFT，其输出为X
14、3(k)，X4(k)，X5(k)，X6(k)，k=0, 1，这由式（5-14）式（5-17）四个式子可以计算出来。例如，由式（5-14）可得：,式中， , 故上式不需要乘法。 类似地, 可求出X4(k)，X5(k)，X6(k)。,（5-16）,（5-17）,5.2 按时间抽取（DIT）的基-2 FFT算法,23,这些两点的DFT都可以用一个蝶形结表示如图5-5所示。由此可得一个按时间抽取运算的完整N=8点DFT流图，如图5-6所示。,5.2 按时间抽取（DIT）的基-2 FFT算法,图5-5 两点的DFT都可以用一个蝶形结,24,由上可见，这种方法的每一步分解，都是按输入序列在时间上的次序是属
15、于偶数还是属于奇数来分解为两个更短的子序列，所以称为“按时间抽取法”。,图 5-6 N=8 按时间抽取的FFT运算流图,5.2 按时间抽取（DIT）的基-2 FFT算法,25,5.3.1 蝶形计算 5.3.2 同址（原位）运算 5.3.3 变址（倒位序）运算 5.5.4 蝶形运算两节点的“距离”和Nr的确定 5.3.5 按时间抽取的FFT算法的其他形式流图,5.3 蝶形、同址、变址运算,26,5.3.1 蝶形计算,由按时间抽取法FFT的流图可见，当N=2M时，都可以进行M次分解，故而共有M级蝶形，每级都由N/2个蝶形运算组成，每个蝶形完成下述基本迭代运算：,其流图的一般形式如图5-7所示。,（
16、5-18a）（5-18b）,Xm(q),Xm(p),1,Xm+1(p),Xm+1(q),图5-7 按时间抽取蝶形的一般流图形式,27,完成每个蝶形需要一次复乘和二次复加，因而每级运算都需N/2次复乘和N次复加，这样M级运算总共需要,复乘数,复加数,5.3.1 蝶形计算,(5-19),(5-20),实际计算量与这个数字稍有不同，因为W0N=1，NN/2=-1，NN/4= -j，这几个系数都不用乘法运算，但是这些情况在直接计算DFT中也是存在的。此外，当N较大时，这些特例相对而言就很少。所以，为了统一作比较起见，下面都不考虑这些特例。 由于计算机上乘法运算所需的时间比加法运算所需的时间多得多，故以
17、乘法为例，直接DFT复数乘法次数是N2，FFT复数乘法次数是(N/2) log2N。 直接计算DFT与FFT算法的计算量之比为,28,由此可得，直接计算DFT与FFT算法所需运算量分别与点数N的关系曲线如图5-8所示，可以更直观地看出FFT算法的优越性，特别当点数N越大时，FFT的优点更为明显。,（5-20）,5.3.1 蝶形计算,例如，5.1节中举例处理一幅NN （N=1024）点的二维图像的DFT直接计算，如用每秒可做10万次复数乘法的计算机（不考虑加法运算时间），当时直接计算DFT所需时间为3000小时，用FFT算法复数乘法约为，仅为原来的10万分之一，现在只需要2 分钟。,图5-8 直
18、接计算DFT与FFT算法所需运算量分别与点数N的关系曲线,29,从图5-6可以看出这种运算是很有规律的，其每级（每列）计算都是由N/2 个蝶形运算构成的，每一个蝶形结构完成式（5-18a）、式（5-18b）的基本迭代运算。 蝶形运算如图5-7所示，由一次复乘和两次复加（减）组成。由图5-6的流图看出，某一列的任何两个节点p和q的节点变量进行蝶形运算后，得到结果为下一列p、q两节点的节点变量，而和其他节点变量无关，因而可以采用原位运算，即某一列的N个数据送到存储器后，经蝶形运算，其结果为下一列数据，它们以蝶形为单位仍存储在这同一组存储器中，直到最后输出，中间无需其他存储器。也就是蝶形的两个输出值
19、仍放回蝶形的两个输入所在的存储器中。每列的N/2 个蝶形运算全部完成后，再开始下一列的蝶形运算。这样存储器数据只需N个存储单元。下一级的运算仍采用这种原位方式，只不过进入蝶形结的组合关系有所不同。这种原位运算结构可以节省存储单元，降低设备成本。,5.3.2 同址（原位）运算,30,观察图5-6的同址计算结构，发现当运算完成后，FFT的输出X(k)按正常顺序排列在存储单元中，即按X(0)，X(1)，X(7)的顺序排列，但是这时输入x(n)却不是按自然顺序存储的，而是按x(0)，x(4)， ， x(7)的顺序存入存储单元，看起来好像是“混乱无序”的，实际上是有规律的，我们称之为倒位序。 造成倒位序
20、的原因是输入x(n)按标号n的偶奇的不断分组。 如果n用二进制数表示为(n2n1n0)2(当N=8=23时，二进制为三位)， 第一次分组，由图5-2看出，n为偶数（相当于n的二进制数的最低位n0=0）在上半部分，n为奇数（相当于n的二进制数的最低位 n0=1）在下半部分。 下一次则根据次最低位n1为“0”或是“1”来分偶奇（而不管原来的子序列是偶序列还是奇序列）， 如此继续分下去，直到最后N个长度为1的子序列。 图5-9的树状图描述了这种分成偶数子序列和奇数子序列的过程。,5.3.3 变址（倒位序）运算,31,图5-9 倒位序的形成,5.3.3 变址（倒位序）运算,32,蝶形、同址、变址运算一
21、般实际运算中， ，总是先按自然顺序将输入序列存入存储单元， 为了得到倒位序的排列，我们通过变址运算来完成。如果输入序列的自然顺序号I用二进制数（例如n2n1n0）表示，则其倒位序J对应的二进制数就是（n0n1n2），这样，在原来自然顺序时应该放x(I)的单元，现在倒位序后应放x(J)。 例如, N=8时，x(3)的标号是I=3，它的二进制数是011，倒位序的二进制数是110，即J=6，所以原来存放在x(011)单元的数据现在应该存放在x(110)内。表5-1列出了N=8时的自然顺序二进制数以及相应的倒位序二进制数。,5.3.3 变址（倒位序）运算,表5-1 N=8时的自然顺序二进制数和相应的倒
22、位序二进制数,33,由表5-1可见，自然顺序数I增加1，是在顺序数的二进制数最低位加1，向左进位。而倒序数J则是在二进制数最高位加1， 逢2向右进位。 例如，在（000）最高位加1，则得（100），再在（100）最高位加1，向右进位，则得（010）。因（100）最高位为1， 所以最高位加1要向次高位进位， 其实质是将最高位变为0，再在次高位加1。用这种算法，可以从当前任一倒序值求得下一个倒序值。对于N=2M，M为二进制数最高位，其权值为N/2，且从左向右二进制位的权值依次为N/4，N/8， 2，1。 因此，最高位加1相当于十进制运算J+N/2。如果最高位是0（JN/2）， 则直接由J+N/2得
23、下一个倒序值；如果最高位是1（JN/2），则要将最高位变为0（J J-N/2），次高位加1（J+N/4）。但次高位加1时，同样要判断次高位0、1值，如果为0（JN/4），则直接加1（J J+N/4）; 否则，将次高位变为0（J J-N/4)，再判断下一位；以此类推，直到完成最高位加1，逢2向右进位的运算。,5.3.3 变址（倒位序）运算,34,把按自然顺序存放在存储单元中的数据，换成FFT原位运算流图所要求的倒位序的变址功能如图5-10所示，当I=J时，不必调换，当IJ时， 必须将原来存放数据x(I)的存储单元内调入数据x(J)，而将存放x(J)的存储单元内调入x(I)。为了避免把已调换过的数
24、据再次调换，保证只调换一次（否则又回到原状），我们只需看J是否比I小。若J比I小，则意味着此x(I)在前边已和x(J)互相调换过，不必再调换了; 只有当JI时，才将原存放x(I)及存放x(J)的存储单元内的内容互换。这样就得到输入所需的倒位序列的顺序。可以看出，其结果与图5-5的要求是一致的。,5.3.3 变址（倒位序）运算,图5-10 N=8 倒位序的变址处理,35,以图 5-6 的8点FFT为例，其输入是倒位序的，输出是自然顺序的。 其第一级（第一列）每个蝶形的两节点间“距离”为1， 第二级每个蝶形的两节点“距离”为2， 第三级每个蝶形的两节点“距离”为4。 由此类推得，对N=2M点FFT
25、，当输入为倒位序， 输出为正常顺序时，其第m级运算，每个蝶形的两节点“距离”为2m-1。由于对第m级运算，一个DFT蝶形运算的两节点“距离”为2m-1， 因而式（5-18）可写成：,5.3.4 蝶形运算两节点的“距离”和Nr的确定,为了完成上两式运算，还必须知道系数Nr的变换规律。仔细观察图5-6的流图可以发现r的变换规律是： 把式（5-22）中蝶形运算两节点中的第一个节点标号值，即k值，表示成M位（注意N=2M）二进制数； 把此二进制数乘上2M-m，即将此M位二进制数左移M-m位（注意m是第m级运算），把右边空出的位置补零，此数即为所求r的二进制数。从图5-6看出，Nr因子最后一列有N/2种
26、，顺序为WN0， WN1， 其余可类推。,36,显然，对于任何流图，只要保持各节点所连的支路及传输系数不变，则不论节点位置怎么排列所得流图总是等效的，所得最后结果都是x(n)的DFT的正确结果，只是数据的提取和存放的次序不同而已。这样就可得到按时间抽取的FFT算法的若干其他形式流图。 将图5-6中和x(4)水平相连的所有节点和x(1)水平相连的所有节点位置对调，再将和x(6)水平相连的所有节点与和x(3)水平相连的所有节点对调，其余诸节点保持不变，可得图5-11的流图。图5-11与图5-6的蝶形相同，运算量也一样，不同点是： 数据存放的方式不同，图5-6是输入倒位序、输出自然顺序，图5-11是
27、输入自然顺序、输出倒位序；取用系数的顺序不同,图5-6的最后一列是按N0,N1,N2,N3 的顺序取用系数，且其前一列所用系数是后一列所用系数中具有偶数幂的那些系数（例如N0,N2,.）；图5-11是最后一列是按N0,N1,N2,N3 的顺序取用系数，且其前一列所用系数是后一列所用系数的前一半，这种流图最初由库利和图基给出的时间抽取法。经过简单变换，也可得输入与输出都是按自然顺序排列的流图以及其他各种形式的流图，可以参见其他参考文献。,5.3.5 按时间抽取的FFT算法的其他形式流图,37,图 5-11 时间抽取、 输入自然顺序、 输出倒位序的FFT流图,5.3.5 按时间抽取的FFT算法的其
28、他形式流图,38,5.4 按频率抽取（DIF）的基 -2 FFT算法,5.4.1 算法原理5.4.2 按频率抽取法的运算特点,39,5.4.1 算法原理,这种按频率抽取（DIF）的基-2 FFT算法推导过程遵循两个规则：对时间进行前后分解；对频率进行偶奇分解。设序列点数为N=2M，M为正整数。同按时间抽取（DIT）的基-2 FFT算法相同，若不满足这个条件，可以人为地加上若干个零值点，达到这一要求。按规则把输入序列按前一半、后一半分开（不是按偶数、 奇数分开），把N点DFT写成两部分。,k=0, 1, , N-1,40,式中，用的是 ，而不是 ，因而这并不是N/2点DFT。 由于 ，故 ，可得
29、,k=0, 1, , N-1,（5-23）,当k为偶数时，(-1)k=1；k为奇数时，(-1)k=-1。因此，按k的奇偶可将X(k)分为两部分:,（5-24）,5.4.1 算法原理,41,（5-25）,式（5-24）为前一半输入与后一半输入之和的N/2点DFT， 式（5-25）为前一半输入与后一半输入之差再与WNn之积的N/2点DFT。 令:,（5-27）,5.4.1 算法原理,42,图 5-12 频率抽取法蝶形运算单元,由此构成频率抽取法蝶形运算单元，如图5-12所示。,5.4.1 算法原理,这样，我们就把一个N点DFT按k的奇偶分解为两个N/2点的DFT了。 N=8时，上述分解过程示于图5
30、-13。与时间抽取法的推导过程一样，由于N=2M，N/2仍是一个偶数，因而可以将每个N/2点DFT的输出再分解为偶数组与奇数组，这就将N/2点DFT进一步分解为两个N/4 点DFT。 这两个N/4点DFT的输入也是先将N/2点DFT的输入上下对半分开后通过蝶形运算而形成的，图5-14示出了这一步分解的过程。这样的分解可以一直进行到第M次（N=2M），第M次实际上是做两点DFT，它只有加减运算。 然而，为了有统一的运算结构，仍然用一个系数为WN0的蝶形运算来表示， 这N/2个两点DFT的N个输出就是x(n)的N点DFT的结果X(k)。 图5-15表示一个N=8 完整的按频率抽取的基-2 FFT运
31、算结构。,43,图 5-13 按频率抽取的第一次分解,5.4.1 算法原理,44,图 5-14 按频率抽取的第二次分解,5.4.1 算法原理,45,图 5-15 按频率抽取的FFT（N=8）信号流图,5.4.1 算法原理,46,按频率抽取法的运算特点与时间抽取法基本相同。从图5-15可以看出，它也是通过(N/2)M个蝶形运算完成的。每一个蝶形结构完成下述基本迭代运算：,式中，m表示m列迭代，k,j为整数所在行数，此式的蝶形运算如图5-16所示，也需要一次复乘和两次复加。,5.4.2 按频率抽取法的运算特点,图 5-16 频率抽取法蝶形运算一般流图形式,47,由图 5-15 与图 5-6 相比较
32、，DIF法与DIT法的区别是： 图5-15的DIF输入是自然顺序， 输出是倒位序的，这与图5-6的DIT法正好相反。但这不是实质性的区别，因为DIF法与DIT法一样，都可将输入或输出进行重排，使二者的输入或输出顺序变成自然顺序或倒位序顺序。DIF的基本蝶形（图5-16）与DIT的基本蝶形（图5-7）则有所不同，这才是实质的不同，DIF的复数乘法只出现在减法之后，DIT则是先作复乘后再作加减法。就运算量来说DIF与DIT则是相同的，即都有M级（列）运算，每级运算需N/2个蝶形运算来完成，总共需要mF=(N/2) logN次复乘与aF=N logN次复加，DIF法与DIT法都可进行原位运算。频率抽
33、取FFT算法的输入是自然顺序，输出是倒位序的。因此运算完毕后，要通过变址计算将倒位序转换成自然位序，然后再输出。转换方法与时间抽取法相同。,5.4.2 按频率抽取法的运算特点,48,由按时间抽取法与按频率抽取法的基本蝶形（图5-6与图5-16）运算看出，如果将DIT的基本蝶形加以转置，就得到DIF的基本蝶形; 反过来, 将DIF的基本蝶形加以转置， 就得到DIT的基本蝶形， 因而DIT法与DIF法的基本蝶形是互为转置的。 按照转置定理， 两个流图的输入-输出特性必然相同。转置就是将流图的所有支路方向都反向，并且交换输入与输出，但节点变量值不交换，这样即可从图5-6得到图5-16或者从图5-16
34、得到图5-6，因而对每一种按时间抽取的FFT流图都存在一个按频率抽取的FFT流图。这样把图5-5、图5-11的流图分别加以转置，就可得到不同DIF的FFT流图。因此可以说，有多少种按时间抽取的FFT流图就存在多少种按频率抽取的FFT流图。频率抽取法与时间抽取法是两种等价的FFT运算。,5.4.2 按频率抽取法的运算特点,49,以上所讨论的FFT的运算方法同样可用于IDFT的运算，称为快速傅里叶反变换，简称为IFFT。由DFT和IDFT的定义（见式（5-28），将两定义式进行比较，可以看出，只要把DFT运算中的每一个系数WN nk改为WN -nk，并乘以系数1/N，就可以用FFT算法来计算IDF
35、T，也就得到了IFFT的算法。当把时间抽取FFT算法用于IFFT计算时，由于 原来输入的时间序列x(n)现在变成了频率序列X(k)，原来是将x(n)奇偶分组的，而现在变成对X(k)进行奇偶分组了，从而这种算法改称为频率抽取IFFT算法。类似地，当把频率抽取FFT算法用于计算IFFT时，应该称为时间抽取IFFT算法。在IFFT计算中经常把常量1/N分解为M各1/2连乘，即1/N=(1/2)M，并且在M级的迭代运算中，每级的运算都分别乘以1/2因子。具体的流图可以类似FFT的流图形式画出。,5.5 快速傅里叶反变换,（5-28）,50,上面这种IFFT算法虽然编程很方便，但要改变FFT的程序和参数
36、才能实现。现介绍另一种计算IFFT的方法，可以完全不必改动FFT程序。推导如下：可见：只要将X(k)取共轭，就可以直接利用FFT子程序，最后再将运算结果取一次共轭，并除以N就得到x(n)的值。即将X(k)的虚部乘以-1，即先取X(k)的共轭，得X*(k)；将X*(k)直接送入FFT程序即可得出Nx*(n)；最后再对运算结果取一次共轭变换，并乘以常数1/N，即可以求出IFFT变换的x(n)的值。因此，FFT运算和IFFT运算就可以共用一个子程序块，这样就很方便进行FFT和IFFT的计算。,5.5 快速傅里叶反变换,51,5.6.1 利用FFT对信号进行谱分析 5.6.2 线性卷积的FFT算法,5
37、.6 快速傅里叶变换FFT的应用,52,所谓谱分析就是计算信号的频谱，包括幅度谱、相位谱和功率谱。,1. 谱分析步骤 对连续信号进行谱分析的框图如图5-17所示。在图5-17中，前置低通滤波器 LPF（预滤波器）的引入，是为了消除或减少时域连续信号转换成序列时可能 出现的频谱混叠的影响。在实际工作中，离散时间信号x(n)的时宽是很长的， 甚至是无限长的（例如语音或音乐信号）。由于DFT的需要(实际应用FFT计 算)，必须把x(n)限制在一定的时间区间之内，即进行数据截断。数据的截断 相当于加窗处理（窗函数见7.3节）。 因此， 在计算FFT之前，用一个时域 有限的窗函数w(n)加到x(n)上是
38、非常必要的。,图 5-17 时域连续信号离散傅里叶分析的处理步骤,5.6.1 利用FFT对信号进行谱分析,53,xc(t)通过A/D变换器转换(忽略其幅度量化误差)成取样序列x(n)，其频谱用X(ej)表示，它是频率的周期函数，即,(5-29),式中，Xc(j)或 为xc(t)的频谱。,5.6.1 利用FFT对信号进行谱分析,在实际应用中，前置低通滤波器的阻带不可能是无限衰减的， 故由Xc(j) 周期延拓得到的X(ej)有非零重叠，即出现频谱混叠现象。由于进行FFT的需要，必须对序列x(n)进行加窗处理，即v(n)=x(n)w(n)，加窗对频域的影响，用卷积表示如下:,54,(5-30),最后
39、是进行FFT运算。 加窗后的DFT是,0kN-1,(5-31),式中，假设窗函数长度L小于或等于DFT长度N，为进行FFT运算，这里选择N为2的整数幂次即N=2m。,5.6.1 利用FFT对信号进行谱分析,55,有限长序列v(n)=x(n)w(n)的DFT相当于v(n)傅里叶变换的等间隔取样。,(5-32),V(k)便是sc(t)的离散频率函数。因为DFT对应的数字域频率间隔为=2/N，且模拟频率和数字频率间的关系为=， 其中=2f。所以，离散的频率函数第k点对应的模拟频率为:,(5-33),(5-34),5.6.1 利用FFT对信号进行谱分析,56,由式（5-34）很明显可看出，数字域频率间
40、隔=2/N对应的模拟域谱线间距为,(5-35),谱线间距，又称频谱分辨率（单位：Hz）。所谓频谱分辨率是指可分辨两频率的最小间距。它的意思是，如设某频谱分析的F=5Hz，那么信号中频率相差小于5Hz的两个频率分量在此频谱图中就分辨不出来。,5.6.1 利用FFT对信号进行谱分析,57,总结可见，谱分析中参数主要有： T为取样时间间隔（取样周期），（单位：s）； fs为取样频率（单位：Hz），fs=1/N； fh为信号频率的最高频率分量（单位：Hz）； tp为截取连续时间信号的样本长度（又称最小记录长度，单位：s）；(5-36) F为谱线间距，又称频率分辨率（单位：Hz）。(5-37)在实际应用
41、中, 要根据信号最高频率fh和频谱分辨率F的要求，来确定T、tp和N的大小。,5.6.1 利用FFT对信号进行谱分析,58,由以上讨论得到谱分析的步骤如下： （1）数据准备首先，由取样定理，为保证取样信号不失真，fs2fh(fh为信号频率的最高频率分量，也就是前置低通滤波器阻带的截止频率)， 即应使取样周期T满足然后，由频谱分辨率F和T确定N, ；为了使用FFT运算，这里选择N为2的幂次即N=2m，由式(5-37)可知，N大，分辨率好，但会增加样本记录时间tp。 最后，由N, T确定最小记录长度，tp=NT。 这样在记录长度中进行N点取样就得到要求的离散序列x(n)。,5.6.1 利用FFT对
42、信号进行谱分析,59,（2）用FFT计算信号的频谱应用前面讨论的FFT算法计算x(n)的频谱，得到X(k)。 （3）由X(k) 求幅度谱、相位谱和功率谱一个长度为N的时域离散序列x(n)，其离散傅里叶变换X(k)（离散频谱）是由实部和虚部组成的复数，即（5-38） 将X(k)写成极坐标形式 （5-39） 式中，|X(k)|称为幅度谱，argX(k)称为相频谱。 实际中常常用信号的功率谱表示，功率谱是幅频谱的平方，功率谱PSD定义为 （5-40） 功率谱具有突出主频率的特性，在分析带有噪声干扰的信号时特别有用。将 式（5-39）绘成的图形称为频谱图。由频谱图可以知道信号存在哪些频率分 量，它们就
43、是谱图中峰值对应的点。谱图中最低频率为k=0，对应实际频率为 0（即直流）；最高频率为k=N/2，对应实际频率为f= fs/2；对处于0，1，2， N/2上的任意点k，对应的实际频率为f=kF=k fs/N。,5.6.1 利用FFT对信号进行谱分析,60,例 5-1 有一频谱分析用的FFT处理器，其取样点数必须是2的整数幂， 假定没有采用任何特殊的数据处理措施，已给条件为： 频率分辨率10 Hz； 信号最高频率4kHz。试确定以下参量最小记录长度tp、最大取样间隔T（即最小取样频率）、在一个记录中的最少点数N。,5.6.1 利用FFT对信号进行谱分析,解： 由分辨率的要求确定最小长度tp,从信
44、号的最高频率确定最大可能的取样间隔T（即最小取样频率fs=1/T）。 按取样定理 即,从而最小记录点数N应满足,取,61,2. 谱分析的误差分析 利用FFT对连续信号进行傅里叶分析时可能造成一定的误差，从而产生如下现象。 （ 1） 频谱混叠失真在图5-17 画出的基本步骤中，A/D变换前利用前置低通滤波器进行预滤波，使xc(t)频谱中最高频率分量不超过fh。假设A/D变换器的取样频率为fs，按照奈奎斯特取样定理，为了不产生混叠， 必须满足,也就是取样间隔T满足,5.6.1 利用FFT对信号进行谱分析,（5-41）,一般应取,fs=(2.53.0)fh,如果不满足fs2fh，就会产生频谱混叠失真
45、。对于FFT来说，频率函数也要取 样，变成离散的序列，其取样间隔为F(即频率分辨率)。满足,62,从以上可以看出，信号的最高频率分量fh与频率分辨率F存在矛盾关系，要想fh增加，则时域取样间隔T就一定减小，从而fs就增加，若是固定N，必然要增加F， 即分辨率下降。反之，要提高分辨率（减小F），就要增加tp，当N给定时， 必然导致T的增加（fs减小）。 要不产生混叠失真，则必然会减小高频容量（信号的最高频率分量）fh。 要想兼顾高频容量fh与频率分辨率F，即一个性能提高而另一个性能不变（或也得以提高）的惟一办法就是增加记录长度的点数N， 即要满足 这个公式是未采用任何特殊数据处理（例如加窗处理）
46、的情况下，为实现基本FFT算法所必须满足的最低条件。如果加窗处理，相当于时域相乘，则频域周期卷积，必然加宽频谱分量， 频率分辨率就可能变差，为了保证频率分辨率不变，则须增加数据长度tp。,5.6.1 利用FFT对信号进行谱分析,63,（2） 栅栏效应利用FFT计算频谱，只给出离散点k=2k/N或k=2k/(NT)上的频谱取样值，也就是只给出了基频F的整数倍的频谱值，而不可能得到连续频谱函数， 这就像通过一个“栅栏”观看信号频谱，所以只能在离散点上看到信号频谱， 这种现象称之为“栅栏效应”, 如图5-18 所示。这时，如果在两个离散的谱线之间有一个特别大的频谱分量，就无法检测出来了。减小栅栏效应
47、的一个方法就是要使频域取样更密，即增加频域取样点数N，在不改变时域数据的情况下，必然是在数据末端添加一些零值点，使一个周期内的点数增加，但并不改变原有的记录数据。频谱取样为2k/N，N增加，必然使样点间距更近（单位圆上样点更多），谱线更密，谱线变密后原来看不到的谱分量就有可能看到了。 注意，补零以改变计算FFT的周期时，所用窗函数的宽度不能改变。换句话说，必须按照数据记录的原来的实际长度选择窗函数， 而不能按照补了零值点后的长度来选择窗函数。补零不能提高频率分辨率，这是因为数据的实际长度仍为补零前的数据长度。,5.6.1 利用FFT对信号进行谱分析,64,图 5-18 栅栏效应,5.6.1 利
48、用FFT对信号进行谱分析,65,（3） 频谱泄漏对信号进行FFT计算，首先必须使其变成有限时宽的信号，这就相当于信号在时域乘一个窗函数如矩形窗，窗内数据并不改变。时域相乘即v(n)=x(n)w(n)， 加窗对频域的影响，可用式（5-30）卷积公式表示,卷积的结果，造成所得到的频谱V(ej)与原来的频谱X(ej)不相同，有失真。这种失真最主要的是造成频谱的“扩散”（拖尾、 变宽），这就是所谓的“频谱泄漏”。由上可知，泄漏是由于我们截取有限长信号所造成的。对具有单一谱线的正弦波来说，它必须是无限长的。也就是说，如果我们输入信号是无限长的，那么FFT就能计算出完全正确的单一线频谱。可是我们不可能这么
49、做，而只能取有限长记录样本。如果在该有限长记录样本中，正弦信号又不是整数个周期时，就会产生泄漏。,5.6.1 利用FFT对信号进行谱分析,66,例如，一个周期为N=16的余弦信号x(n)=cos(6n16）截取一个周期长度的信号即x1(n)=cos(6n/16)R16(n), 其16点FFT的频谱图见 5-19（a)所示，若截取的长度为13，则其16点FFT的频谱图见 5-19(b)所示。 由此可见，频谱不再是单一的谱线，它的能量散布到整个频谱的各处。 这种能量散布到其他谱线位置的现象即为“频谱泄漏”。,5.6.1 利用FFT对信号进行谱分析,图5-19 16点FFT的频谱图,（a）,（b）,
50、67,应该说明，泄漏也会造成混叠，因为泄漏将会导致频谱的扩展， 从而使最高频率有可能超过折叠频率（fs/2），造成频率响应的混叠失真。 泄漏造成的后果是降低频谱的分辨率。此外，由于在主谱线两边形成很多旁瓣，引起不同频率分量间的干扰(简称谱间干扰)， 特别是强信号谱的旁瓣可能湮没弱信号的主谱线，或者把强信号谱的旁瓣误认为是另一信号的谱线，从而造成假信号，这样就会使谱分析产生较大偏差。在进行FFT运算时，时域截断是必然的，从而频谱泄漏和谱间干扰也是不可避免的。为尽量减小泄漏和谱间干扰的影响，需增加窗的时域宽度(频域主瓣变窄)，但这又导致运算量及存储量的增加；其次，数据不要突然截断，也就是不要加矩形窗，而是加各种缓变的窗（例如：三角形窗、升余弦窗、改进的升余弦窗等），使得窗谱的旁瓣能量更小，卷积后造成的泄漏减小，这个问题在FIR滤波器设计一章（第7章）中将会讨论到。,
