chapter-5--武汉工程大学课件

X2024/5/4信息学科立体化教材第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变换变换1X2024/5/4信息学科立体化教材5.1 直接计算DFT算法存在的问题及改进途径5.1.1 直接计算直接计算DFT的运算量的运算量 5.1.2 改进途径改进途径2X2024/5/4信息学科立体化教材5.1 直接计算直接计算DFT算法存在的问题及改进途径算法存在的问题及改进途径快速傅里叶变换（FastFourierTransform，简称为FFT）并不是一种新的变换，而是离散傅里叶变换（DFT）的一种快速算法。由于有限长序列在其频域也可离散化为有限长序列（DFT），因此离散傅里叶变换（DFT）在数字信号处理中是非常有用的。但是，在相当长的时间里，由于DFT的计算量太大，即使采用计算机也很难对问题进行实时处理，所以并没有得到真正的运用。直到1965年库利（J.W.Cooley）和图基（J.W.Tukey）首次提出了DFT运算的一种快速算法，后来又有桑德（G.Sande）和图基（J.W.Tukey）的快速算法相继出现以后，情况才发生了根本的变化。人们开始认识到DFT运算的一些内在规律，从而很快地发展和完善了一套高速有效的运算方法，这就是现在人们普遍称之为快速傅里叶变换（FFT）的算法。FFT出现后使DFT的运算大大简化，运算时间一般可缩短一、二个数量级之多，从而使DFT的运算在实际中真正得到了广泛的应用。3X2024/5/4信息学科立体化教材5.1.1 直接计算直接计算DFT的运算量的运算量设x(n)为N点有限长序列，其DFT为k=0,1,N1（51）反变换（IDFT）为n=0,1,N1（52）二者的差别只在于WN的指数符号不同，以及差一个常数乘因子1/N，所以IDFT与DFT具有相同的运算工作量。下面我们只讨论DFT的运算量。4X2024/5/4信息学科立体化教材一般来说，x(n)和WNnk都是复数，X(k)也是复数，因此每计算一个X(k)值，需要N次复数乘法和N1次复数加法。而X(k)一共有N个点（k从0取到N1），所以完成整个DFT运算总共需要N2次复数乘法及N(N1)次复数加法。在这些运算中乘法运算要比加法运算复杂，需要的运算时间也多一些。因为复数运算实际上是由实数运算来完成的，这时DFT运算式可写成5.1.1 直接计算直接计算DFT的运算量的运算量（53）5X2024/5/4信息学科立体化教材由此可见，一次复数乘法需用四次实数乘法和二次实数加法；一次复数加法需二次实数加法。因而每运算一个X(k)需4N次实数乘法和2N+2(N-1)=2(2N-1)次实数加法。所以，整个DFT运算总共需要4N 2次实数乘法和2N(2N-1)次实数加法。上述统计与实际需要的运算次数稍有出入，因为某些WN nk可能是1或j，就不必相乘了，例如WN0=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小时。这对实时性很强的信号处理来说，要么提高计算机的计算速度，而这样，对计算速度的要求太高了。另外，就只能通过改进对DFT的计算方法，以大大减少运算次数。5.1.1 直接计算直接计算DFT的运算量的运算量6X2024/5/4信息学科立体化教材下面讨论减少运算量的途径。仔细观察DFT的运算就可看出，利用系数（旋转因子）的以下固有特性，就可减少运算量：（1）WNnk的对称性（2）WNnk的周期性5.1.2 5.1.2 改进途径改进途径（3）WNnk的可约性另外7X2024/5/4信息学科立体化教材这样，利用这些特性，使DFT运算中有些项可以合并，并能将长序列DFT分解为短序列的DFT进行运算。而前面已经说到，DFT的运算量是与N 2成正比的，所以N越小越有利，因而小点数序列的DFT比大点数序列的DFT的运算量要小。快速傅里叶变换算法正是基于这样的基本思想而发展起来的。它的算法形式有很多种，但基本上可以分成两大类，即按时间抽取（DecimationinTime，缩写为DIT）法和按频率抽取（DecimationinFrequency，缩写为DIF）法。下面分别进行探讨。5.1.2 5.1.2 改进途径改进途径8X2024/5/4信息学科立体化教材5.2 按时间抽取（DIT）的基-2 FFT算法 设序列点数为N2M，M为整数。若不满足这个条件，可以人为地加上若干个零值点，达到这一要求。这种N为2的整数幂的FFT称为基2FFT。按时间抽取（DIT）的基2FFT算法的基本出发点是，利用旋转因子WN nk的对称性和周期性，将一个长序列的DFT分解为一些逐次变小的DFT来计算。分解过程遵循两条规则：对时间进行奇偶分解；对频率进行前后分解。设序列x(n)长度为N，且满足N=2M，M为正整数。按n的奇偶把x(n)分解为两个N/2点的子序列：（54）9X2024/5/4信息学科立体化教材则可将DFT化为由于，故上式可表示成（55）5.2 按时间抽取（DIT）的基-2 FFT算法 10X2024/5/4信息学科立体化教材式中，X1(k)与X2(k)分别是x1(r)及x2(r)的N/2点DFT：（56）（57）由此，我们可以看到，一个N点DFT已分解成两个N/2点的DFT。这两个N/2点的DFT再按照式（55）组合成一个N点DFT。这里应该看到X1(k)，X2(k)只有N/2个点，即k=0,1,N/21。而X(k)却有N个点，即k=0,1,N1，故用式（55）计算得到的只是X(k)的前一半的结果，要用X1(k)，X2(k)来表达全部的X(k)值，还必须应用系数的周期性，即5.2 按时间抽取（DIT）的基-2 FFT算法 11X2024/5/4信息学科立体化教材这样可得到（58）同理可得（59）式（58）、式（59）说明了后半部分k值（N/2kN1）所对应的X1(k)，2(k)分别等于前半部分k值（0kN/21）所对应的X1(k)，X2(k)。5.2 按时间抽取（DIT）的基-2 FFT算法 12X2024/5/4信息学科立体化教材再考虑到WkN 的以下性质:这样，把式（58）、式（59）、式（510）代入式（55），就可将X(k)表达为前后两部分：（510）（511）（512）5.2 按时间抽取（DIT）的基-2 FFT算法 13X2024/5/4信息学科立体化教材图51时间抽取法蝶形运算流图符号式（511）、式（512）的运算可以用图51的蝶形运算流图表示。流图表示法中，当支路没有系数时，该支路的传输系数为1。5.2 按时间抽取（DIT）的基-2 FFT算法 14X2024/5/4信息学科立体化教材图52按时间抽取将一个N点DFT分解为两个N/2点DFT（N=8）5.2 按时间抽取（DIT）的基-2 FFT算法 采用此法，可将上面讨论的分解过程表示于图52中。此图为N238的情况，其中输出值X(0)到X(3)由式（511）给出，而输出值X(4)到X(7)由式（512）给出。15X2024/5/4信息学科立体化教材可以看出，每个蝶形运算需要一次复数乘法X2(k)WN k及两次复数加（减）法。据此，一个N点DFT分解为两个N/2点DFT，每一个N/2点DFT只需(N/2)2=N2/4次 复 数 乘 法，N/2(N/21)次 复 数 加 法。两 个 N/2点 DFT共 需2(N/2)2=N2/2次复数乘法和N(N/21)次复数加法。此外，把两个N/2点DFT合成为N点DFT时，有N/2个蝶形运算，还需要N/2次复数乘法及2N/2=N次复数加法。因而通过第一步分解后，总共需要(N2/2)+(N/2)=N(N+1)/2N2/2次复数乘法和N(N/21)+N=N2/2次复数加法。由此可见，通过这样分解后运算工作量差不多减少了一半。5.2 按时间抽取（DIT）的基-2 FFT算法 16X2024/5/4信息学科立体化教材（513）5.2 按时间抽取（DIT）的基-2 FFT算法 既然如此，由于N=2M，因而N/2仍是偶数，可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。17X2024/5/4信息学科立体化教材且式中：（514）（515）图53给出N=8时，将一个N/2点DFT分解成两个N/4点DFT，由这两个N/4点DFT组合成一个N/2点DFT的流图。5.2 按时间抽取（DIT）的基-2 FFT算法 18X2024/5/4信息学科立体化教材图53N/2点DFT分解为两个N/4点DFT5.2 按时间抽取（DIT）的基-2 FFT算法 19X2024/5/4信息学科立体化教材X2(k)也可进行同样的分解：式中：5.2 按时间抽取（DIT）的基-2 FFT算法 将系数统一为WN/2 kWN 2k，则一个N8点DFT就可以分解为四个N/42点的DFT，如图54所示。而根据上面同样的分析可得，这样利用四个N/4点的DFT及两级蝶形组合运算来计算N点DFT，比只用一次分解蝶形组合方式的计算量又减少了大约一半。20X2024/5/4信息学科立体化教材图54一个N=8点DFT分解为四个N/4点DFT5.2 按时间抽取（DIT）的基-2 FFT算法 21X2024/5/4信息学科立体化教材最后，剩下的是2点DFT，对于此例N=8，就是四个N/4=2点DFT，其输出为X3(k)，X4(k)，X5(k)，X6(k)，k=0,1，这由式（514）式（517）四个式子可以计算出来。例如，由式（514）可得：式中，,故上式不需要乘法。类似地,可求出X4(k)，X5(k)，X6(k)。（516）（517）5.2 按时间抽取（DIT）的基-2 FFT算法 22X2024/5/4信息学科立体化教材这些两点的DFT都可以用一个蝶形结表示如图55所示。由此可得一个按时间抽取运算的完整N=8点DFT流图，如图56所示。5.2 按时间抽取（DIT）的基-2 FFT算法 图55两点的DFT都可以用一个蝶形结 x(0)x(1)X(0)X(1)123X2024/5/4信息学科立体化教材由上可见，这种方法的每一步分解，都是按输入序列在时间上的次序是属于偶数还是属于奇数来分解为两个更短的子序列，所以称为“按时间抽取法”。图56N=8按时间抽取的FFT运算流图5.2 按时间抽取（DIT）的基-2 FFT算法 24X2024/5/4信息学科立体化教材5.3.1蝶形计算5.3.2同址（原位）运算5.3.3变址（倒位序）运算5.5.4蝶形运算两节点的“距离”和Nr的确定5.3.5按时间抽取的FFT算法的其他形式流图5.3 蝶形、同址、变址运算25X2024/5/4信息学科立体化教材5.3.1 蝶形计算由按时间抽取法FFT的流图可见，当N=2M时，都可以进行M次分解，故而共有M级蝶形，每级都由N/2个蝶形运算组成，每个蝶形完成下述基本迭代运算：其流图的一般形式如图57所示。（518a）（518b）NrXm(q)Xm(p)1Xm+1(p)Xm+1(q)图57按时间抽取蝶形的一般流图形式26X2024/5/4信息学科立体化教材完成每个蝶形需要一次复乘和二次复加，因而每级运算都需N/2次复乘和N次复加，这样M级运算总共需要复乘数复加数5.3.1 蝶形计算(5-19)(520)实际计算量与这个数字稍有不同，因为W0N=1，NN/2=1，NN/4=j，这几个系数都不用乘法运算，但是这些情况在直接计算DFT中也是存在的。此外，当N较大时，这些特例相对而言就很少。所以，为了统一作比较起见，下面都不考虑这些特例。由于计算机上乘法运算所需的时间比加法运算所需的时间多得多，故以乘法为例，直接DFT复数乘法次数是N2，FFT复数乘法次数是(N/2)log2N。直接计算DFT与FFT算法的计算量之比为27X2024/5/4信息学科立体化教材由此可得，直接计算DFT与FFT算法所需运算量分别与点数N的关系曲线如图58所示，可以更直观地看出FFT算法的优越性，特别当点数N越大时，FFT的优点更为明显。（520）5.3.1 蝶形计算例如，5.1节中举例处理一幅NN（N=1024）点的二维图像的DFT直接计算，如用每秒可做10万次复数乘法的计算机（不考虑加法运算时间），当时直接计算DFT所需时间为3000小时，用FFT算法复数乘法约为，仅为原来的10万分之一，现在只需要2分钟。图图5-8 直接计算直接计算DFT与与FFT算法所需运算量分别与点数算法所需运算量分别与点数N的关系曲线的关系曲线28X2024/5/4信息学科立体化教材从图56可以看出这种运算是很有规律的，其每级（每列）计算都是由N/2个蝶形运算构成的，每一个蝶形结构完成式（518a）、式（518b）的基本迭代运算。蝶形运算如图57所示，由一次复乘和两次复加（减）组成。由图56的流图看出，某一列的任何两个节点p和q的节点变量进行蝶形运算后，得到结果为下一列p、q两节点的节点变量，而和其他节点变量无关，因而可以采用原位运算，即某一列的N个数据送到存储器后，经蝶形运算，其结果为下一列数据，它们以蝶形为单位仍存储在这同一组存储器中，直到最后输出，中间无需其他存储器。也就是蝶形的两个输出值仍放回蝶形的两个输入所在的存储器中。每列的N/2个蝶形运算全部完成后，再开始下一列的蝶形运算。这样存储器数据只需N个存储单元。下一级的运算仍采用这种原位方式，只不过进入蝶形结的组合关系有所不同。这种原位运算结构可以节省存储单元，降低设备成本。5.3.2 同址（原位）运算29X2024/5/4信息学科立体化教材观察图56的同址计算结构，发现当运算完成后，FFT的输出X(k)按正常顺序排列在存储单元中，即按X(0)，X(1)，X(7)的顺序排列，但是这时输入x(n)却不是按自然顺序存储的，而是按x(0)，x(4)，x(7)的顺序存入存储单元，看起来好像是“混乱无序”的，实际上是有规律的，我们称之为倒位序。造成倒位序的原因是输入x(n)按标号n的偶奇的不断分组。如果n用二进制数表示为(n2n1n0)2(当N=8=23时，二进制为三位)，第一次分组，由图52看出，n为偶数（相当于n的二进制数的最低位n0=0）在上半部分，n为奇数（相当于n的二进制数的最低位n0=1）在下半部分。下一次则根据次最低位n1为“0”或是“1”来分偶奇（而不管原来的子序列是偶序列还是奇序列），如此继续分下去，直到最后N个长度为1的子序列。图59的树状图描述了这种分成偶数子序列和奇数子序列的过程。5.3.3 变址（倒位序）运算30X2024/5/4信息学科立体化教材图59倒位序的形成5.3.3 变址（倒位序）运算31X2024/5/4信息学科立体化教材蝶形、同址、变址运算一般实际运算中，总是先按自然顺序将输入序列存入存储单元，为了得到倒位序的排列，我们通过变址运算来完成。如果输入序列的自然顺序号I用二进制数（例如n2n1n0）表示，则其倒位序J对应的二进制数就是（n0n1n2），这样，在原来自然顺序时应该放x(I)的单元，现在倒位序后应放x(J)。例如,N=8时，x(3)的标号是I=3，它的二进制数是011，倒位序的二进制数是110，即J=6，所以原来存放在x(011)单元的数据现在应该存放在x(110)内。表51列出了N=8时的自然顺序二进制数以及相应的倒位序二进制数。5.3.3 变址（倒位序）运算自然顺序（I）二进制数倒位序二进制数倒位序（J）0123456700000101001110010111011100010001011000110101111104261537表51N=8时的自然顺序二进制数和相应的倒位序二进制数32X2024/5/4信息学科立体化教材由表51可见，自然顺序数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得下一个倒序值；如果最高位是1（JN/2），则要将最高位变为0（JJN/2），次高位加1（J+N/4）。但次高位加1时，同样要判断次高位0、1值，如果为0（JI时，才将原存放x(I)及存放x(J)的存储单元内的内容互换。这样就得到输入所需的倒位序列的顺序。可以看出，其结果与图55的要求是一致的。5.3.3 变址（倒位序）运算图510N=8倒位序的变址处理34X2024/5/4信息学科立体化教材以图56的8点FFT为例，其输入是倒位序的，输出是自然顺序的。其第一级（第一列）每个蝶形的两节点间“距离”为1，第二级每个蝶形的两节点“距离”为2，第三级每个蝶形的两节点“距离”为4。由此类推得，对N=2M点FFT，当输入为倒位序，输出为正常顺序时，其第m级运算，每个蝶形的两节点“距离”为2m1。由于对第m级运算，一个DFT蝶形运算的两节点“距离”为2m1，因而式（518）可写成：5.3.4 蝶形运算两节点的“距离”和Nr的确定为了完成上两式运算，还必须知道系数Nr的变换规律。仔细观察图56的流图可以发现r的变换规律是：把式（522）中蝶形运算两节点中的第一个节点标号值，即k值，表示成M位（注意N=2M）二进制数；把此二进制数乘上2Mm，即将此M位二进制数左移Mm位（注意m是第m级运算），把右边空出的位置补零，此数即为所求r的二进制数。从图56看出，Nr因子最后一列有N/2种，顺序为WN0，WN1，其余可类推。35X2024/5/4信息学科立体化教材显然，对于任何流图，只要保持各节点所连的支路及传输系数不变，则不论节点位置怎么排列所得流图总是等效的，所得最后结果都是x(n)的DFT的正确结果，只是数据的提取和存放的次序不同而已。这样就可得到按时间抽取的FFT算法的若干其他形式流图。将图56中和x(4)水平相连的所有节点和x(1)水平相连的所有节点位置对调，再将和x(6)水平相连的所有节点与和x(3)水平相连的所有节点对调，其余诸节点保持不变，可得图511的流图。图511与图56的蝶形相同，运算量也一样，不同点是：数据存放的方式不同，图56是输入倒位序、输出自然顺序，图511是输入自然顺序、输出倒位序；取用系数的顺序不同,图56的最后一列是按N0,N1,N2,N3的顺序取用系数，且其前一列所用系数是后一列所用系数中具有偶数幂的那些系数（例如N0,N2,.）；图511是最后一列是按N0,N1,N2,N3的顺序取用系数，且其前一列所用系数是后一列所用系数的前一半，这种流图最初由库利和图基给出的时间抽取法。经过简单变换，也可得输入与输出都是按自然顺序排列的流图以及其他各种形式的流图，可以参见其他参考文献。5.3.5 按时间抽取的FFT算法的其他形式流图36X2024/5/4信息学科立体化教材图511时间抽取、输入自然顺序、输出倒位序的FFT流图5.3.5 按时间抽取的FFT算法的其他形式流图37X2024/5/4信息学科立体化教材5.4 按频率抽取（DIF）的基-2 FFT算法 5.4.1算法原理5.4.2按频率抽取法的运算特点38X2024/5/4信息学科立体化教材5.4.1 算法原理这种按频率抽取（DIF）的基2FFT算法推导过程遵循两个规则：对时间进行前后分解；对频率进行偶奇分解。设序列点数为N=2M，M为正整数。同按时间抽取（DIT）的基2FFT算法相同，若不满足这个条件，可以人为地加上若干个零值点，达到这一要求。按规则把输入序列按前一半、后一半分开（不是按偶数、奇数分开），把N点DFT写成两部分。k=0,1,N139X2024/5/4信息学科立体化教材式中，用的是，而不是，因而这并不是N/2点DFT。由于，故，可得k=0,1,N1（523）当k为偶数时，(1)k=1；k为奇数时，(1)k=1。因此，按k的奇偶可将X(k)分为两部分:（524）5.4.1 算法原理40X2024/5/4信息学科立体化教材（525）式（524）为前一半输入与后一半输入之和的N/2点DFT，式（525）为前一半输入与后一半输入之差再与WNn之积的N/2点DFT。令:（527）5.4.1 算法原理41X2024/5/4信息学科立体化教材图512频率抽取法蝶形运算单元由此构成频率抽取法蝶形运算单元，如图512所示。5.4.1 算法原理这样，我们就把一个N点DFT按k的奇偶分解为两个N/2点的DFT了。N=8时，上述分解过程示于图513。与时间抽取法的推导过程一样，由于N=2M，N/2仍是一个偶数，因而可以将每个N/2点DFT的输出再分解为偶数组与奇数组，这就将N/2点DFT进一步分解为两个N/4点DFT。这两个N/4点DFT的输入也是先将N/2点DFT的输入上下对半分开后通过蝶形运算而形成的，图514示出了这一步分解的过程。这样的分解可以一直进行到第M次（N=2M），第M次实际上是做两点DFT，它只有加减运算。然而，为了有统一的运算结构，仍然用一个系数为WN0的蝶形运算来表示，这N/2个两点DFT的N个输出就是x(n)的N点DFT的结果X(k)。图515表示一个N=8完整的按频率抽取的基2FFT运算结构。42X2024/5/4信息学科立体化教材图513按频率抽取的第一次分解5.4.1 算法原理43X2024/5/4信息学科立体化教材图514按频率抽取的第二次分解5.4.1 算法原理44X2024/5/4信息学科立体化教材图515按频率抽取的FFT（N=8）信号流图5.4.1 算法原理45X2024/5/4信息学科立体化教材按频率抽取法的运算特点与时间抽取法基本相同。从图515可以看出，它也是通过(N/2)M个蝶形运算完成的。每一个蝶形结构完成下述基本迭代运算：式中，m表示m列迭代，k,j为整数所在行数，此式的蝶形运算如图516所示，也需要一次复乘和两次复加。5.4.2 按频率抽取法的运算特点图516频率抽取法蝶形运算一般流图形式46X2024/5/4信息学科立体化教材由图515与图56相比较，DIF法与DIT法的区别是：图515的DIF输入是自然顺序，输出是倒位序的，这与图56的DIT法正好相反。但这不是实质性的区别，因为DIF法与DIT法一样，都可将输入或输出进行重排，使二者的输入或输出顺序变成自然顺序或倒位序顺序。DIF的基本蝶形（图516）与DIT的基本蝶形（图57）则有所不同，这才是实质的不同，DIF的复数乘法只出现在减法之后，DIT则是先作复乘后再作加减法。就运算量来说DIF与DIT则是相同的，即都有M级（列）运算，每级运算需N/2个蝶形运算来完成，总共需要mF=(N/2)logN次复乘与aF=NlogN次复加，DIF法与DIT法都可进行原位运算。频率抽取FFT算法的输入是自然顺序，输出是倒位序的。因此运算完毕后，要通过变址计算将倒位序转换成自然位序，然后再输出。转换方法与时间抽取法相同。5.4.2 按频率抽取法的运算特点47X2024/5/4信息学科立体化教材由按时间抽取法与按频率抽取法的基本蝶形（图56与图516）运算看出，如果将DIT的基本蝶形加以转置，就得到DIF的基本蝶形;反过来,将DIF的基本蝶形加以转置，就得到DIT的基本蝶形，因而DIT法与DIF法的基本蝶形是互为转置的。按照转置定理，两个流图的输入输出特性必然相同。转置就是将流图的所有支路方向都反向，并且交换输入与输出，但节点变量值不交换，这样即可从图56得到图516或者从图516得到图56，因而对每一种按时间抽取的FFT流图都存在一个按频率抽取的FFT流图。这样把图55、图511的流图分别加以转置，就可得到不同DIF的FFT流图。因此可以说，有多少种按时间抽取的FFT流图就存在多少种按频率抽取的FFT流图。频率抽取法与时间抽取法是两种等价的FFT运算。5.4.2 按频率抽取法的运算特点48X2024/5/4信息学科立体化教材以上所讨论的FFT的运算方法同样可用于IDFT的运算，称为快速傅里叶反变换，简称为IFFT。由DFT和IDFT的定义（见式（528），将两定义式进行比较，可以看出，只要把DFT运算中的每一个系数WN nk改为WN-nk，并乘以系数1/N，就可以用FFT算法来计算IDFT，也就得到了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 快速傅里叶反变换（528）49X2024/5/4信息学科立体化教材上面这种IFFT算法虽然编程很方便，但要改变FFT的程序和参数才能实现。现介绍另一种计算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 快速傅里叶反变换50X2024/5/4信息学科立体化教材5.6.1利用FFT对信号进行谱分析5.6.2线性卷积的FFT算法5.6 快速傅里叶变换FFT的应用51X2024/5/4信息学科立体化教材所谓谱分析就是计算信号的频谱，包括幅度谱、相位谱和功率谱。1.谱分析步骤 对连续信号进行谱分析的框图如图5-17所示。在图5-17中，前置低通滤波器LPF（预滤波器）的引入，是为了消除或减少时域连续信号转换成序列时可能出现的频谱混叠的影响。在实际工作中，离散时间信号x(n)的时宽是很长的，甚至是无限长的（例如语音或音乐信号）。由于DFT的需要(实际应用FFT计算)，必须把x(n)限制在一定的时间区间之内，即进行数据截断。数据的截断相当于加窗处理（窗函数见7.3节）。因此，在计算FFT之前，用一个时域有限的窗函数w(n)加到x(n)上是非常必要的。图517时域连续信号离散傅里叶分析的处理步骤5.6.1 利用FFT对信号进行谱分析52X2024/5/4信息学科立体化教材 xc(t)通过A/D变换器转换(忽略其幅度量化误差)成取样序列x(n)，其频谱用X(ej)表示，它是频率的周期函数，即(529)式中，Xc(j)或为xc(t)的频谱。5.6.1 利用FFT对信号进行谱分析在实际应用中，前置低通滤波器的阻带不可能是无限衰减的，故由Xc(j)周期延拓得到的X(ej)有非零重叠，即出现频谱混叠现象。由于进行FFT的需要，必须对序列x(n)进行加窗处理，即v(n)=x(n)w(n)，加窗对频域的影响，用卷积表示如下:53X2024/5/4信息学科立体化教材(530)最后是进行FFT运算。加窗后的DFT是0kN1(531)式中，假设窗函数长度L小于或等于DFT长度N，为进行FFT运算，这里选择N为2的整数幂次即N=2m。5.6.1 利用FFT对信号进行谱分析54X2024/5/4信息学科立体化教材有限长序列v(n)=x(n)w(n)的DFT相当于v(n)傅里叶变换的等间隔取样。(532)V(k)便是sc(t)的离散频率函数。因为DFT对应的数字域频率间隔为=2/N，且模拟频率和数字频率间的关系为=，其中=2f。所以，离散的频率函数第k点对应的模拟频率为:(533)(534)5.6.1 利用FFT对信号进行谱分析55X2024/5/4信息学科立体化教材由式（534）很明显可看出，数字域频率间隔=2/N对应的模拟域谱线间距为(535)谱线间距，又称频谱分辨率（单位：Hz）。所谓频谱分辨率是指可分辨两频率的最小间距。它的意思是，如设某频谱分析的F=5Hz，那么信号中频率相差小于5Hz的两个频率分量在此频谱图中就分辨不出来。5.6.1 利用FFT对信号进行谱分析56X2024/5/4信息学科立体化教材总结可见，谱分析中参数主要有：T为取样时间间隔（取样周期），（单位：s）；fs为取样频率（单位：Hz），fs=1/N；fh为信号频率的最高频率分量（单位：Hz）；tp为截取连续时间信号的样本长度（又称最小记录长度，单位：s）；(536)F为谱线间距，又称频率分辨率（单位：Hz）。(537)在实际应用中,要根据信号最高频率fh和频谱分辨率F的要求，来确定T、tp和N的大小。5.6.1 利用FFT对信号进行谱分析57X2024/5/4信息学科立体化教材由以上讨论得到谱分析的步骤如下：（1 1）数据准备）数据准备首先，由取样定理，为保证取样信号不失真，fs2fh(fh为信号频率的最高频率分量，也就是前置低通滤波器阻带的截止频率)，即应使取样周期T满足然后，由频谱分辨率F和T确定N,；为了使用FFT运算，这里选择N为2的幂次即N=2m，由式(537)可知，N大，分辨率好，但会增加样本记录时间tp。最后，由N,T确定最小记录长度，tp=NT。这样在记录长度中进行N点取样就得到要求的离散序列x(n)。5.6.1 利用FFT对信号进行谱分析58X2024/5/4信息学科立体化教材（2）用用FFT计算信号的频谱计算信号的频谱应用前面讨论的FFT算法计算x(n)的频谱，得到X(k)。（3）由由X(k)求幅度谱、相位谱和功率谱求幅度谱、相位谱和功率谱一个长度为N的时域离散序列x(n)，其离散傅里叶变换X(k)（离散频谱）是由实部和虚部组成的复数，即（538）将X(k)写成极坐标形式（539）式中，|X(k)|称为幅度谱，argX(k)称为相频谱。实际中常常用信号的功率谱表示，功率谱是幅频谱的平方，功率谱PSD定义为（540）功率谱具有突出主频率的特性，在分析带有噪声干扰的信号时特别有用。将式（539）绘成的图形称为频谱图。由频谱图可以知道信号存在哪些频率分量，它们就是谱图中峰值对应的点。谱图中最低频率为k=0，对应实际频率为0（即直流）；最高频率为k=N/2，对应实际频率为f=fs/2；对处于0，1，2，N/2上的任意点k，对应的实际频率为f=kF=k fs/N。5.6.1 利用FFT对信号进行谱分析59X2024/5/4信息学科立体化教材例例 5-1有一频谱分析用的FFT处理器，其取样点数必须是2的整数幂，假定没有采用任何特殊的数据处理措施，已给条件为：频率分辨率10Hz；信号最高频率4kHz。试确定以下参量最小记录长度tp、最大取样间隔T（即最小取样频率）、在一个记录中的最少点数N。5.6.1 利用FFT对信号进行谱分析解：解：由分辨率的要求确定最小长度tp从信号的最高频率确定最大可能的取样间隔T（即最小取样频率fs=1/T）。按取样定理即从而最小记录点数N应满足取60X2024/5/4信息学科立体化教材2.谱分析的误差分析谱分析的误差分析利用FFT对连续信号进行傅里叶分析时可能造成一定的误差，从而产生如下现象。（1）频谱混叠失真频谱混叠失真在图517画出的基本步骤中，A/D变换前利用前置低通滤波器进行预滤波，使xc(t)频谱中最高频率分量不超过fh。假设A/D变换器的取样频率为fs，按照奈奎斯特取样定理，为了不产生混叠，必须满足也就是取样间隔T满足5.6.1 利用FFT对信号进行谱分析（541）一般应取fs=(2.53.0)fh如果不满足fs2fh，就会产生频谱混叠失真。对于FFT来说，频率函数也要取样，变成离散的序列，其取样间隔为F(即频率分辨率)。满足61X2024/5/4信息学科立体化教材从以上可以看出，信号的最高频率分量fh与频率分辨率F存在矛盾关系，要想fh增加，则时域取样间隔T就一定减小，从而fs就增加，若是固定N，必然要增加F，即分辨率下降。反之，要提高分辨率（减小F），就要增加tp，当N给定时，必然导致T的增加（fs减小）。要不产生混叠失真，则必然会减小高频容量（信号的最高频率分量）fh。要想兼顾高频容量fh与频率分辨率F，即一个性能提高而另一个性能不变（或也得以提高）的惟一办法就是增加记录长度的点数N，即要满足这个公式是未采用任何特殊数据处理（例如加窗处理）的情况下，为实现基本FFT算法所必须满足的最低条件。如果加窗处理，相当于时域相乘，则频域周期卷积，必然加宽频谱分量，频率分辨率就可能变差，为了保证频率分辨率不变，则须增加数据长度tp。5.6.1 利用FFT对信号进行谱分析62X2024/5/4信息学科立体化教材（2）栅栏效应栅栏效应利用FFT计算频谱，只给出离散点k=2k/N或k=2k/(NT)上的频谱取样值，也就是只给出了基频F的整数倍的频谱值，而不可能得到连续频谱函数，这就像通过一个“栅栏”观看信号频谱，所以只能在离散点上看到信号频谱，这种现象称之为“栅栏效应”,如图518所示。这时，如果在两个离散的谱线之间有一个特别大的频谱分量，就无法检测出来了。减小栅栏效应的一个方法就是要使频域取样更密，即增加频域取样点数N，在不改变时域数据的情况下，必然是在数据末端添加一些零值点，使一个周期内的点数增加，但并不改变原有的记录数据。频谱取样为2k/N，N增加，必然使样点间距更近（单位圆上样点更多），谱线更密，谱线变密后原来看不到的谱分量就有可能看到了。注意，补零以改变计算FFT的周期时，所用窗函数的宽度不能改变。换句话说，必须按照数据记录的原来的实际长度选择窗函数，而不能按照补了零值点后的长度来选择窗函数。补零不能提高频率分辨率，这是因为数据的实际长度仍为补零前的数据长度。5.6.1 利用FFT对信号进行谱分析63X2024/5/4信息学科立体化教材图图 5-18 栅栏效应栅栏效应5.6.1 利用FFT对信号进行谱分析64X2024/5/4信息学科立体化教材（3）频谱泄漏频谱泄漏对信号进行FFT计算，首先必须使其变成有限时宽的信号，这就相当于信号在时域乘一个窗函数如矩形窗，窗内数据并不改变。时域相乘即v(n)=x(n)w(n)，加窗对频域的影响，可用式（530）卷积公式表示卷积的结果，造成所得到的频谱V(ej)与原来的频谱X(ej)不相同，有失真。这种失真最主要的是造成频谱的“扩散”（拖尾、变宽），这就是所谓的“频谱泄漏”。由上可知，泄漏是由于我们截取有限长信号所造成的。对具有单一谱线的正弦波来说，它必须是无限长的。也就是说，如果我们输入信号是无限长的，那么FFT就能计算出完全正确的单一线频谱。可是我们不可能这么做，而只能取有限长记录样本。如果在该有限长记录样本中，正弦信号又不是整数个周期时，就会产生泄漏。5.6.1 利用FFT对信号进行谱分析65X2024/5/4信息学科立体化教材例如，一个周期为N=16的余弦信号x(n)=cos(6n16）截取一个周期长度的信号即x1(n)=cos(6n/16)R16(n),其16点FFT的频谱图见519（a)所示，若截取的长度为13，则其16点FFT的频谱图见519(b)所示。由此可见，频谱不再是单一的谱线，它的能量散布到整个频谱的各处。这种能量散布到其他谱线位置的现象即为“频谱泄漏”。5.6.1 利用FFT对信号进行谱分析图图5-19 16点点FFT的频谱图的频谱图（a）（b）66X2024/5/4信息学科立体化教材应该说明，泄漏也会造成混叠，因为泄漏将会导致频谱的扩展，从而使最高频率有可能超过折叠频率（fs/2），造成频率响应的混叠失真。泄漏造成的后果是降低频谱的分辨率。此外，由于在主谱线两边形成很多旁瓣，引起不同频率分量间的干扰(简称谱间干扰)，特别是强信号谱的旁瓣可能湮没弱信号的主谱线，或者把强信号谱的旁瓣误认为是另一信号的谱线，从而造成假信号，这样就会使谱分析产生较大偏差。在进行FFT运算时，时域截断是必然的，从而频谱泄漏和谱间干扰也是不可避免的。为尽量减小泄漏和谱间干扰的影响，需增加窗的时域宽度(频域主瓣变窄)，但这又导致运算量及存储量的增加；其次，数据不要突然截断，也就是不要加矩形窗，而是加各种缓变的窗（例如：三角形窗、升余弦窗、改进的升余弦窗等），使得窗谱的旁瓣能量更小，卷积后造成的泄漏减小，这个问题在FIR滤波器设计一章（第7章）中将会讨论到。5.6.1 利用FFT对信号进行谱分析67X2024/5/4信息学科立体化教材 例例5-2 已知信号x(t)=0.15sin(2f1t)+sin(2f2t)0.1sin(2f3t)，其中,f1=1Hz，f2=2Hz，f3=3Hz。x(t)包含三个正弦波，但从时域波形图526（a）来看，似乎是一个正弦信号，很难看到小信号的存在，因为它被大信号所掩盖。取fs=32Hz作频谱分析。解解 因fs=32Hz，故该信号为周期信号，其周期为N=32。现对x(n)作32点的离散傅里叶变换(DFT)，其幅度特性|X(k)|如图520(b)所示。图中仅给出了k=0,1,15的结果。k=16,17,32的结果可由|X(Nk)|=|X(k)|得出。因N=32,故频率分辨率F=fs/N=1Hz。由图520(b)可知,k=1,2,3所对应的频谱即为频率f1=1Hz，f2=2Hz,f3=3Hz的正弦波所对应的频谱。幅频图如图520（b）,该图中小信号成分清楚显示出来。可见小信号成分在时域中很难辨识而在频域中容易识别。5.6.1 利用FFT对信号进行谱分析68X2024/5/4信息学科立体化教材图520已知信号的时域图和幅频图5.6.1 利用FFT对信号进行谱分析69X2024/5/4信息学科立体化教材y(n)也是有限长序列，其点数为L+M1点。下面首先讨论直接计算线性卷积的运算量。由于每一个x(n)的输入值都必须和全部的h(n)值相乘一次，因而总共需要LM次乘法，这就是直接计算的乘法次数，以md表示为md=LM用FFT算法也就是用圆周卷积来代替这一线性卷积时，为了不产生混叠，其必要条件是使x(n)，h(n)都补零值点，补到至少N=M+L1，即:（543）0nL1LnN10nM1MnN15.6.2 线性卷积的FFT算法1.FFT的快速卷积法的快速卷积法以FIR滤波器为例，因为它的输出等于有限长单位脉冲响应h(n)与有限长输入信号x(n)的离散线性卷积。设x(n)为L点，h(n)为M点，输出y(n)为70X2024/5/4信息学科立体化教材然后计算圆周卷积N这时，y(n)就能代表线性卷积的结果。用FFT计算y(n)的步骤如下：求H(k)=DFTh(n)，N点；求X(k)=DFTx(n),N点；计算Y(k)=X(k)H(k);求y(n)=IDFTY(k)，N点。5.6.2 线性卷积的FFT算法步骤、都可用FFT算法来完成。此时的工作量如下：三次FFT运算共需要次相乘，还有步骤的N次相乘，因此共需要相乘次数为（544）比较直接计算线性卷积(简称直接法)和FFT法计算线性卷积(简称FFT法)这两种方法的乘法次数。设式（543）与式（544）的比值为Km，则（545）71X2024/5/4信息学科立体化教材分两种情况讨论如下：（1）x(n)与h(n)点数差不多。例如，M=L，则N=2M12M，则这样可得下表：M=L8163264128256512102420484096Km0.572 0.9411.62.785.928.821629.2453.999.9当M=8时，FFT法的运算量大于直接法；当M=32时，二者相当；当M=512时，FFT法运算速度可快16倍；当M=4096时，FFT法约快100倍。可以看出，当M=L且M超过32以后，M越长，FFT法的好处越明显。因而将圆周卷积称为快速卷积。5.6.2 线性卷积的FFT算法72X2024/5/4信息学科立体化教材（2）当x(n)的点数很多时，即当LM。通常不允许等x(n)全部采集齐后再进行卷积;否则，使输出相对于输入有较长的延时。此外，若N=L+M1太大，h(n)必须补很多个零值点，很不经济，且FFT的计算时间也要很长。这时FFT法的优点就表现不出来了，因此需要采用分段卷积或称分段过滤的办法。即将x(n)分成点数和h(n)相仿的段，分别求出每段的卷积结果，然后用一定方式把它们合在一起，便得到总的输出，其中每一段的卷积均采用FFT方法处理。有两种分段卷积的办法:重叠相加法和重叠保留法。5.6.2 线性卷积的FFT算法73X2024/5/4信息学科立体化教材 2重叠相加法重叠相加法设h(n)的点数为M，信号x(n)为很长的序列。我们将x(n)分解为很多段，每段为L点，L选择成和M的数量级相同，用xi(n)表示x(n)的第i段：iLn(i+1)L1其他ni=0,1,（547）则输入序列可表示成（548）5.6.2 线性卷积的FFT算法这样，x(n)和h(n)的线性卷积等于各xi(n)与h(n)的线性卷积之和，即（549）74X2024/5/4信息学科立体化教材每一个xi(n)*h(n)都可用上面讨论的快速卷积办法来运算。由于xi(n)*h(n)为L+M1点，故先对xi(n)及h(n)补零值点，补到N点。为便于利用基2FFT算法，一般取N=2mL+M1，然后作N点的圆周卷积：N由于xi(n)为L点，而yi(n)为（L+M1）点（设N=L+M1），故相邻两段输出序列必然有（M1）个点发生重叠，即前一段的后（M1）个点和后一段的前（M1）个点相重叠，如图521所示。按照式（549），应该将重叠部分相加再和不重叠的部分共同组成输出y(n)。5.6.2 线性卷积的FFT算法和上面的讨论一样，用FFT算法实现重叠相加法的步骤如下：计算N点FFT，H(k)=DFTh(n)；将x(n)分段为xi(n)，计算N点FFT，Xi(k)=DFTxi(n)；相乘，Yi(k)=Xi(k)H(k)；计算N点IFFT，yi(n)=IDFTYi(k)；将各段yi(n)（包括重叠部分）相加得到。由于线性卷积为各输出段的重叠部分相加而得的，故称此法为重叠相加法。75X2024/5/4信息学科立体化教材图521重叠相加法图形76X2024/5/4信息学科立体化教材图521重叠相加法图形77X2024/5/4信息学科立体化教材 3 重叠保留法重叠保留法如果将重叠相加法中分段序列中补零的部分改为保留原来输入序列的值，即将x(n)分为L=NM+1个点的段，在每一段的前边补上前一段保留下来的（M1）个输入序列值，组成L+M1点序列xi(n)，如图522（a）所示。如果L+M11时，随着k的增加螺线内缩;W01时，则随k的增加螺线外伸;W0=1时，表示是半径为A0的一段圆弧。若又有A0=1,则这段圆弧是单位圆的一部分。5.7.1 算法基本原理88X2024/5/4信息学科立体化教材图524螺线取样5.7.1 算法基本原理89X2024/5/4信息学科立体化教材当M=N，A=A0ej0=1，W=W0ej0=（W0=1，0=2/N）这一特殊情况时，各zk就均匀等间隔地分布在单位圆上，这就是求序列的DFT。将式（555）的zk代入变换表达式（554），可得0kM1(559)5.7.1 算法基本原理90X2024/5/4信息学科立体化教材直接计算这一公式，与直接计算DFT相似，总共算出M个取样点，需要NM次复数乘法与(N1)M次复数加法。当N，M很大时，这个量很大，这就限制了运算速度。但是，下面我们将看到，通过一定的变换，以上运算可以转换为卷积形式，从而可以采用FFT算法，这样就可以大大提高运算速度。nk可以用布鲁斯坦所提出的表达式来替换5.7.1 算法基本原理（560）将式（560）代入式（559），可得（561）令（562）（563）n=0,1,N-191X2024/5/4信息学科立体化教材则它们的卷积为式中，k=0,1,M1。式（561）又可以用卷积的形式表示为k=0,1,M15.7.1 算法基本原理（564）（565）由式（565）可看出，如果我们对信号按式（562）先进行一次加权处理，加权系数为;然后，通过一个单位脉冲响应为h(n)的线性系统即求g(n)与h(n)的线性卷积;最后，对该系统的前M点输出再做一次加权，这样就得到了全部M点螺线取样值X(zn)(n=0,1,M1)。这个过程可以用图525表示。从图中我们看到，运算的主要部分是由线性系统来完成的。由于系统的单位脉冲响应可以想象为频率随时间(n)呈线性增长的复指数序列。在雷达系统中,这种信号称为线性调频信号（ChirpSignal），因此，这里的变换称为线性调频Z变换。92X2024/5/4信息学科立体化教材图525ChirpZ变换的线性系统表示5.7.1 算法基本原理93X2024/5/4信息学科立体化教材由式（563）可看出，线性系统h(n)是非因果的，当n的取值为0到N1，k的取值为0,1,M1时，h(n)是在n=(N1)到n=M1取值。也就是说，h(n)是一个有限长序列，点数为N+M1，见图526（a）。输入信号g(n)也是有限长序列，点数为N。g(n)*h(n)的点数为2N+M2，因而用圆周卷积代替线性卷积且不产生混叠失真条件是圆周卷积的点数应大于或等于2N+M2。但是，由于我们只需要前M个值X（zk）(k=0,1,M1)，对以后的其他值是否有混叠失真并不感兴趣，这样可将圆周卷积的点数缩减到最小为N+M1。当然，为了进行基2FFT运算，圆周卷积的点数应取为LN+M1，同时又满足L=2m的最小L。这样可将h(n)先补零值点，补到点数等于L，也就是从n=M开始补L(N+M1)个零值点，补到n=L-N处，或补L(N+M1)个任意序列值，然后将此序列以L为周期而进行周期延拓，再取主值序列，从而得到进行圆周卷积的一个序列，如图526（b）所示。进行圆周卷积的另一个序列只需要将g(n)补上零值点，使之成为L点序列即可，如图526（c）所示。5.7.2 Chirp-Z5.7.2 Chirp-Z变换（变换（CZTCZT）的实现步骤）的实现步骤94X2024/5/4信息学科立体化教材图526ChirpZ变换的圆周卷积图（MnL1时，h(n)和g(n)的圆周卷积不代表线性卷积）5.7.2 Chirp-Z5.7.2 Chirp-Z变换（变换（CZTCZT）的实现步骤）的实现步骤95X2024/5/4信息学科立体化教材归纳起来，CZT运算的实现步骤如下：（1）选择一个最小的整数L，使其满足LN+M1，同时满足L=2m，以便采用基2FFT算法。（2）将（见图526（c）补上零值点，变为L点序列，因而有：0nN10nL1（566）并利用FFT法求此序列的L点DFT0rL1（541）（567）5.7.2 Chirp-Z5.7.2 Chirp-Z变换（变换（CZTCZT）的实现步骤）的实现步骤（3）形成L点序列h(n)，如上所述构造长为L的序列h(n)，0nM1LN+1nL10nM1（568）96X2024/5/4信息学科立体化教材如