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

第四章,快速傅氏变换,引言：,离散傅里叶变换在实际应用中是非常重要的，利用它可以计算信号的频谱、功率谱和线性卷积等。 但是，如果使用定义式来直接计算DFT，当N很大时，即使使用高速计算机，所花的时间也太多。因此，如何提高计算DFT的速度，便成了重要的研究课题。 1965年库利 (Cooley)和图基(Tukey)在前人的研究成果的基础上提出了快速计算DFT的算法，之后，又出现了各种各样快速计算DFT的方法，这些方法统称为快速傅里叶变换(Fast Fourier Transform)，简称为FFT。,快速傅里叶变换(FFT)是离散傅里叶变换(DFT)的快速算法,它是DSP领域中的一项重大突破，它考虑了计算机和数字硬件实现的约束条件、研究了有利于机器操作的运算结构，使DFT的计算时间缩短了12个数量级，还有效地减少了计算所需的存储容量。 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来计算。 设N2M，M为正整数。为推导方便，取N238，即离散时间信号为 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流程图,这样就把一个N点的DFT分解成了两个N/2点的DFT，只是最后求和的符号不同。,式4.1与4.2的运算可用蝶形信号流图来表示,图1 N点的DFT分解成两个N/2点DFT的的时间抽选流程图N=8,按照式4.1和式4.2可画出下图所示的信号流程图。,式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),图2 N/2点的DFT分解成两个N/4点DFT的的时间抽选流程图N=8,这样 ，用式(4.3)(4.6)就可计算图1中的两组N/2点DFT。,图2所示的是其中一组G(k)的计算。,将图1与图2合并，便得到图3所示的信号流程图。,G(0) G(1) G(2) G(3),H(0) H(1) H(2) H(3),将2点DFT的信号流程图移入图3，得到图4所示的8点时间抽选的完整的FFT流程图。,2点的DFT流程图,N=8，图3中N/4点的DFT就是2点的DFT。,图4 基2时间抽选法8点FFT流程图,从图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个蝶形组成。,下图表示了蝶形的一般形式表示。其输入和输出之间满足下列关系：,大多数情况下复数乘法所花的时间最多，因此下面仅以复数乘法的计算次数为例来与直接计算进行比较。,直接计算DFT需要的乘法次数为aD=N2，于是有,即直接计算DFT所需复数乘法次数约为FFT的205倍。 显然，N越大，FFT的速度优势越大。,复数乘法次数：,例：当N=1024时，,因此，完成N点的时间抽选FFT计算的总运算量为,复数加法次数：,不同N值所对应的两种计算方法的复数乘法次数和它们的比值,三、同址(原位)计算,图4包含M级迭代运算，每级由N/2个蝶形计算构成。蝶形计算的优点是可进行所谓同址(原位)计算。 首先每一个蝶形运算都需要两个输入数据，计算结果也是两个数据，与其它结点的数据无关，其它蝶形运算也与这两结点的数据无关。 因此，一个蝶形运算一旦计算完毕，原输入数据便失效了。这就意味着输出数据可以立即使用原输入数据结点所占用的内存。原来的数据也就消失了。这样输出、输入数据利用同一内存单元的这种蝶形计算称为同址计算。,类似地，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个存储存储单元,四、变址计算 （倒序重排）,从图4所示的流程图看出，同址计算要求输入x(n)是“混序”排列的。 所谓输入为“混序”，并不是说输入是杂乱无章的，实际上它是有规律的。 如果输入x(n)的序号用二进制码来表示，就可以发现输入的顺序恰好是正序输入的“码位倒置”。 下表2列出了这种规律。,表2,存储单元 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表示码位倒置的标号,当ln时，x(n)和x(l)不必互相调换。 当ln时， 必须将x(l)和x(n)互相调换，但只能调换一次，为此必须规定每当ln时，要将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流程图,图6 基2时间抽选法8点FFT流程图（输入正序）,另一种形式的流程图是将节点排列成输入和输出两者都是正序排列，但这类流程图不能进行同址计算，因而需要两列长度为N的复数存储器。,其中 r=0, 1, , N/2-1,4.2 基2频率抽选法（桑德图基算法）,基2频率抽选法是在频域内将X(k)逐次分解成奇偶序列，再进行DFT计算 设N2M，M为正整数。为推导方便，取N=238。 将频率偶奇分， 即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所示。,图7 N点的DFT分解成两个N/2点DFT的的频率抽选流程图N=8,时间抽选 蝶形单元,频率抽选 蝶形单元,与时间抽选的FFT算法一样，图7所示的流程图的基本运算也是蝶形运算，但是频率抽选与时间抽选中的蝶形单元有所不同。,N是2的整数幂，所以N/2仍然是偶数。这样可以将N/2点DFT的输出再分为偶数组和奇数组，即将N/2点的DFT计算进一步分解为两个N/4点的DFT计算，以此类推直到不能分解。如图8,图8 基2频率抽选法8点FFT流程图,频率抽选与时间抽选的异同点：,比较图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,根据上式可画出蝶形图,以此类推可求出x(n)的各点，反变换流程如图9,根据式4.1与4.2可得到,图9 8点IFFT流程图,比较上面两式，可以看出，只要把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求IFFT,FFT算法同样可以应用于IDFT的计算。 已知DFT和IDFT公式为,图10 8点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 且MN 式中,如图：,一、算法基本原理,已知x(n)长度为N，则,W0、A0为正实数,由图可知：,1、A0表示起始抽样点Z0的矢量半径长度，通常A0 1，否则Z0将在单位圆之外； 2、q0表示起始抽样点Z0相角，可正可负; 3、f0表示两相邻抽样点间角度差， f0为正时，表示路径逆时针， f0为负时表示路径顺时针。,4、W0表示螺线的伸展率，W01表示随k增大螺线内缩，W01表示随k增大螺线外伸。W0=1表示半径为A0的一段圆弧，又如果A0=1,则圆弧是单位圆一部分。,与直接计算DFT相似，当N、M很大时，运算量很大； 于是可采用布鲁斯坦（Bluestein）等式，将其转换为卷积和形式，从而可采用FFT算法。,Bluestein等式,令,可得：,0kM-1,将Zk代入X(z)，有,0kM-1,计算过程如图：,令,则,或,可想象为频率随时间n成线性增长的复指数序列，在雷达系统中，这种信号称为线性调频信号(Chirp signal)，因此在这里称为线性调频Z变换算法。,序列,0kM-1,0kM-1,0kM-1,0kM-1,由,可知h(n)的长度为n=-(N-1)M-1 点数为L=N+M-1，y(n)的长度为N，要利用L点循环卷积来计算，y(n)需补M-1个零。,二、实现步骤,考虑h(n)的长度 ： 无限长,然后计算h(n)与y(n)的L点循环卷积 其中g(0)到g(M-1)这M个值正与所需要的线性卷积值相等，后N-1个值不需要。 当L=M+N-1=2m时，可利用基2FFT算法。,而h(n)选取值如图：,2、将,3、形成L点h(n)，并利用FFT法求H(k),0nM-1,MnL-N,L-N+1nL-1,CZT运算步骤：,1、选择L值使LM+N-1，且L=2m,补零，变为L点序列。