《快速傅里叶变换》PPT课件.ppt

第4章 快速傅里叶变换,4.1 引言 4.2 直接计算DFT的问题及改进的途径 4.3 按时间抽取（DIT）的基2-FFT算法 4.4 按频率抽取（DIF）的基2-FFT算法 4.5 离散傅里叶反变换（IDFT）的快速计算方法 4.10 线性卷积的FFT算法,4.1 引 言,快速傅里叶变换（FFT）并不是一种新的变换， 而是离散傅里叶变换（DFT）的一种快速算法。 由于有限长序列在其频域也可离散化为有限长序列（DFT），因此离散傅里叶变换（DFT）在数字信号处理中是非常有用的。例如，在信号的频谱分析、 系统的分析、 设计和实现中都会用到DFT的计算。 但是，在相当长的时间里， 由于DFT的计算量太大，即使采用计算机也很难对问题进行实时处理，所以并没有得到真正的运用。 直到1965年首次发现了DFT运算的一种快速算法以后，情况才发生了根本的变化。人们开始认识到DFT运算的一些内在规律，从而很快地发展和完善了一套高速有效的运算方法， 这就是现在人们普遍称之为快速傅里叶变换（FFT）的算法。 FFT出现后使DFT的运算大大简化，运算时间一般可缩短一二个数量级之多，从而使DFT的运算在实际中真正得到了广泛的应用。,4.2 直接计算DFT的问题及改进的途径,一、直接计算DFT的运算量 设x(n)为N点有限长序列，其DFT为,k=0, 1, , N-1,（4-1）,反变换（IDFT）为,n=0, 1, , N-1,（4-2）,二者的差别只在于WN的指数符号不同，以及差一个常数乘因子1/N，所以IDFT与DFT具有相同的运算工作量。 下面我们只讨论DFT的运算量。 一般来说，x(n)和WNnk都是复数，X(k)也是复数，因此每计算一个X(k)值，需要N次复数乘法和N-1次复数加法。而X(k)一共有N个点（k从0取到N-1），所以完成整个DFT运算总共需要N2次复数乘法及N(N-1)次复数加法。 在这些运算中乘法运算要比加法运算复杂，需要的运算时间也多一些。因为复数运算实际上是由实数运算来完成的， 这时DFT运算式可写成,（4-3）,由此可见，一次复数乘法需用四次实数乘法和二次实数加法；一次复数加法需二次实数加法。 因而每运算一个X(k)需4N次实数乘法和2N+2(N-1)=2(2N-1)次实数加法。 所以，整个DFT运算总共需要4N2次实数乘法和2N(2N-1)次实数加法。,当然，上述统计与实际需要的运算次数稍有出入，因为某些WNnk可能是1或j，就不必相乘了，例如W0N=1，W N=-1， WNN/4=-j等就不需乘法。 但是为了便于和其他运算方法作比较， 一般都不考虑这些特殊情况，而是把WNnk都看成复数，当N很大时，这种特例的影响很小。 从上面的统计可以看到，直接计算DFT，乘法次数和加法次数都是和N2成正比的，当N很大时，运算量是很可观的，有时是无法忍受的。,例3-1 根据式（3-1），对一幅NN点的二维图像进行DFT变换，如用每秒可做10万次复数乘法的计算机，当N=1024时，问需要多少时间（不考虑加法运算时间）？ 解 直接计算DFT所需复乘次数为（N2)21012次，因此用每秒可做10万次复数乘法的计算机，则需要近3000小时。 这对实时性很强的信号处理来说，要么提高计算速度，而这样，对计算速度的要求太高了。另外，只能通过改进对DFT的计算方法，以大大减少运算次数。,二、 改善途径 能否减少运算量，从而缩短计算时间呢？仔细观察DFT的运算就可看出， 利用系数Wnk的以下固有特性，就可减少运算量： （1） WNnk的共轭对称性,（2） WNnk的周期性,（3） WNnk的可约性,另外,这样，利用这些特性，使DFT运算中有些项可以合并，并能使DFT分解为更少点数的DFT运算。由于DFT的运算量是与N2成正比的，所以N越小越有利，因而小点数序列的DFT比大点数序列的DFT的运算量要小。 快速傅里叶变换算法正是基于这样的基本思想而发展起来的。它的算法形式有很多种，但基本上可以分成两大类： 按时间抽取(ecimation inTime，缩写为DIT)法 按频率抽取(Decimationin Frequency，缩写为DIF)法,4.3 按时间抽取(DIT)的基-2 FFT算法 (库利-图基算法),一、 算法原理 设序列x(n)长度为N，且满足N=2L，L为正整数。按n的奇偶把x(n)分解为两个N/2点的子序列：,（4-4）,则可将DFT化为,由于 ， 则上式可表示成,（4-5）,式中，X1(k)与X2(k)分别是x1(r)及x2(r)的N/2点DFT：,（4-6）,（4-7）,X1(k)，X2(k)只有N/2个点，即k=0, 1, , N/2-1。 而X(k)却有N个点，即k=0, 1, , N-1， 故用式（4-5）计算得到的只是X(k)的前一半（ ）的结果。,所以,（4-8）,同理可得,（4-9）,式（4-8）、式（4-9）说明了后半部分k值（N/2kN-1）所对应的X1(k),2(k)分别等于前半部分k值（0kN/2-1）所对应的X1(k), X2(k)。,（周期性）,由于,再考虑到WkN 的以下性质:,这样，把式（4-8）、式（4-9）、式（4-10）代入式（4-5），就可将X(k)表达为前后两部分：,（4-10）,（4-11）,（4-12）,图 4-1 时间抽选法蝶形运算流图符号,.蝶形运算,图 4-2 按时间抽选将一个N点DFT分解 为两个N/2点DFT（N=8）,(1)N/2点的DFT运算量: 复乘次数: 复加次数: (2)两个N/2点的DFT运算量:复乘次数: 复加次数: (3)N/2个蝶形运算的运算量:复乘次数: 复加次数:,运算量,复乘:,复加:,总共运算量:,*N点DFT的复乘为N2 ;复加N(N-1);与分解后相比可知,计算工作点差不多减少 一半。,由于N=2L，因而N/2仍是偶数，可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。,（4-13）,例如，n为偶数时的 N/2点分解为:,且,式中：,（4-14）,（4-15）,图4-3给出N=8时，将一个N/2点DFT分解成两个N/4点DFT， 由这两个N/4点DFT组合成一个N/2点DFT的流图。,图 4-3 N/2点DFT分解为两个N/4点DFT,X2(k)也可进行同样的分解：,式中：,同样对n为奇数时，N/2点分为两个N/4点 的序列进行DFT,图 4-4 一个N=8点DFT分解为四个N/4点DFT,根据上面同样的分析知道，利用四个N/4点的DFT及两级蝶形组合运算来计算N点DFT，比只用一次分解蝶形组合方式的计算量又减少了大约一半。,对于N=8时的DFT, N/4点即为两点DFT,因此,式中， , 故上式不需要乘法。,可得：,这种方法的每一步分解，都是按输入序列在时间上的次序是属于偶数还是属于奇数来分解为两个更短的子序列，所以称为“按时间抽取法”。,图 4-5 N=8 按时间抽取的FFT运算流图,二、 运算量 由按时间抽取法FFT的流图可见，当N=2L时，共有L级蝶形， 每级都由N/2个蝶形运算组成，每个蝶形需要一次复乘、 二次复加，因而每级运算都需N/2次复乘和N次复加，这样L级运算总共需要,由于计算机上乘法运算所需的时间比加法运算所需的时间多得多，故以乘法为例，直接DFT复数乘法次数是N2，FFT复数乘法次数是(N/2) log2N。 直接计算DFT与FFT算法的计算量之比为,当N=2048时，这一比值为372.4，即直接计算DFT的运算量是FFT运算量的372.4倍。当点数N越大时，FFT的优点更为明显(见书中P150.表4-1)。,（4-20）,三、按时间抽取的FFT算法的特点,1. 原位运算（同址运算） 从图4-5可以看出这种运算是很有规律的，其每级（每列）计算都是由N/2 个蝶形运算构成的，每一个蝶形结构完成下述基本迭代运算：,式中，m表示第m列迭代，k, j为数据所在行数。式（4-21）的蝶形运算如图4-7所示，由一次复乘和两次复加（减）组成。,图 4-7 按时间抽选的蝶形运算结构,在某列进行蝶形运算的任意两个节点(行)k和j的节点变量 就完全可以确定蝶形运算的结果 ,与其它行(节点)无关。 这样,蝶形运算的两个输出值仍可放回蝶形运算的两个输入所在的存储器中,即实现所谓原位运算。每一组(列)有N/2个蝶形运算,所以只需N个存储单元，可以节省存储单元。,由图4-5的流图看出，某一列的任何两个节点k和j的节点变量进行蝶形运算后，得到结果为下一列k, j两节点的节点变量，而和其他节点变量无关，因而可以采用原位运算，即某一列的N个数据送到存储器后，经蝶形运算，其结果为下一列数据，它们以蝶形为单位仍存储在这同一组存储器中，直到最后输出，中间无需其他存储器。也就是蝶形的两个输出值仍放回蝶形的两个输入所在的存储器中。每列的N/2 个蝶形运算全部完成后，再开始下一列的蝶形运算。这样存储器数据只需N个存储单元。下一级的运算仍采用这种原位方式，只不过进入蝶形结的组合关系有所不同。 这种原位运算结构可以节省存储单元，降低设备成本。,2. 倒位序规律 FFT的输出X(k)按正常顺序排列在存储单元中，即按X(0)，X(1)，X(7)的顺序排列，但是这时输入x(n)却不是按自然顺序存储的，而是按x(0)，x(4)， ， x(7)的顺序存入存储单元，看起来好像是“混乱无序”的，实际上是有规律的，我们称之为倒位序。,这是由奇偶分组造成的,以N=8为例 说明如下:,造成倒位序的原因是输入x(n)按标号n的偶奇的不断分组。 如果n用二进制数表示为(n2n1n0)2（当N=8=23时，二进制为三位）， 第一次分组，由图4-2看出，n为偶数（相当于n的二进制数的最低位n0=0）在上半部分，n为奇数（相当于n的二进制数的最低位 n0=1）在下半部分。 下一次则根据次最低位n1为“0”或是“1”来分偶奇（而不管原来的子序列是偶序列还是奇序列）， 如此继续分下去，直到最后N个长度为1的子序列。 图4-8的树状图描述了这种分成偶数子序列和奇数子序列的过程。,图4-8 描述倒位序的树状图,3.倒位序实现 输入序列先按自然顺序存入存储单元,然后经变址运算来实现倒位序排列。 设输入序列的序号为n,二进制为(n2 n1 n0 )2 , 倒位序顺序用 表示,其倒位序二进制为(n0 n1 n2 )2 , 例如 ，N=8时如下表：,表4-2 码位的倒位序（N=8）,存储单元,自然顺序,变址,倒位序,图 4-9 倒位序的变址处理 (N=8),变址处理方法,4. 蝶形运算两节点的“距离” 以图4-5的8点FFT为例，其输入是倒位序的，输出是自然顺序的。 其第一级（第一列）每个蝶形的两节点间“距离”为1， 第二级每个蝶形的两节点“距离”为2， 第三级每个蝶形的两节点“距离”为4。 由此类推得，对N=2L点FFT，当输入为倒位序， 输出为正常顺序时，其第m级运算，每个蝶形的两节点“距离”为2m-1。,5. WNr的确定 由于对第m级运算，一个DFT蝶形运算的两节点“距离”为2m-1， 因而式（4-21）可写成:,（4-22a）,（4-22b）,为了完成上两式运算，还必须知道系数Nr的变换规律： 把式（4-22）中蝶形运算两节点中的第一个节点标号值， 即k值，表示成L位（注意N=2L）二进制数；, 把此二进制数乘上2L-m，即将此L位二进制数左移M-m位（注意m是第m级运算），把右边空出的位置补零， 此数即为所求r的二进制数。 从图4-5看出，WNr因子最后一列有N/2种，顺序为WN0， WN1， ， 其余可类推。,6.存储单元 存输入序列 需N个存储单元 存放系数 需N/2个存储单元； 共计需(N+N/2)个存储单元。,四、 按时间抽取的FFT算法的其他形式流图 显然，对于任何流图，只要保持各节点所连的支路及传输系数不变，则不论节点位置怎么排列所得流图总是等效的，所得最后结果都是x(n)的DFT的正确结果，只是数据的提取和存放的次序不同而已。这样就可得到按时间抽取的FFT算法的若干其他形式流图。 将图4-5中和x(4)水平相连的所有节点和x(1)水平相连的所有节点位置对调，再将和x(6)水平相连的所有节点与和x(3)水平相连的所有节点对调，其余诸节点保持不变，可得图4-11的流图。 图4-11与图4-5的蝶形相同，运算量也一样，不同点是：, 数据存放的方式不同，图4-5是输入倒位序、输出自然顺序，图4-11是输入自然顺序、输出倒位序； 取用系数的顺序不同，图4-5的最后一列是按 的顺序取用系数，且其前一列所用系数是后一列所用系数中具有偶数幂的那些系数（例如 ）；图4-11是最后一列是按 的顺序取用系数，且其前一列所用系数是后一列所用系数的前一半， 这种流图是最初由库利和图基给出的时间抽取法。 经过简单变换，也可得输入与输出都是按自然顺序排列的流图以及其他各种形式的流图 。,图 4-10 时间抽取、 输入自然顺序、 输出倒位序的FFT流图,4.4 按频率抽取（DIF）的基 -2 FFT算法（桑德-图基算法）,一、 算法原理 仍设序列点数为N=2L，L为正整数。在把输出X(k)按k的奇偶分组之前，先把输入序列按前一半、后一半分开（不是按偶数、 奇数分开）， 把N点DFT写成两部分。,k=0, 1, , N-1,由于 ，故 ，可得,k=0, 1, , N-1,当k为偶数时，(-1)k=1；k为奇数时，(-1)k=-1。因此，按k的奇偶可将X(k)分为两部分:,为前一半输入与后一半输入之和的N/2点DFT,为前一半输入与后一半输入之差再与WNn之积的N/2点DFT。,图 4-14 按频率抽取蝶形运算流图符号,这样可把一个N点DFT按k的奇偶分解为两个N/2点的DFT。 N=8时，上述分解过程示于图4-15。,图 4-15 按频率抽取，将N点DFT分解为两个N/2点DFT的组合,由于N=2L，N/2仍是一个偶数，因而可以将每个N/2点DFT的输出再分解为偶数组与奇数组，这就将N/2点DFT进一步分解为两个N/4 点DFT。 图4-16示出了这一步分解的过程。,图 4-16 按频率抽取的第二次分解,这样的分解可以一直进行到第L次（N=2L），第L次实际上是做两点DFT，它只有加减运算。 然而，为了有统一的运算结构，仍然用一个系数为WN0的蝶形运算来表示， 这N/2个两点DFT的N个输出就是x(n)的N点DFT的结果X(k)。 图4-16表示一个N=8 完整的按频率抽取的基-2 FFT运算结构。,图 4-16 按频率抽取的FFT流图 （N=8）,4.5 离散傅里叶反变换的快速计算方法(IFFT),一.稍微变动FFT程序和参数可实现IFFT,比较两式可知,只要DFT的每个系数 换成 ,最后再乘 以常数1/N就可以得到IDFT的快速算法-IFFT。 另外,可以将常数1/N分配到每级运算中, , 也就是每级蝶形运算均乘以1/2。,二.不改(FFT)的程序直接实现IFFT,这就是说,先将X(k)取共轭,即将X(k)的虚部乘-1，直接利用FFT程序计算DFT；然后再取一次共轭；最后再乘1/N,即得x(n)。所以FFT,IFFT可用一个子程序。,一.线性卷积的长度 设一离散线性移不变系统的冲激响应为 ,其输入信号为 . 其输出为 . 并且 的长度为L点, 的 长度为M点,则：,4.10 线性卷积的FFT算法,二.用FFT算 的步骤：,流程图,三.几点说明 1. 当 M=L 时,用圆周卷积计算线性 卷积的速度快,点数越多,速度越快, N64时,速度增加明显. 2. LM 时,速度增加不太明显,为了 增加速度,采用 (1)重叠相加法 (2) 重叠保留法(从略).,