数字信号处理 第4章 快速傅里叶变换(FFT)-

第4章快速傅里叶变换(FFT)4.1引言4.2基2FFT算法问题的提出解决问题的思路与方法基2时间抽取FFT算法基2频率抽取FFT算法4.3进一步减少运算量的措施

4.1引言

DFT是信号分析与处理中的一种重要变换。因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换(简称FFT)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。直到1965年发现了DFT的一种快速算法以后，情况才发生了根本的变化。4.2基2FFT算法

问题的提出4点序列{2，3，3，2}DFT的计算复杂度复数加法N(N-1)复数乘法N

2如何提高DFT的运算效率?时间抽取FFT解决问题的思路1.将长序列DFT分解为短序列的DFT2.利用旋转因子的周期性、对称性、可约性。时间抽取FFT旋转因子的性质(1)周期性(2)对称性(3)可约性时间抽取FFT解决问题的方法将时域序列逐次分解为一组子序列，利用旋转因子的特性，由子序列的DFT来实现整个序列的DFT。基2时间抽取(Decimationintime)FFT算法基2频率抽取(Decimationinfrequency)FFT算法时间抽取FFT基2时间抽取FFT算法基2时间抽取FFT算法推导基2时间抽取FFT算法流图基2时间抽取FFT算法的计算复杂度基2时间抽取FFT算法流图规律时间抽取FFT基2时间抽取FFT算法推导时间抽取FFT基2时间抽取FFT算法推导因此有:时间抽取FFT基2时间抽取FFT算法流图N=2x[k]={x[0],x[1]}时间抽取FFT4点基2时间抽取FFT算法流图x[0]x[2]x[1]x[3]X1[0]X1[1]X2[0]X2[1]2点DFT2点DFT-1-1-1-1X

[0]X

[1]X

[2]X

[3]时间抽取FFT4点基2时间抽取FFT算法流图时间抽取FFT8点基2时间抽取FFT算法流图4点DFT4点DFTx[0]x[2]x[4]x[6]x[1]x[3]x[5]x[7]X1[0]X1[1]X1[2]X1[3]X2[0]X2[1]X2[2]X2[3]X

[0]X

[1]X

[2]X

[3]X

[4]X

[5]X

[6]X

[7]-1-1-1-1时间抽取FFT4点DFT4点DFTx[0]x[2]x[4]x[6]x[1]x[3]x[5]x[7]X1[0]X1[1]X1[2]X1[3]X2[0]X2[1]X2[2]X2[3]X

[0]X

[1]X

[2]X

[3]X

[4]X

[5]X

[6]X

[7]-1-1-1-18点基2时间抽取FFT算法流图时间抽取FFT第一级第二级第三级8点基2时间抽取FFT算法流图时间抽取FFT算法的计算复杂度复乘次数复乘次数NN2时间抽取FFT基2时间抽取FFT算法流图第一级第二级第三级时间抽取FFTFFT算法流图旋转因子规律第二级的蝶形系数为，蝶形节点的距离为2。第一级的蝶形系数均为，蝶形节点的距离为1。第三级的蝶形系数为，蝶形节点的距离为4。第M级的蝶形系数为，蝶形节点的距离为N/2。时间抽取FFT倒序k0k1k2x[k2k1k0]x[000]x[100]x[010]01011]12x[kk0]x[k2k101x[110]x[001]x[101]x[011]x[111]01010101时间抽取FFT例：试利用N=4基2时间抽取的FFT流图计算8点序列x[k]={1,-1,1,-1,2,1,1,2}的DFT。解：根据基2时间抽取FFT算法原理，8点序列的DFTX[m]可由两个4点序列的DFTX1[m]和X2[m]表达。如果按照序列x[k]序号的奇偶分解为x1[k]和x2[k]，则存在其中x1[k]={1,1,2,1}，x2[k]={-1,-1,1,2}，X1[m]和X2[m]可通过4点的FFT来计算。时间抽取FFT例：试利用N=4基2时间抽取的FFT流图计算8点序列x[k]={1,-1,1,-1,2,1,1,2}的DFT。解：

X1[m]={5,-1,1,-1},X2[m]={1,-2+3j,1,-2-3j}利用上述公式，可得序列x[k]的DFTX[m]为X[m]={6,-0.293+3.535j,1+j,-1.707+3.535j,4,-1.707-3.535j,1-j,-0.293-3.535j}时间抽取FFT则x(n)的DFT为由于所以序列x(n)的长度为N,M为自然数按n的奇偶把x(n)分解为两个N/2点的子序列N/2点DFTX(K)前半部分X(K)后半部分基2时间抽取FFT算法流图N=2x[k]={x[0],x[1]}图4.2.2N点DFT的一次时域抽取分解图(N=8)与第一次分解相同，将x1(r)按奇偶分解成两个N/4长的子序列x3(l)和x4(l)，即那么，X1(k)又可表示为(4.2.9)式中同理，由X3(k)和X4(k)的周期性和WmN/2的对称性Wk+N/4

N/2=-WkN/2最后得到：(4.2.10)用同样的方法可计算出(4.2.11)其中图4.2.3N点DFT的第二次时域抽取分解图(N=8)图4.2.4N点DIT―FFT运算流图(N=8)4.2.3DIT―FFT算法与直接计算DFT运算量的比较每一级运算都需要N/2次复数乘和N次复数加(每个蝶形需要两次复数加法)。所以，M级运算总共需要的复数乘次数为复数加次数为例如，N=210=1024时图4.2.5FFT算法与直接计算DFT所需乘法次数的比较曲线4.2.4DIT―FFT的运算规律及编程思想FFT的每级（列）计算都是由N个复数数据（输入）两两构成一个蝶型（共N/2个蝶形）运算而得到另外N个复数数据（输出）。当数据输入到存储器以后，每一组运算的结果，仍然存放在这同一组存储器中直到最后输出。例：将x(0)放在单元A(0)中，将x(4)放在单元A(1)中，W80

放在一个暂存器中。将x(0)+W80x(4)→送回A(0)单元将x(0)-W80x(4)→送回A(1)单元X3(0)X3(1)x(0)x(4)1）

原位运算(亦称同址计算)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)回顾：N点DIT―FFT运算流图(N=8)输入数据、中间运算结果和最后输出均用同一存储器。如上所述，N点DIT―FFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子WNP，称其为旋转因子，p称为旋转因子的指数。2)旋转因子的变化规律观察FFT运算流图发现，第L级共有2L-1个不同的旋转因子。N=23=8时的各级旋转因子表示如下：L=1时，WNp=WN/4J，N/4=21=2L，J=0L=2时，WNp=WN/2J，N/2=22=2L，J=0，1L=3时，WNp=WNJ，N=23=2L，J=0，1，2，3对N=2M的一般情况，第L级的旋转因子为：设序列x(n)经时域抽选(倒序)后，存入数组X中。如果蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶形运算可表示成如下形式：下标L表示第L级运算，XL(J)则表示第L级运算后数组元素X(J)的值。蝶形运算两节点的距离B=2L-1,L=1…,M例如N=8=23,第一级(列)距离为21-1=1,第二级(列)距离为22-1=2，

第三级(列)距离为23-1=4。3)编程思想及流程图开始送入x(n)和N=2M调整输入x(n)的顺序for(L=1;L<=M;L++)B=2L-1for(J=0;J<=B-1;J++)p=J·2M-Lfor(k=J;k<=N-1;k=k+2L)输出结果结束4）码位倒序由N=8蝶形图看出：原位计算时，FFT输出的X(k)的次序正好是顺序排列的，即X(0)…X(7),但输入x(n)都不能按自然顺序存入到存储单元中，而是按x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)的顺序存入存储单元，即为乱序输入，顺序输出。这种顺序看起来相当杂乱，然而它是有规律的。即码位倒读规则。n=00n=10n=01n=11n=01n=1101010101

(n2)(000)0(100)4(010)2(110)6(001)1(101)5(011)3(111)7（偶）（奇）码位倒读规则由奇偶分组造成的,以N=8为例说明如下:输入序列的序号n二进制为(n2n1n0)码位倒序二进制为(n0n1n2)

00

0

00

0

00

10

0

11

0

04

20

1

00

1

02

30

1

11

1

06

41

0

00

0

11

51

0

11

0

15

61

1

00

1

13

71

1

11

1

17自然顺序n二进制nnn

倒位序二进制nnn倒位顺序n^210012自然顺序n二进制码表示码位倒读码位倒置顺序n’以N=8为例：0123456700000101001110010111011100010001011000110101111104261537看出：码位倒读后的顺序刚好是数据送入计算机内的顺序。A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(8)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)变址处理方法存储单元自然顺序变址倒位序对于数N，在其二进制最高位加1，等于加N/2。若已知某个反序号为J，为求下一个反序号，可先判J的最高位：

1)若为0，则把该位变成1（即加N/2）就得到下一个反序号，2)若为1，则需判断次高位：

①若次高位为0，则把最高位变0（相当减去N/2）后，再把次高位变1（即加N/4）。

②若次高位为1，则需判断次次高位……分析：倒序排列算法的流程图正序序列已在数组A[]中，输入NLH=N/2,J=LH,N1=N-2J=J-kk=k/2k=LHJ<kJ=J+kT=A(I)A(I)=A(J)A(J)=Tfor(i=1;i<=N1;i++)i≥JNYYN第四节按频率抽选的基2-FFT算法在基2快速算法中，频域抽取法FFT也是一种常用的快速算法，简称DIF―FFT。设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式DIF―FFT一次分解运算流图(N=8)4点DFT4点DFTx(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)x1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)DIF―FFT二次分解运算流图(N=8)DIF―FFT运算流图(N=8)时间抽取算法与频率抽取算法的比较1)频率抽选法和时间抽选法总的计算量是相同的复乘：复加：2)频率抽取法和时间抽取法一样，都适用于原位运算，即蝶形的输入和输出占用同一个存储单元。3)均存在码位倒序问题。4)频率抽选法和时间抽选法一样，基本运算也是蝶形运算。但两者的蝶形形式略有不同。第五节IDFT的快速算法－IFFT上述FFT算法流图也可以用于离散傅里叶逆变换(InverseDiscreteFourierTransform，简称IDFT)。比较DFT和IDFT的运算公式：1)旋转因子：2)系数：DIT―IFFT运算流图DIT―IFFT运算流图(防止溢出)如果希望直接调用FFT子程序计算IFFT，则可用下面的方法：对上式两边同时取共轭，得：例1、如果通用计算机的速度为平均每次复乘需要5s，每次复加需要0.5s，用它来计算512点的DFT[x(n)]，问：1）直接计算需要多少时间？2）用FFT需要多少时间？解：1）用DFT进行运算：复乘：T1=N2×5×10-6=1.31072秒复加：T2=N(N-1)×0.5×10-6=0.130816秒总共：T=T1+T2=1.441536秒2）用FFT进行运算：复乘：T1’=(N/2)log2N×5×