数字信号处理-第4章

第四章离散傅里叶变换快速算法

FastFourierTransform

(FFT)1

快速傅里叶变换（FFT）是计算DFT的一种快速有效方法，

并不是新的傅立叶变换式。

从前面的讨论中看到，有限长序列在数字技术中占有很重要的地位。有限长序列的一个重要特点是其频域也可以离散化，即离散傅里叶变换（DFT）。2DFT是信号分析与处理中的一种重要变换。但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。1965年库里和图基发现了DFT的一种快速算法，使DFT的运算效率提高1－2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了条件，推动了数字处理技术的发展。1984年，提出了分裂基快速算法，使运算效率进一步提高；3一、直接计算DFT的问题及改进的方法DFT的定义两者形式类似，差别只在于的指数符号不同，及常数因子运算量是相同的41、运算量(a+jb)(c+jd)=(ac-bd)+j(ad+bc)∑

-=10)(NnnkNWnx562、减少运算量的途径之一具有如下特性：使DFT运算中的有些项可以合并。对称性周期性可约性knNNkNnNWW=--)()(=kNNkNWW-=+)2/(NNW-=2/,17减少运算量的途径之二可将长序列的DFT分解为短序列的DFT。一个长序列DFT乘加运算比多个短序列DFT乘加运算多长序列DFT可由多个短序列DFT组合而成8二、按时间抽取(DIT)的基-2FFT算法

设M为自然数将长度为N的序列x(n)按n的奇偶分成两组)12/(,1,0),12()()12/(,1,0),2()(21-=+=-==NrrxrxNrrxrx……1、算法原理9则x(n)的DFT为)()(21kXWkXkN+=∑∑

--+=12/2/212/2/1)()(NkrNkNNkrNWrxWWrxr=0r=01、算法原理10其中是x(2r)与x(2r+1)的N/2点DFT。1、算法原理11由于得12,,1,0-=Nk…1、算法原理12信号流图

将X1

(k)和X2

(k)合成X(k)运算可归结为：蝶形运算的简化Wa+bWa-bW-W-1a+bWa-bWWabab(a)(b)X1(k)X2(k)X1(k)+WNkX2(k)X1(k)WNkX2(k)WNk图：蝶形运算符号134点基2时间抽取FFT算法流图x[0]x[2]x[1]x[3]X1[0]X1[1]X2[0]X2[1]2点DFT2点DFTX

[0]X

[1]X

[2]X

[3]14最后剩下的是2点DFT，它可以用一个蝶形结表示：x(0)x(1)X(0)X(1)158点基2FFT算法N/2点DFTN/2点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)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)｛N=8点的DIT-2FFT(时域抽取图)一次分解图16

(3)第二次分解：

将x1(r)按r取奇、偶可分解成2个长度为N/4的子序列

x3(l)=x1(2l)、

x4(l)=x1(2l+1)，根据上面推导可得：X1(k)=X3(k)+WN/2kX4(k)，k=0,1,…,N/2-1将x2(r)按r取奇、偶可分解成2个长N/4的子序列

x5(l)=x2(2l)，l=0,1,…，N/4x6(l)=x2(2l+1)；同理得l=0,1,…，N/4-1；X1(k+N/2)=X3(k)WN/2kX4(k)，k=0,1,…,N/4-1；X1(k)=X3(k)+WN/2kX4(k)，k=0,1,…,N/4-1；X2(k)=X5(k)+WN/2kX6(k)，k=0,1,…,N/4-1；X2(k+N/2)=X5(k)WN/2kX6(k),k=0,1,…,N/4-1；17N/4点DFTx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)N/4点DFTN/4点DFTN/4点DFTX3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN/20WN/21WN/21WN/20N=8点DFT的二次时域抽取分解图X1(k+N/2)=X3(k)WN/2kX4(k)X1(k)=X3(k)+WN/2kX4(k)X2(k)=X5(k)+WN/2kX6(k)X2(k+N/2)=X5(k)WN/2kX6(k)k=0,1,…,N/4-1；18再次分解，对N=8点，可分解三次。X(5)N=8点DIT-FFT运算流图x(0)x(4)x(2)x(6)x(1)x(5)x(3)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(6)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN/20WN/21WN/20WN/40WN/40WN/40x(7)WN/21WN/40L=1级L=2L=3X(7)X3(0)X1(0)192点DFT的蝶形结表示：x(0)x(1)X(0)X(1)208点基2DIT-FFT算法N=8点DIT-FFT运算流图x(0)x(4)x(2)x(6)x(1)x(5)x(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)WN0WN2WN0WN0WN0WN0x(7)WN2WN0L=1级L=2L=3X(7)21基2

DIT-FFT算法

由于这种方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的，所以称为“按时间抽取法”或“时间抽取法”。22

DIT―FFT算法与直接计算DFT运算量的比较1、直接DFT运算N点运算：复数乘次数：N×N

复数加次数：N×(N-1)2、用DIT-FFT作N点运算：复数乘次数：M×N/2=N/2×log2N；复加次数：2×N/2×M=N×log2N。整个运算流图中有M级蝶形，每一级运算流图中有N/2个蝶形，每个蝶形需一次复乘和两次复数加运算。2324时间抽取法FFT的运算特点（1）蝶形运算（2）原位计算（3）序数重排（4）蝶形类型随迭代次数成倍增加25（1）蝶形运算对于N=2M，总是可以通过M次分解最后成为2点的DFT运算。这样构成从x(n)到X(k)的M级运算过程。从上面的流图可看到，每一级运算都由N/2个蝶形运算构成。因此每一级运算都需要N/2次复乘和N次复加，这样经过时间抽取后M级运算总共需要的运算：复乘复加

而直接运算时则与N2成正比。26（2）原位计算当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器，这叫原位计算。每一级运算均可在原位进行，这种原位运算结构可节省存储单元，降低设备成本，还可节省寻址的时间。27（3）序数重排对按时间抽取FFT的原位运算结构，当运算完毕时，正好顺序存放着x(0)，x(1)，x(2)，…x(7)，因此可直接按顺序输出，但这种原位运算的输入x（n）却不能按这种自然顺序存入存储单元中，而是按x(0)，x(4)，x(2)，x(6)，…，x(7)的顺序存入存储单元，这种顺序看起来相当杂乱，然而它也是有规律的。当用二进制表示这个顺序时，它正好是“码位倒置”的顺序。例如，原来的自然顺序应是x(1)的地方，现在放着x(4)，用二进制码表示这一规律时，则是在

x（001）处放着x（100），

x（011）处放着x（110）。28自然顺序二进码表示码位倒置码位倒序0000000010011004201001023011110641000011510110156110010371111117

在实际运算中，一般直接将输入数据x（n）按码位倒置的顺序排好输入很不方便，总是先按自然顺序输入存储单元，然后再通过变址运算将自然顺序的存储转换成码位倒置顺序的存储，然后进行FFT的原位计算。目前有许多通用DSP芯片支持这种码位倒置的寻址功能。29造成倒位序的原因是输入序列x(n)，按标号(序号)n的奇偶不断地分组造成的。倒位序的形成30（4）蝶形类型随迭代次数成倍增加观察8点FFT的三次迭代运算:第一级迭代，有一种类型的蝶形运算系数W80，两个数据点间隔为1第二级迭代，有二种类型的蝶形运算系数W80

、W82

，参加运算的两个数据点间隔为2。第三级迭代，有四类蝶形运算系数W80

、W81

、W82

、W83

，参加运算的两个数据点间隔为4。结论：每迭代一次，蝶形类型增加一倍，数据点间隔也增大一倍。每一级的取数间隔和蝶形类型种类均为2i-1，i=1,2,…M。31N点DIT—FFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子，称其为旋转因子，p称为旋转因子的指数。但各级的旋转因子和循环方式都有所不同，为了编写计算程序，应先找出旋转因子与运算级数的关系。用L表示从左到右的运算级数（L=1，2，、、，M）第L级共有个不同的旋转因子。N=8时各级旋转因子表示如下：32对的一般情况，第L级的旋转因子为从而可以确定第L级运算的旋转因子。33FFT算法流图3435倒位序的实现（a）只要将顺序二进制数（n2n1n0）的二进制位倒置，得（n0n1n2）。根据这种规律，容易用硬件电路和汇编语言产生倒位序数。（b）用高级语言程序实现时，必须找出产生倒序数的十进制运算规律。36000000001

001

4

100201020103

011

6

1104

100

1

001510151016

110

3

01171117111左边为按自然顺序排列的二进制数，下面的一个数是上面一个数的最低位上加上1，且向高位进位。右边为倒位序数，下面的一个数是上面一个数的最高位上加上1，且由高位向低位进位。称为反向进位加法顺序与倒序二进制对照表可由当前任一倒序值求得下一个倒序值37

N=2M，用M位二进制数表示，则从左至右的十进制权值为：

N/2、N/4、N/8,…、2,1

对倒序数J，其下一个序数是在该序数J的二进制首位码加1，相当于十进制运算J+N/2。计算机上倒序数的实现过程为：J>N/2？J>N/4输入当前倒序数十进制数值JNYNYJ>N/2MNY结束J=J-N/2倒序数的十进制运算规律38N=8倒位序的变址处理分析上图N=8点倒序规律，顺序数I与倒序数J

存在关系：

I=J时，不交换；

I<J时，交换存储器内容。

I>J时，不交换，直接更新序数值；39LH=N/2J=LHN1=N-2I=1,N1I≥JT=A(I)A(I)=A(J)A(J)=TK=LHJ<KJ=J+KJ=J-KK=K/2NN

YY倒序程序框图计算倒序值的程序流图4041三、按频率抽取（DIF）的基-2FFT算法1、算法原理仍设序列点数为N=2M，M为正整数。将X(k)按k的奇偶分组之前，先将输入序列按前一半、后一半分开。k=0,1,…,N-142

k=0,1,…,N-1按k的奇偶可将X(k)分为两部分:12,,1,0-=Nr…4312,,1,0-=Nr…令则x(n)x1(n)=x(n)+x(n+N/2)x2(n)=[x(n)x(n+N/2)]WNnWNnx(n+N/2)DIF-FFT蝶形运算流图符号44X(3)N/2点DFTN/2点DFTx(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)WN0WN1WN2WN3X(0)X(2)X(4)X(6)X(1)X(5)X(7)N=8时第1次分解的运算流图45N=2M，N/2仍是偶数，继续将N/2点进行分解。将输入序列x1(n)、x2(n)分别按前、后对半分解成4个长N/4的子序列,其n=0,1,…，N/4-1｛x3(n)=x1(n)+x1(n+N/4)

x4(n)=[x1(n)-x1(n+N/4)]WN/2n

x5(n)=x2(n)+x2(n+N/4)

x6(n)=[x2(n)-x2(n+N/4)]WN/2n｛N/4点DFTWN0WN1WN2WN3x(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)WN0WN2WN0WN2N/4点DFTN/4点DFTN/4点DFTx1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)x3(0)x3(1)x4(0)x4(1)x5(0)x5(1)x6(0)x6(1)DIF―FFT二次分解运算流图(N=8)

46经过三次分解后，DIF―FFT运算流图(N=8)为：x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)WN0WN1WN2WN3X(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)x3(0)x3(1)x4(0)x4(1)x5(0)x5(1)x6(0)x6(1)WN0WN2WN0WN2WN0WN0WN0WN0472、DIF-FFT运算规律

DIF-FFT算法也可采用原位计算；N=2M时，共有M级运算，每级共有N/2个蝶形，DIT与DIF算法的运算次数相同。DIT与DIF不同的是：

DIF-FFT算法输入序列为自然序列，输出为倒序序列。因此，在M级运算完成后，需对输出数据进行倒序才能得到自然顺序的X(k)。

蝶形运算符号不同：DIT-FFT蝶形是先相乘，后加/减；而DIF-FFT蝶形是先加/减，后相乘。DIT-FFT蝶形运算符号X1(k)X2(k)X1(k)+WNkX2(k)X1(k)WNkX2(k)WNkx2(n)=[x(n)x(n+N/2)]WNnx(n)x1(n)=x(n)+x(n+N/2)WNnx(n+N/2)DIF-FFT蝶形运算流图符号48FFT的变形运算流图(一)49FFT的变形运算流图(二)50FFT减少运算量的方法旋转因子的几个特殊点取值存储旋转因子表提高运算速度51五、FFT应用中的几个问题1)IDFT的运算方法

以以上所讨论的FFT算法可用于IDFT运算——简称为IFFT比比较IDFT的定义式：IDFT：

DFT：X(k)=DFT[x(n)]=52IDFT与DFT的差别：1）把DFT中的每一个系数改为，2）再乘以常数1/N，则以上所讨论的时间抽取或频率抽取的FFT运算均可直接用于IDFT运算，当然，蝶形中的系数应改为。53DIT―IFFT运算流图实际中，为了避免发生溢出，将1/N分配到每一级蝶形运算中，因为，所以每一级的每个蝶形输出到乘以1/2.

54

b）第二种方法，完全不需要改动FFT程序，而是直接利用它作IFFT。考虑到故

IFFT计算分三步：①将X（k）取共轭（虚部乘以-1）

②对直接作FFT③对FFT的结果取共轭并乘以1/N，得x（n）。虽然2次用了取共轭运算，但可以和FFT共用一个子程序，实现方便。5556求某实信号x(n)的复谱,可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用FFT求其离散傅里叶变换。把实序列变成复序列,存储器要增加一倍,且计算机运行时,即使虚部为零,也要进行涉及虚部的运算。合理的解决方法是利用复数据FFT对实数据进行有效计算。2）实数序列的FFT

57

（1）用

一个N点FFT同时计算两个N点实序列的DFT

设x

(n)、y

(n)是彼此独立的两个N点实序列,且X

(k)=DFT[x

(n)]，Y

(k)=DFT[y(n)]则X

(k)、Y(k)可通过一次FFT运算同时获得。首先将x

(n)、y(n)分别当作一复序列的实部及虚部,令g(n)=x

(n)+jy(n)对g(n)进行N点FFT运算，得到则：58

(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)=