数字信号处理4 fft

第四章快速傅里叶变换（FFT)4.2、直接计算DFT的问题和改进途径4.3、按时间抽取的FFT算法(DIT—FFT)4.4、按频率抽取的FFT算法(DIF—FFT)4.5、N为复合数的FFT算法—统一FFT算法4.6、分裂基FFT算法4.7、FFT应用4.8、线性调频z变换（CZT)算法频谱分析和DFT运算很重要，但在很长一段时间里，由于DFT运算复杂，并没有得到真正的运用。频谱分析以前大多采用模拟信号滤波的方法解决，直到1965年首次提出DFT运算的一种快速算法以后，情况才发生了根本变化。4.1引言2有限长序列在数字技术中占有很重要的地位。有限长序列的一个重要特点是其频域也可以离散化表示，即离散傅里叶变换（DFT）。

认识到DFT运算的一些内在规律，发展和完善了高速有效的运算方法——快速傅里叶变换（FFT）算法。FFT的出现，使DFT的运算大大简化，运算时间缩短1～2个数量级，使DFT的运算在实际中得到广泛应用。4.1引言31

DFT运算的特点

分析有限长序列x(n)进行一次DFT运算所需的运算量。

一般，和

都是复数，每计算一个值，要进行N次复数相乘，N-1次复数相加，一共有N个点，完成全部DFT运算，需要N2次复数相乘和N(N-1)次复数相加。44.2

直接计算DFT的问题和改进途径乘法比加法复杂，需要的运算时间多，尤其是复数相乘；每个复数相乘包括4个实数相乘和2个实数相加：1

DFT运算的特点51

DFT运算的特点整个DFT运算（N点）需要4N2实数相乘和2N（2N-1）次实数相加。

每个复数相加包括2个实数相加；每计算一个要进行4N次实数相乘和2N+2（N-1）=2（2N-1）次实数相加；61

DFT运算的特点

例，计算N=10点的DFT，需要100次复数相乘，而N=1024点时，需要1048576（一百多万）次复数乘法。如果要求实时处理，则要求有很高的计算速度才能完成上述计算量。

在DFT计算中，不论是乘法和加法，运算量均与N2成正比。N较大时，运算量十分可观。7

反变换IDFT与DFT的运算结构相同，只是多乘一个常数1/N，所以二者的计算量相同。1

DFT运算的特点82FFT算法的基本思想考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量：1）系数是一个周期函数，它的周期性和对称性可用来改进运算，提高计算效率。92FFT算法的基本思想利用这些周期性和对称性，使DFT运算中有些项可合并；2）把长度为N点的大点数的DFT运算依次分解为若干个小点数的DFT。因为DFT的计算量正比于N2，N小，计算量也就小。10 FFT算法是基于这样的基本思想发展起来的。有多种形式，基本上可分为两类：

时间抽取法（DIT）

频率抽取法（DIF）

2FFT算法的基本思想114.3按时间抽取的FFT算法:DIT--FFT先从一个特殊情况开始，假定N是2的整数次方：

N=2M，M：正整数--基-2FFT首先将序列x(n)分解为两组，一组为偶数项，一组为奇数项，

12将DFT运算也相应分为两组：131算法原理因为故其中1算法原理

注意：H(k)，G(k)只有N/2个点，即k=0，1，…，N/2-1，还必须应用系数

的周期性和对称性求X(k)的N/2～N-1点：1算法原理15得：一个N点的DFT被分解为两个N/2点的DFT，这两个N/2点的DFT再合成为一个N点DFT。1算法原理16得

依此类推，G(k)和H(k)可以继续拆分下去，这种按时间抽取算法是在输入序列分成越来越小的子序列上执行DFT运算，最后再合成为Ｎ点的DFT。1算法原理17蝶形信号流图：将G(k)和H(k)合成X(k)运算可归结为：蝶形运算流图符号181算法原理

流图需一次乘法、两次加减法。

运算结构像一蝴蝶，称作蝶形运算结构，简称蝶形结；采用这种表示法，可以将FFT过程用流图表示。19蝶形信号流图：1算法原理N=23=8

20蝶形信号流图：1算法原理N=23=8

21蝶形信号流图：1算法原理N=23=8

22蝶形信号流图：1算法原理N=23=8

23蝶形信号流图：1算法原理N=23=8

24蝶形信号流图：1算法原理N=23=8

25蝶形信号流图：1算法原理N=23=8

26蝶形信号流图：1算法原理N=23=8

27蝶形信号流图：1算法原理

按照这个办法，继续把N/2用２除；

N=2M，N/2可以被２整除，可以对两个N/2点的DFT再分别作进一步的分解，即对{G(k)}和{H(k)}的计算；{G(k)}和{H(k)}又可以分别通过计算两个长度为N/4=2点的DFT，进一步节省计算量。

一个８点的DFT就可以分解为四个２点的DFT。28蝶形信号流图：1算法原理29x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)第三级抽取过程：1算法原理301,3,5,70,2,4,6x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)第二级第三级抽取过程：1算法原理1,3,5,70,2,4,60,42,61,53,7x(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)第一级第二级第三级31抽取过程：1算法原理

由四个2点DFT组成8点DFT32抽取过程：1算法原理

由四个2点DFT组成8点DFT34抽取过程：1算法原理

由四个2点DFT组成8点DFT35抽取过程：1算法原理最后剩下的是2点DFT，它可以用一个蝶形结表示：这样，得到一个8点的完整的按时间抽取运算的流图。36抽取过程：1算法原理由于这种方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的，所以称为“按时间抽取法”或“时间抽取法”。37抽取过程：1算法原理N点DIT―FFT运算流图(N=8)38抽取过程：1算法原理2时间抽取法FFT的运算特点1）蝶形运算2）原位计算3）序数重排391）蝶形运算

对于N=2M，总是可以通过M次分解最后成为2点的DFT运算。这样构成从x(n)到X(k)的M级运算过程。从流图可知，每一级运算都由N/2个蝶形运算构成。每一级运算都需要N/2次复乘和N次复加。402时间抽取法FFT的运算特点1）蝶形运算经过时间抽取后M级运算总共需要的运算：复乘复加直接运算时则与N2成正比。N=2048，N2=4194304，(N/2)log2N=11264，N2/[(N/2)log2N]=392.4。FFT显然要比直接法快得多。412时间抽取法FFT的运算特点422时间抽取法FFT的运算特点432时间抽取法FFT的运算特点2）原位计算当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器，这叫原位计算。

每一级运算均可在原位进行，这种原位运算结构可节省存储单元，降低设备成本，还可节省寻址的时间。442时间抽取法FFT的运算特点452时间抽取法FFT的运算特点2）原位计算N点DIT―FFT运算流图(N=8)382）原位计算2时间抽取法FFT的运算特点3）序数重排对按时间抽取FFT的原位运算结构，当运算完毕时，正好顺序存放着X(0)，X(1)，X(2)，…，X(7)，因此可直接按顺序输出。但原位运算的输入x(n)却不是自然顺序，而是按x(0)，x(4)，x(2)，x(6)，…，x(7)的顺序存入存储单元，这种顺序看起来相当杂乱，然而它也是有规律的。462时间抽取法FFT的运算特点3）序数重排当用二进制表示这个顺序时，它正好是“码位倒置”（位翻转）的顺序。

例如，原来的自然顺序应是x(1)的地方，现在放着x(4)，用二进制码表示这一规律时，则是在

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

x（011）处放着x（110）。472时间抽取法FFT的运算特点码位倒置顺序表（位翻转）自然顺序二进码表示码位倒置码位倒置顺序00000000100110042010010230111106410000115101101561100103711111172时间抽取法FFT的运算特点3）序数重排在实际运算中，直接将输入数据x(n)按码位倒置的顺序排好输入很不方便，总是先按自然顺序输入存储单元，然后再通过变址运算将自然顺序的存储转换成码位倒置顺序的存储，然后进行FFT的原位计算。目前有许多通用DSP芯片支持这种码位倒置的寻址功能。492时间抽取法FFT的运算特点倒序规律502时间抽取法FFT的运算特点4.4按频率抽取[DIF]的FFT算法51DITFFT：52将

分成前后两半：前半子序列：后半子序列：4.4按频率抽取[DIF]的FFT算法53由定义：4.4按频率抽取[DIF]的FFT算法按k的奇偶将X(k)分为两部分：54则4.4按频率抽取[DIF]的FFT算法令554.4按频率抽取[DIF]的FFT算法4.4按频率抽取[DIF]的FFT算法则56DIF―FFT一次分解运算流图(N=8)4.4按频率抽取[DIF]的FFT算法57DIF―FFT二次分解运算流图(N=8)4.4按频率抽取[DIF]的FFT算法58DIF―FFT运算流图(N=8)4.4按频率抽取[DIF]的FFT算法594.5N为复合数的FFT算法(统一的FFT算法)处理方法补零N为素数，直接DFT或CZTN为复合数，用统一的FFT算法，基-2是特例6061算法原理4.5N为复合数的FFT算法(统一的FFT算法)614.5N为复合数的FFT算法(统一的FFT算法)n1n001230012314567289101161算法原理4.5N为复合数的FFT算法(统一的FFT算法)61算法效率4.5N为复合数的FFT算法(统一的FFT算法)运算复数乘法复数加法M个L点DFTML2ML(L-1)乘因子NL个M点DFTLM2LM(M-1)N点复合数FFTN(M+L+1)N(M+L-2)直接计算N点DFTN2N(N-1)加速比例N=60=12*53.33.9原理：序列按模4同余数抽取为4个序列，一直分到仅剩4点。71基-4DIT-FFT算法4.5N为复合数的FFT算法(统一的FFT算法)展开：扩展：72基-4DIT-FFT算法4.5N为复合数的FFT算法(统一的FFT算法)蝶形图734.5N为复合数的FFT算法(统一的FFT算法)基-4DIT-FFT算法规律小结输入反序，输出正序原位运算蝶形级数每级蝶形数复乘复加744.5N为复合数的FFT算法(统一的FFT算法)基-4DIT-FFT算法把x(n)按前后分作4段，求DFT。754.5N为复合数的FFT算法(统一的FFT算法)基-4DIF-FFT算法按模4同余数抽取X(k)764.5N为复合数的FFT算法(统一的FFT算法)基-4DIF-FFT算法蝶形图

先不计算DFT，而继续按先后分4段，直到剩4点为止。774.5N为复合数的FFT算法(统一的FFT算法)基-4DIF-FFT算法综合基-2、基-4DIF而形成的一种高效、实用算法。784.6分裂基FFT算法4.6分裂基FFT算法794.6分裂基FFT算法814.6分裂基FFT算法一个N点的DFT计算可分解为一个N/2点DFT和两个N/4点DFT。804.6分裂基FFT算法L-蝶形图824.6分裂基FFT算法N/4个“L”形，每个L形包含2个乘法运算，其余为加法运算。蝶形图中求和向上者按基-2，向下者按基-4(包含L形)进行分解，直到4点为止。834.6分裂基FFT算法8点分裂基FFT流图844.6分裂基FFT算法85规律结构与基-2DIF相同，仅乘系数不同；运算量较基-2减少；每级包含L形个数：4.6分裂基FFT算法每个L形2次复乘，总共基-2复乘数基-4复乘数864.7FFT应用1IDFT的运算方法

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

87IDFT与DFT的差别：把DFT中的每一个系数

改为再乘以常数1/N

因此，以上所讨论的时间抽取或频率抽取的FFT运算均可直接用于IDFT运算，当然，蝶形中的系数应改为。881IDFT的运算方法

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

1IDFT的运算方法

90

IFFT计算分三步：①将X(k)取共轭，即虚部乘以-1

②对直接作FFT③对FFT的结果取共轭并乘以1/N，得x(n)1IDFT的运算方法

91FFT算法都是复数运算,包括序列x(n)也认为是复数;大多数场合,信号是实数序列,任何实数都可看成虚部为零的复数;例如,求某实信号x(n)的复谱,可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用FFT求其离散傅里叶变换。2实数序列的FFT

这种作法很不经济,把实序列变成复序列,存储器要增加一倍,且计算机运行时,即使虚部为零,也要进行涉及虚部的运算,浪费了运算量。

合理的解决方法是利用复数据FFT对实数据进行有效计算,下面介绍两种方法。2实数序列的FFT92用一个N点FFT同时计算两个N点实序列的DFT

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

(k)=DFT[x

(n)]，Y

(k)=DFT[y(n)]

则X

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

2实数序列的FFT932实数序列的FFT

首先，将x

(n)、y(n)分别当作一复序列的实部及虚部,令94对g(n)进行N点FFT运算，得到则：2实数序列的FFT952实数序列的FFT用一个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)96通过N点FFT运算可得到：

Y(k)=X1(k)+jX2(k) 根据前面的讨论，得到

972实数序列的FFT

为求2N点x(n)所对应X(k)，需求出X(k)与X1(k)、X2(k)的关系：

2实数序列的FFT98所以2实数序列的FFT99计算步骤：1）将序列x(n)按奇偶序拆分成x1(n)及x2(n)，并组成复序列y(n)，经FFT运算求得Y(k)，2）利用共轭对称性求出X1(k)、X2(k)，3）最后利用上式求出X(k),实现用一个N点的FFT计算一个2N点的实序列的DFT。2实数序列的FFT100

线性卷积是求离散系统响应的主要方法之一,许多重要应用都建立在这一理论基础上,如卷积滤波等。以前曾讨论了用圆周卷积计算线性卷积的方法归纳如下:3线性卷积的FFT算法1013线性卷积的FFT算法

将长为N2的序列x(n)延长到L,补L-N2个零，将长为N1的序列h(n)延长到L,补L-N1个零，如果L≥N1+N2-1,则圆周(循环)卷积与线性卷积相等,此时,可用FFT计算线性卷积，方法如下:

1023线性卷积的FFT算法a.计算X(k)=FFT[x(n)]b.求H(k)=FFT[h(n)]c.求Y(k)=H(k)X(k)k=0～L-1d.求y(n)=IFFT[Y(k)]n=0～L-1可见，只要进行二次FFT,一次IFFT就可完成线性卷积计算。103计算表明,L>32时,上述计算线性卷积的方法比直接计算线卷积有明显的优越性,因此,也称上述循环卷积方法为快速卷积法。上述结论适用于x(n)、h(n)两序列长度比较接近或相等的情况。1043线性卷积的FFT算法

如果x(n)、h(n)长度相差较多,例如,h(n)为某滤波器的单位脉冲响应,长度有限,用来处理一个很长的输入信号x(n),或者处理一个连续不断的信号，按上述方法,h(n)要补许多零再进行计算,计算量有很大的浪费，或者根本不能实现。为了保持快速卷积法的优越性,可将x(n)分为许多段,每段的长度与h(n)接近,处理方法有两种：3线性卷积的FFT算法1051)

重叠相加法——由分段卷积的各段相加构成总的卷积输出

h(n)x(n)3线性卷积的FFT算法106计算步骤：a.

事先计算滤波器参数H(k)=DFT[h(n)]，N点b.

用N点FFT计算Xi(k)=DFT[xi(n)]c.

Yi(k)=Xi(k)H(k)d.

用N点IFFT求yi(n)=IDFT[Yi(k)]e.

将重叠部分相加1073线性卷积的FFT算法3线性卷积的FFT算法1082)重叠保留

和第一种方法稍有不同，即将上面分段序列中补零的部分不是补零，而是保留原来的输入序列值，这时，如利用DFT实现h(n)和xi(n)的循环卷积，则每段卷积结果中有N1-1个点不等于线性卷积值需舍去。

重叠保留法与重叠相加法的计算量差不多，但省去了重叠相加法最后的相加运算。

3线性卷积的FFT算法109110111相关的概念很重要，互相关运算广泛应用于信号分析与统计分析，如通过相关函数峰值的检测测量两个信号的时延差等。两个长为N的实离散时间序列x（n）与y（n）的互相关函数定义为

4用FFT计算相关函数则可以证明，rxy(m)的离散傅里叶变换为

Rxy(k)=DFT[rxy(m)] =X*(k)Y

(k)0≤k≤N-1

其中X(k)=DFT[x(n)],Y(k)=DFT[y(n)],1124用FFT计算相关函数113证明：互相关函数定义为

x(n)及y(n)的卷积公式

相比较，可以得到相关和卷积的时域关系：

4用FFT计算相关函数114因得证毕。4用FFT计算相关函数115

当x(n)=y(n)时，得到x(n)的自相关函数为：

维纳--辛钦定理：

自相关函数与信号功率谱互为傅立叶变换对。

4用FFT计算相关函数116推导表明，互相关和自相关函数的计算可利用FFT实现。由于离散傅里叶变换隐含着周期性，所以用FFT计算离散相关函数也是对周期序列而言的。直接做N点FFT相当于对两个N点序列x(n)、y(n)作周期延拓，作相关后再取主值（类似圆周卷积）。

4用FFT计算相关函数117实际一般要求的是两个有限长序列的线性相关，为避免混淆，需采用与循环卷积求线性卷积相类似的方法，先将序列延长补0后再用上述方法。4用FFT计算相关函数118设x(n)和y(n)的长均为N，求线性相关；为了使两个有限长序列的线性相关可用其圆周相关代替而不产生混淆，选择周期L≥2N-1，且L=2m，以使用FFT，将x(n)，y(n)补零至长为L。用FFT计算X(k)，Y(k)（k=0，1…，L-1）。利用FFT求两个有限长序列的线性相关：4用FFT计算相关函数119R（k）=X*（k）Y（k）对R（k）作IFFT，取后N-1项，得取前N项，得4用FFT计算相关函数123

二维信号有图像信号、时空信号、时频信号等。二维离散傅里叶变换可用于处理二维离散信号。二维离散傅里叶变换的定义为：

5用FFT计算二维离散的傅里叶变换124

一个二维DFT可用二个一维DFT计算。

若两次DFT都用快速算法求得,且M和N都是2的幂,则可使用基-2FFT算法,所需要乘法次数为

5用FFT计算二维离散的傅里叶变换125二维离散傅里叶变换所需的乘法次数为（r+q）MN，当M和N比较大时用FFT运算