CHP4快速傅立叶变换

1、第二部分 傅立叶变换及其快速算法之第四章快速傅里叶变换,目录,4.1概述 4.2时间抽取（DIT）基2FFT算法 4.3频率抽取（DIF）基2FFT算法4.5分裂基算法4.6线性调频Z变换4.7与本章有关节MATLAB文件,4.1 概述,快速傅里叶变换(FFT)是求解离散傅里叶变换(DFT)的快速算法。 问题的提出 直接计算N点DFT需要的计算量是多少？,计算一个X(k)需要N次复数乘法和N一1次复数加法。算出全部N点X(k)共需N2次复数乘法和N(N一1)次复数加法. 总运算量近似地正比于N2 。当N值很大（如2-D图像处理），运算量将非常庞大；同时，对于实时性很强的信号处理来说，必将对计算

2、速度有十分苛刻的要求。为此，需要改进对DFT的计算方法，以减少总的运算次数。,4.1 概述,在正交矩阵中，虽然有N2个元素，但只有N个不同的值，且有些取值特别简单，主要由于旋转因子具有如下的特点：,对称性,周期性,下面以四点DFT为例来说明快速算法的思路。,4.1 概述,4.1 概述,交换矩阵第二列和第三列得,从上面的结果可以看出,利用对称性和周期性，求四点DFT只需要一次复数乘法，称为Coolkey-Tukey算法。,4.1 概述,算法分类： N为2的整次幂：按基数分为基-2FFT算法、基-4FFT算法、混合基FFT算法、分裂基FFT算法;当N不是2的整次幂:典型的有Winograd 算法.

3、 按抽取方法分：时间抽取(DecimationinTime，简称DIT)；频率抽取(DecimationinFrequency，简称DIF),4.1 概述,4.2时间抽取（DIT）基2FFT算法,为了将大点数的DFT 分解为小点数的DFT运算，要求序列的长度N为N2M (M为正整数)。该情况下的变换称为基2 FFT。,核心思想是,问题是如何分最有效？可以对时间变量分 （DIT)，也可对频率变量分(DIF）,4.2时间抽取（DIT）基2FFT算法,基本思路：从时域将N点序列x(n)按奇偶项分解为两组，分别计算两组N/2点DFT，然后再合成一个N点DFT，按此方法继续下去，直到2点DFT，从而减少

4、运算量。算法具体步骤：,一、算法的推导,1、序列x(n)按奇偶项分解为两组，将DFT运算也相应分为两组,则,4.2时间抽取（DIT）基2FFT算法,2、两个N/2点的DFT合成一个N点DFT,问题：A(k)，B(k)都只有N/2个点，怎样得到X(k)的后N/2点？利用周期性和对称性得,4.2时间抽取（DIT）基2FFT算法,4.2时间抽取（DIT）基2FFT算法,4.2时间抽取（DIT）基2FFT算法,3、继续分解（一直分解到两点DFT变换）,A(K)和B(K)仍是高复合数(N2)的DFT，我们可按上述方法继续以分解。令r2l，r2l十1，l0，1，N4-1，则A(K)和B(K)可分别表示为,

5、4.2时间抽取（DIT）基2FFT算法,4.2时间抽取（DIT）基2FFT算法,4.2时间抽取（DIT）基2FFT算法,令,则,4.2时间抽取（DIT）基2FFT算法,4.2时间抽取（DIT）基2FFT算法,同理，令,则,按此方法一直分解下去直到2点DFT，当N=8时，如下：,4.2时间抽取（DIT）基2FFT算法,4.2时间抽取（DIT）基2FFT算法,4.2时间抽取（DIT）基2FFT算法,下面通过讨论寻找FFT的一般规律。,二、算法的讨论,1、“级”的概念 在分解过程中，每分一次，称为一级运算。因为M=log2N，所以N点DFT可以分解为M级，按抽取算法的信号流图中来定义，从左到右分别称

6、为0级、1级到M-1级。,4.2时间抽取（DIT）基2FFT算法,2、蝶形单元 在算法的信号流图中，第m级存在这种运算，这种结构几何形状像蝴蝶，称为蝶形单元,p、q是参于本蝶形单元运算的上、下节点的序号。由于第m级序号的两点只参于这一个蝶形单元的运算，其输出在第m十l级。且这一蝶形单元也不再涉及别的点。由于这一特点，在计算机编程时，我们可将蝶形单元的输出仍放在输入数组中，这一特点称为“同址运算”。,4.2时间抽取（DIT）基2FFT算法,4.2时间抽取（DIT）基2FFT算法,4.2时间抽取（DIT）基2FFT算法,由于一级都含有N/2个蝶形单元，每个蝶形单元需要1次复数乘法和两次复数加法，因

7、此完成log2N级共需要的复数乘法和加法分别为,直接计算DFT时所需的复乘数与复加数都是与N2成正比的。所以采用FFT算法使运算量大大减少。显然，N值愈大，节省的运算量愈多。,4.2时间抽取（DIT）基2FFT算法,3、“组”的概念 在分解过程中，每一级的N/2个蝶形单元可以分成若干组，每一组具有相同的结构和W因子分布。第m级可分成N/2m+1组。,例：N=8=23，分3级。第一级的分组及Wr因子如下： m=0级,分成四组：因子为 m=1级,分成二组,因子为 m=2级,分成一组,因子为,4.2时间抽取（DIT）基2FFT算法,4、Wr因子的分布 由上分析可知,结论：每由后向前（m由M-1-0级

8、）推进一级，则此系数为后级系数中偶数序号的那一半。,4.2时间抽取（DIT）基2FFT算法,5、码位倒置 在FFT算法中，输出的频谱依照正常次序排列，但输入的序列x(n)是按奇偶分开的，分开的规律，以N=8为例，按如下方法进行排序 （1）、将x(n)的序号写成二进制 x(000),x(001), x(110) ，x(111)。 （2）将二进制的码进行翻转，得 x(000),x(100), x(011) ， x(111)。 （3）将二进制的翻转码转换为对应的十进制 x(0),x(4), x(3)，x(7)。 这就是按奇偶抽取得到的顺序。,4.2时间抽取（DIT）基2FFT算法,说明： 在上述的基

9、2FFT算法中，由于每一步分解都是按输入序列x(n)在时域上的次序是属于偶数还是奇数来抽取的，所以称为“按时间抽取法”或“时间抽取”。 上述的基2FFT算法中，抽取也可在频域进行，引出频率抽取（DIF）基2FFT算法。,4.3 频率抽取（DIF）基2FFT算法,设输入序列长度为N=2M(M为正整数)，频率抽取法将输入序列不是按奇、偶分组，而是按前后对半分开，这样可将N点DFT写成前后两部分;将该序列的频域的输出序列X(k)(也是N点序列，按其频域顺序的奇偶分解为越来越短的子序列，称为基2按频率抽取的FFT算法。也称为Sander-Tukey算法。,4.3 频率抽取（DIF）基2FFT算法,算法

10、分析 现将输入x(n)按n的顺序分前后两部分: 前半子序列x(n),0nN/2-1; 后半子序列x(n+N/2),0nN/2-1; 例：N=8时，前半序列为：x(0),x(1),x(2),x(3); 后半序列为： x(4),x(5),x(6),x(7); 考虑N点的DFT,由DFT定义得,4.3 频率抽取（DIF）基2FFT算法,4.3 频率抽取（DIF）基2FFT算法,按k的奇偶将X(k)分成奇偶两部分, k=2r和k=2r+1,考虑k为偶数情况,令,4.3 频率抽取（DIF）基2FFT算法,考虑k为奇数情况,令,4.3 频率抽取（DIF）基2FFT算法,结论 一个N点的DFT被分解为两个N

11、/2点;与时间抽取法的推演过程一样，由于N=2M,因此,N/2仍为偶数，所以可以将N/2点DFT的输出X(k)再分为偶数组和奇数组，这样就将一个N/2点的DFT分成两个N/4点DFT的输入，也是将N/2点的DFT的输入上、下对半分后通过蝶形运算而形成，直至最后为2点DFT。,8点DIF基2FFT算法流图,4.3 频率抽取（DIF）基2FFT算法,4.3 频率抽取（DIF）基2FFT算法,4.3 频率抽取（DIF）基2FFT算法,DIT与DIF的相同之处： （1）DIF与DIT两种算法均为原位运算。 （2）DIF与DIT运算量相同。 DIT与DIF的不同之处： （1）DIF与DIT两种算法结构倒

12、过来。 DIF为输入顺序，输出乱序。运算完毕再运行“二进制倒读”程序。 DIT为输入乱序，输出顺序。先运行“二进制倒读”程序，再进行求DFT。 （2）DIF与DIT根本区别：在于蝶形结不同。 DIT的复数相乘出现在减法之前。 DIF的复数相乘出现在减法之后。,4.5 分裂基算法,自从基2快速算法出现以来，人们仍在不断寻求更快的算法。1984年，法国的杜梅尔(P.Dohamel)和霍尔曼(H.Hollmann)将基2分解和基4分解糅合在一起，提出了分裂基FFT算法。其运算量比前几种算法都有所减少，运算流图却与基2FFT很接近，运算程序也很短。它是目前一种实用的高效新算法。,4.5 分裂基算法,分

13、裂基算法又称基2/4算法,算法的基本思路是:对偶号输出使用基2算法,对奇号输出使用基4算法。 下面首先介绍基4算法： 令N=4M，对N点DFT可按下面方法进行频率抽取,分别令k=4r， k=4r+2， k=4r+1， k=4r+3，其中，r=0，1，2，N/4-1，有,4.5 分裂基算法,4.5 分裂基算法,4.5 分裂基算法,4.5 分裂基算法,算法分析 对于N=4M点DFT，当输出项k为偶数时，采用基2算法，即,当输出项k为奇数时，采用基4算法，即,4.5 分裂基算法,4.5 分裂基算法,4.5 分裂基算法,分析： 一个N点DFT可以分解为一个N/2点DFT和两个N/4点DFT。由x(n)

14、 x(n+N/4) x(n+N/2)和x(n+3N/4)求N/2点DFT和N/4的DFT只需要两次乘法，可以减少运算量。 N/2点DFT可进一步分解为一个N/4点DFT和两个N/8的DFT。N/4的点DFT进一步分解为一个N/8点DFT和两个N/16的DFT。 这样一直下去，直到分解为两点或4点DFT为止。,4.5 分裂基算法,结论： 分裂基FFT算法结构同基2FFT算法结构相似，适用于N=2M的场合，并由M级运算实现。运算流图输入为顺序，输出为倒序。,分裂基FFT算法的计算量,以上提出FFT算法，可以很快地求出全部DFT值。即求出有限长序列x(n)的z变换X(z)在单位园上N个等间隔抽样点z

15、k处的抽样值。它要求N为高度复合数。即N可以分解成一些因子的乘积。例N=2L 实际上： （1）也许对其它围线上z变换取样发生兴趣。如语音处理中，常常需要知道某一围线z变换的极点所处的复频率。 （2）只需要计算单位圆上某一段的频谱,即M不等于N。如窄带信号，希望在窄带频率内频率抽样能够非常密集，提高分辨率，带外则不考虑。 （3）若N是大素数时，不能加以分解，又如何有效计算这种序列DFT。例N=311，若用基2则须补N=28=512点，要补211个零点。,4.6线性调频Z变换,4.6线性调频Z变换,问题提出 为了提高DFT的灵活性，须用新的方法。线性调频z变换(CZT)就是适用这种更为一般情况下，

16、由x(n)求X(zk)的快速变换。 CZT 来自于雷达专业的专用词汇。Z 变换采用螺线抽样,可计算单位圆上任一段曲线的Z变换，适用于更一般情况下（M不等于N）由x(n)求X(zr)的快速算法, 达到频域细化的目的，这种变换称为线性调频Z变换(简称CZT )。,为适应z可以沿平面内更一般的路径取值，我们沿z平面上的一段螺线作等分角的抽样，则z的取样点Zr可表示为：,已 知 N点序列x(n) ，0nN-1,其z变换为,其中M：表示欲分析的复频谱的点数。M不一定等于N。A，W 都为任意复数,令,4.6线性调频Z变换,一、CZT的定义,4.6线性调频Z变换,上式即为CZT的定义.现在讨论A0,W0,0

17、,0的含义: 为输出M点的变换域值.r=0时的A0ej0是CZT的起点,随着r的变化,r0,r1,RM-1构成CZT的变化路径，对于M-1点其极坐标为,4.6线性调频Z变换,4.6线性调频Z变换,CZT在现Z平面上的变换路径是一条螺旋线。,（1）A为起始样点位置,起点半径，大于1时，表示螺旋线在单位圆外，反之，在单位圆内。,起点半相角。,（2）当W01,螺旋线内旋，反之外旋。,（3）当A0=W0=1时，CZT的变换路径为单位圆上的一段弧，起点为P，终点为Q，且M不一定等于N。,4.6线性调频Z变换,4.6线性调频Z变换,考虑A0=W0=1时，在单位圆上CZT，且M不一定等于N。,4.6线性调频

18、Z变换,4.6线性调频Z变换,CZT的线性滤波计算步骤,4.6线性调频Z变换,二、CZT的计算方法 分析：从上面的推导过程可以看出，计算CZT关键是计算一个线性卷积,其中，g(n)应为N点序列，h(n)应为偶对称的无限长序列，考虑到输出M点序列，h(n)的实际长度应为M点。因此，可用DFT来实现两者的卷积，具体步骤如下：,4.6线性调频Z变换,（1）计算并设置g(n),4.6线性调频Z变换,（2）计算并设置h(n),将h(n)设置成长度为L的序列，考虑到其为偶对称序列，且取M点（输出M点序列），取,如下图所示,4.6线性调频Z变换,4.6线性调频Z变换,（3）计算h(n)和g(n)的DFT，得到L点序列H(k)和G(k)。 （4）令Y(k) =H(k) G(k)（乘积）后作Y(k) 的IDFT（反变换）得到时域输出序列y(r) 。 （5）取y(r)的前M点，并乘以W-r2/2,则得最后的输出X(zr)，即,与标准FFT算法相比，CZT算法有以下特点： （1）输入序列长度N及输出序列长度M不需要相等，且N及M不必是高度合成数，二者均可为素数。 （2）Zk的角间隔是任意的，说明其频率分辨率也是任意可控的，角间隔小，分辨率高，反之，分辨率低。 （3）周线不必是z平面上的圆，在语音分析中螺旋周线具有某些优点。 （4）由于起始点z0可任意选定，因此可以从任