1、第四章 快速傅里叶变换（FFT）,快速傅里叶变换 （FFT）,快速傅里叶变换（FFT）是计算DFT的一种快速高效的算法。有限长序列在数字技术中占有很重要的地位，主要原因是由于其频谱可以离散化。有限长序列的DFT本身可以完全表达序列的频谱，所以DFT也可以直接对信号进行频谱分析。谱分析在通信技术、图象传输、语音压缩、生物医学等领域都得到应用。,虽然频谱分析和DFT运算很重要，但在很长一段时间里，由于DFT运算量太大，因此它并没有得到真正的运用，频谱分析也大多采用模拟信号滤波的方法解决。直到1965年美国IBM公司的库利和图基这两位科学家首次提出DFT运算的一种快速算法以后，情况才发生了根本变化，
2、人们开始认识到DFT运算的一些内在规律，从而很快地发展和完善了一套高速有效的运算方法快速付里叶变换(FFT)算法。FFT的出现，使DFT的运算大大简化，运算时间缩短1-2个数量级，FFT的问世使得DFT在实际中得到广泛应用。,4.1 DFT运算的特点：,有限长序列x(n)进行一次DFT运算所需的运算量:,一. DFT的运算量：,一般x(n)和 都是复数，X(k)也是复数，所以计算一次DFT的运算量是：,在这些运算中，乘法比加法运算复杂，尤其是复数相 乘，每个复数相乘包括4个实数相乘和2个实数相加，例：又因每个复数相加包括2个实数相加,所以,每计算一个 X(k)要进行4N次实数相乘和 2N+2(
3、N-1)=2(2N-1)次实数相 加，因此，整个DFT运算需要4N2实数相乘和2N(2N-1)次实 数相加。,计算一个X(k)值需： N次复数相乘和(N-1)次复数相加,计算N点X(k)值需： 次复数相乘和 N(N-1)次复数相加,从上面的分析看到，在DFT计算中，不论是乘法和加法，运算量均与N2成正比。因此，N较大时，运算量十分可观。例：计算N=10点的DFT，需要100次复数相乘，而N=1024点时，需要1048576（一百多万）次复数乘法，如果要求实时处理，则要求有很高的计算速度才能完成上述计算量。反变换IDFT与DFT的运算结构相同，只是多乘一个常数1/N，所以二者的计算量相同。,考察
4、DFT的运算特点发现,利用以下两个特性可减少运算量： 1）系数 是一个周期函数，利用它的周期性和对称性可改进运算，提高计算效率。,二. DFT的运算特点,周期性,对称性,我们利用系数 的周期性和对称性，考察它是如何简化DFT运算的过程。,以N4为例， 的矩阵形式为：,DFT的矩阵运算,2) 因为DFT的计算量正比于N2，N小计算量也就小。因此可以把长度为N点的大点数的DFT运算依次分解为若干个小点数的DFT来运算。FFT算法正是基于以上两点基本思想来提高DFT的 运算速度。FFT算法基本上可分为两大类：按时间抽取FFT算法和按频率抽取FFT算法。,结论：,1) 利用系数 的周期性和对称性可以提
5、高DFT的 运算速度。上例中作一次DFT需 N2=16次乘法运 算，而FFT只需6次乘法运算。,4.2 按时间抽取的FFT,为了保证N点DFT的运算分解，首先假设N是2的整幂次方，即 N=2M ，M是正整数。 一.算法原理(1)将序列 按序号n的奇、偶分解为两个 点子序列,DFT运算也相应分为两组：,故：,显然， 和 都只含有N/2点DFT, 所以X(k) 也只包含k=0,1,N/2-1点的DFT，还必须利用系数的周期性和对称性，求出X(k)的后N/2点的DFT。,(2)求：,可见，一个N点的DFT被分解为前后两个N/2点的DFT，这两个N/2点的DFT再合成为一个N点DFT。,同理：,再利用
6、W系数的对称性：,以上两个表达式说明：只要我们计算出N/2个点 的所有 和 的值，就等于求出N点的X(k)值, 这使得DFT的运算量减少了约一半。上述两个公式的运算过程可用专用蝶形图表示：,图(a)需要两次乘法、两次加减法,将图(a)简化成图(b)，此时仅需一次乘法、两次加减法。,按时间抽取的蝶形运算流图：,(3) 和 可以继续分解下去，可将N/2点的 子序列再按奇、偶项分解，一直到最后分解成 两两点的DFT为止。,偶序列中的偶序列,偶序列中的奇序列,奇序列中的偶序列,奇序列中的奇序列,四个 点 DFT,例：N=23=8 按时间抽取的DFT分解过程,由于N/2 仍是偶数，可以被整除，因此可以对
7、两个N/2点的DFT再分别作进一步的分解。将一个点的DFT可以分解成四个点的DFT,直到最后得到两两点的DFT为止。,由于这种方法每一步分解都是按输入序列是属于偶数还是奇数来抽取的，所以称为“按时间抽取的FFT算法”。,N8点按时间抽取的FFT运算流图,第一级蝶形,第二级蝶形,第三级蝶形,时间抽取法FFT的运算特点：,（1）蝶形运算,（2）原位运算结构,（3）码位倒置变换,（4）蝶形类型随迭代次数成倍增加,（1）蝶形运算,对于N=2M，总可以通过M次分解最后分解成2点的DFT运算。这构成从x(n)到X(k)的M级运算过程。从上面流图可看到，每一级运算都由N/2个蝶形运算构成。因此每一级运算都需
8、要N/2次复乘和N次复加，经过M级运算后总共需要的运算量为： 复乘： 复加：而直接进行DFT运算时则与N2 成正比。,算法的运算速度,例: N=2048，N2=4194304 ， ，显然,FFT要比直接DFT运算快得多。,（2）原位运算结构（同址运算）,FFT运算结构很有规律，它具有强烈的对称性，它的所有运算都可以由左下图所示的基本运算构成。 由于运算结构规律，所以硬件实现起来简单，软件编程方便。 FFT的蝶形运算结构很经济，它可以采用原位运算。原位运算 指当数据输入到存储器中以后， 每一级运算的结果仍然储存在同一组存储器中，直到最后输出， 中间无需其它任何存储器，称为原位运算。,例如：输入数
9、据A、B相应的存放在寄存器A(1)和A(2)中，当进行完蝶形运算后，输出和存放在A(1)中，输出差存放在A(2)中。下一级蝶形运算的输入数据又分别存放在A(1)和A(2)中，而把前一级输出的结果覆盖掉，直到最后输出，中间无需任何其它存储器 。原位运算结构节省存储单元，可以降低设备成本，还可节省寻址的时间。,（3）码位倒置变换对按时间抽取FFT的原位运算结构，如果输出按顺序输出,即顺序存放x(0),x(1),x(2),x(7), 但它的输入x(n)却不是按自然顺序输入的,而是按x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)存入 ，这种顺序看起来很杂乱，实际上它很有
10、规律。当用二进制表示这个顺序时，它正好是“码位倒置”的顺序。,例如：顺序输出,码位倒置指序号用二进制码形式表示，然后将二进制码的首尾颠倒。,输入,在实际运算中，输入数据总是按自然顺序输入到存储单元，再通过变址运算将自然顺序转换成码位倒置顺序存储，然后再进行FFT的原位计算。目前有许多通用DSP芯片支持这种码位倒置的寻址功能。,下表中以N=8为例给出码位倒置顺序,I,J,码位倒置的规律(1) 当I=J时，A(I)中存储的数不变，A(I)=A(J);(2) 当IJ 时，意味着A(I)和A(J)的内容已经互换 过了，为了避免将已互换过的数据再次调换，所以A(I)和A(J)的内容不变。,N8点按时间抽
11、取的FFT运算流图,第一级蝶形,第二级蝶形,第三级蝶形,观察N=23=8点FFT的蝶形系数 ： 第一级：有一种类型的蝶形运算系数 第二级: 有二种类型的蝶形运算系数 、 第三级: 有四种类型的蝶形运算系数、 、 、 第L级: 有 种蝶形运算系数,（4）蝶形图的系数,将x(n)按n的顺 序前后对半分开：,4.3 按频率抽取的FFT,频率抽取法是N=2M情况下的另外一种FFT算法,按时间抽取FFT算法的基本思路,按频率抽取FFT算法的基本思路,x(n)按奇、偶一级级分开， X(k)按前一半、后一半一级级分开。,X(k)按奇、偶一级级分开， x(n)按前后对半一次次分开。,一.按频率抽取法的FFT算
12、法原理,按k的奇、偶把X(k)进一步分解为两部分：,a) k为偶数，令k=2r,b) k为奇数，令k=2r+1,按频率抽取的蝶形运算流图：,例：以N=8 按频率抽取的FFT,把8点DFT分解成两个4点的DFT,K为偶数 X(2r),K为奇数 X(2r1),将两个4点DFT再分解成四个2点DFT,由于N/2仍是偶数可以继续分解，直到分解为两两点的DFT为止。,下图是N=8按频率抽取的FFT流图,频率抽取法的分解过程： (1)把输入序列前后对半分； (2)输出按k值的奇、偶将X(k)分解为奇数组和偶数组两部分；(3)输入如果按顺序输入，输出则是码位倒置形式。,按频率抽取法FFT的运算特点： （1）
13、蝶形运算的运算量其运算量与按时间抽取运算量相同。即：复乘： 复加： （2）输入如果是自然顺序，输出则是码位倒置形式。如果要求输出为自然顺序，则需要进行变址运算。按频率抽取的FFT流图可以由按时域抽取的FFT流图转置 而得到。转置指将流图中所有支路的方向反 ,并将输入变输出，输出变输入。,频率抽取法的基本蝶形图与时间抽取法的基本蝶形图有所不同。时间抽取法是先乘后加、减，频率抽取法是先加、减，再乘系数,按频率抽取的 蝶形运算流图,按时间抽取的 蝶形运算流图,按时间抽取的FFT算法以及按频率抽取的FFT算法是DFT的两种快速算法，这些算法同样可以用於IDFT，我们称为快速傅立叶反变换(IFFT)。我
14、们把IDFT公式与DFT公式作以比较：,(1) DFT中的系数 (2) IDFT中要乘以常数,4.4 IDFT的快速算法IFFT,另外一种IFFT算法的实现,以上讨论的时间抽取或频率抽取的FFT运算均可直接用于IFFT运算。但需要将蝶形系数 改为 ，每级乘以系数1/2 。,IFFT计算分为三步： 将X(k)取共轭(虚部乘以-1) ; 对 直接作FFT ； 对FFT的结果取共轭并乘以1/N，就得到x(n)。,按上述方法求得的IFFT，完全不需要改动FFT程序，FFT和 IFFT可以共用一个程序。,总结FFT和IFFT蝶形图的表示方法： 按时间抽取的FFT 变到IFFT时，对应的流图应该是按频率抽
15、取的IFFT。因为：,按时间抽取的FFT输入变量是x(n) x(n)按奇、偶分解,按频率抽取的IFFT输入变量是X(k) X(k)按奇、偶分解,同理，按频率抽取的FFT 变到IFFT时，对应的流图应该是按时间抽取的IFFT。,（1）输入按奇、偶分解：,（2）输入按前、后分解：,4.5 任意基数的FFT算法,前面讨论的FFT算法，序列的长度均满足N=2M ，这种情况实际使用得最多，因为它程序简单，效率 高，使用方便。如果不满足N=2M, 处理方法有两种: 补零: 将x(n)补零,使N增加到临近的一个值， 满足N=2M 。例如：N=30, 补上x(30)=x(31)=0 两个零点, 使得满足N=2
16、5=32, 这样可直接调用以2为基数的FFT运算。有限长序列补零后并不影响其频谱 ， 补零只是频谱的采样点数增加了。, 如 ，如果要求N点DFT值,几乎要补一半0值(2044个0)，显然用第一种方法处理很不经济，因此，如果要求求准确的N点DFT值,可采用任意基数的FFT算法,但其运算效率低于以2为基数的FFT算法。 任意基数的FFT算法的基本思想首先要求N可以分解成因子乘积的形式，N不能是素数。任意基数的FFT算法的基本思想仍是把一个大N点的FFT尽量分解为小点数的FFT。首先将N分解为两个整数p与q的乘积,即：，其步骤为：,(1)将序列x(n)分成p组,p组 r=0,1,q-1,(2)将N点
17、DFT分成p组q点的DFT来运算,例：N6 p3 q2,P组,4.6 Chirp-z变换（简称CZT),采用FFT可以算出全部N点DFT值，而DFT值是Z变换X(z)在Z平面单位圆上的等间隔取样值，但有些情况下不需要求出全部的DFT值。例如： 对于窄带信号，只需要对信号所在的一段频带进 行谱分析，这时,希望频谱的采样集中在这一段频带内，以获得较高的分辨率，而频带以外部分可不考虑。,语音信号处理中，需要知道Z变换的极点所在的复频率位置，即共振峰的位置。如果极点位置离单位圆较远，则其单位圆上的频谱会很平滑，如果采样不是沿单位圆而是沿一条接近这些极点的弧线进行，则在极点所在的频率上将出现明显的尖峰，
18、这些尖峰就是共振峰，它可较准确地测定极点对应的频率。,当N是素数时，不能采用任意基数的FFT算法。鉴于以上三种情况，我们提出一种最有效的变换螺线采样，或称为Chirp-z变换。螺线采样是FFT的另外一种快速算法，它是沿着一条螺旋线的采样。,一.算法原理：,首先对Z平面单位圆内的一段螺线作M等分采样。令采样点 zk=AW-k ，k=0,M-1，M为采样点数，,A0 表示起始取样点的 半径长度，通常A01,0表示起始取样点的相角 0表示两相邻点之间的等分角,A、W为任意复数，其中：,W0 表示螺旋线的伸展率， 则随K的增大螺线外伸， 则螺线内缩(反时针)，W0=1表示半径为A0的一段圆弧,若A0=
19、1表示这段圆弧是单位圆的一部分。把,代入zk=AW-k ，计算Z变换在采样点zk的值：,为了能够利用FFT运算，想办法把Z变换转换为线性卷积形式运算，从而可采用FFT进行，这样可大大提高计算速度。 Z变换在采样点zk上的值为：,zk=AW-k,令：,则,上式说明，如对信号x(n)先进行一次加权处理，加权系数为 ，然后通过一个单位脉冲响应为h(n)的线性系统，最后对该系统的前M点输出再作一次 加权，就可得到全部M点的螺旋线采样值。,系统的单位脉冲响应 是具有二次相位的复指数信号，h(n)可以想象为频率随时间n成线性增长的复指数序列，它和雷达系统中的线性调频信号相似，因此称为Chirp-z变换。,
20、二.CZT的算法实现,输入信号 是有限长序列，长度为N，,但 是无限长序列,而计算0M-1点卷积中h(n)取值只需要n=-(N-1)M-1那一部分的值。因此，可以认为h(n)是一个有限长序列,其长度为L=N+M-1。所以，把Chirp-z变换可以看成两个有限长序列的线性卷积。可利用圆圈卷积通过FFT来实现。h(n)的主值序列 可由h(n)作周期延拓后取0nL-1部分值获得，将 与g(n)作圆周卷积后，其输出的前M个值就是Chirp-z变换的M个值。这个圆周卷积的过程可在频域用FFT实现。,Chirp-z变换的计算步骤：,(1)求h(n)的主值序列(2)求 的L点的FFT: H(k)=FFT (
21、3) 对x(n)加权并补零,(4) G(k)= FFTg(n) , L点 (5) Y(k)= G(k)H(k), L点 (6) y(n)= IFFTY(k) , L点(7) , 0kM-1,利用FFT计算Chirp-z变换,Chirp-z变换的特点：,1）输入序列长度N与输出序列长度M不需要相等； 2）N及M不必是高度复合数，二者均可为素数； 3）相邻采样点zk之间的角间隔0是任意的，即频 率分辨率是任意的； 4）围线是任意的，不必是Z平面上的圆； 5）起始点z0可任意选定，即可从任意频率上开始对输入数据进行窄带高分辨率分析； 6）若 A=1, M=N , ,可用Chirp-z变换计算DFT（
22、即使N为素数）。,4.7 线性卷积的FFT算法线性卷积是求离散系统响应的主要方法之一,许多重要应用都建立在这一理论基础上。第二章中已讨论过，线性卷积可以用圆周卷积的方法替代，为了不产生混叠，序列长度将:长度为N2的序列x(n)延长到L ,补L-N2个零，长度为N1的序列h(n)延长到L ,补L-N1个零，如果LN1+N2-1,则圆周卷积与线性卷积相等,此时,可用FFT运算替代线性卷积，方法如下:,a.计算X(k)=FFTx(n) b. 求H(k)=FFTh(n) c. 求Y(k)=H(k)X(k) k=0L-1 d. 求y(n)=IFFTY(k) n=0L-1可见，只要进行二次FFT,一次IF
23、FT就可完成 线性卷积计算。计算表明, L32时, 上述计算线 性卷积的方法比直接计算线卷积有明显的优越性。 因此,也称圆周卷积的方法为快速卷积法。,上述方法适用于x(n)、h(n)两序列长度比较接近或相等的情况。如果x(n)、h(n)长度相差较多,例如：h(n)为某滤波器的单位脉冲响应,长度有限,用来处理一个很长很长的输入信号x(n),按上述方法,h(n)要补许多零值后再进行计算,这降低了运算速度，发挥不出圆周卷积的优点。为了保持快速卷积法的优越性,可将x(n)分为许多段,每段的长度与h(n)接近, 处理方法有两种： (1)重叠相加法 (2)重叠保留法,(1)重叠相加法由分段卷积的各段相加构
24、 成总的卷积输出,h(n),x(n),序列长度为,序列长度为,假定 表示x(n)序列的第i段 ：则输入序列可表为：于是输出可分解为：其中,d.重叠部分相加构成最后的输出序列。,重叠相加法计算步骤： a. 先对h(n)及 补零，补到具有N点长度， N=N1+N2-1，并且满足N=2M 。 b. 用N点FFT计算 c. 用N点FFT计算,由于 的长度为N，而 的长度为N2，因此相邻两段序列 必然有N-N2=N1-1 点发生重叠，最后的输出应该是这些重叠部分相加起来，再和不重叠部分共同组成输出序列 。重叠相加法对处理一个无限长的语音信号或地震信号都是十分有效的。,有N1-1个点发生重叠,(2)重叠保
25、留法,这种方法和第一种方法稍有不同，即将上面分段序列中补零的部分不是补零，而是保留原来的输入序列值，如果利用FFT实现h(n)和xi(n)的圆周卷积，则每段卷积结果中有N1-1个点不等于线性卷积值需舍去。重叠保留法与重叠相加法的计算量差不多，但省去了重叠相加法最后的相加运算。,h(n),x(n),序列长度为,序列长度为,FFT应用中的几个问题,实数序列的FFT以上讨论的FFT算法都是复数运算,包括序列x(n)也认为是复数,但大多数场合,信号是实数序列,任何实数都可看成虚部为零的复数。例如,求某实信号x(n)的复谱,可认为是实信号加上数值为零的虚部组成复信号(x(n)+j0),再用FFT求其离散
26、傅里叶变换。这种作法很不经济,因为把实序列变成复序列,存储器要增加一倍,且计算机运行时,即使虚部为零,也要进行涉及虚部的运算,造成运算速度下降。合理的解决方法是利用复数据FFT对实数据进行有效计算,下面介绍两种方法。,（1）用一个N点FFT同时计算两个N点实序列的DFT设x(n)、y(n)是彼此独立的两个N点实序列,且X(k)=DFTx(n)， Y(k)=DFTy(n)则X(k)、Y(k)可通过一次FFT运算同时获得。 首先将x(n)、y(n)分别当作一复序列的实部及 虚部,令：g(n)=x(n)+jy(n)通过FFT运算可获得g(n)的DFT值,利用离散傅里叶变换的共轭对称性通过g(n)的F
27、FT运算结果G(k),由上式可得到 X(k)的值。Y(k)的值也可以通过g(n)的FFT运算结果G(k)得到。,作一次点复序列的FFT，再通过加、减法运算就可以将X(k)与Y(k)分离出来。显然，这将使运算效率提高一倍。,（2）用一个N点的FFT获得一个2N点实序列的DFT设x(n)是2N点的实序列,现人为地将x(n)分为偶数组x1(n)和奇数组x2(n)x1(n)=x(2n) n=0,1,N-1x2(n)=x(2n+1) n=0,1,N-1然后将x1(n)及x2(n)组成一个复序列：y(n)=x1(n)+jx2(n)通过N点FFT运算可得到： Y(k)=X1(k)+jX2(k),为求2N点x(n)所对应X(k)，需求出X(k)与 X1(k)、X2(k)的关系 ：,所以 1）由x1(n)及x2(n)组成复序列，经FFT运算求得 Y(k); 2）利用共轭对称性求出 X1(k)、X2(k); 3）最后利用上式求出X(k), 达到用一个N点的FFT计算一个2N点的实序列的DFT的目的。,X(k)=X1(k)+W2Nk X2(k),
