数字信号处理：第3章3.5-3.7 快速傅里叶变换

1、13.5 快速傅里叶变换(FFT) 2 引言 DFT是信号分析与处理中的一种重要变换。因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换(简称FFT)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。直到1965年发现了DFT的一种快速算法以后，情况才发生了根本的变化。 33.5 基2FFT算法 3.5.1 直接计算DFT的特点及减少运算量的基本途径 长度为N的有限长序列x(n)的DFT为 考虑x(n)为复数序列的一般情况，对某一个k值，直接按(3.2.1)式计算X(k)值需要N次复数乘法、(N-1)次复数加法。共需要N2次复数乘

2、法, N(N-1)次复数加法。(3.2.1) 4在导出FFT算法之前，首先来估计一下直接计算DFT所需的计算量。将式（3.20a）展开成方程组将方程组写成矩阵形式5用向量表示为 式中，X为N维列向量，称为变换列向量，为复数；W为 方阵，称为系数矩阵，为复数；x为N维列向量，表示离散时间信号，可以是复数。即使x是实数，也可用复数表示，于是6从式(3.59)看出，由于计算一个X(k)值需要N次复乘法和(N-1)次复数加法，因而计算N个X(k)值，共需 次复乘法和N(N-1)次复数加法。式(3.61)说明，每次复乘法包括4次实数乘法和2次实数加法，每次复加法包括2次实数加法，因此计算N点DFT共需4

3、 次实数乘法和 次实数加法。当N很大时，这是一个非常大的计算量。FFT算法利用了 的两个性质：1.对称性，即 ；2.周期性，即 ，r为任意整数。图3.14所示的是N=8情况下这两个性质的示意图。 FFT算法有许多种，本书主要介绍时间抽选 (Decimation-in-Time)基2FFT算法和频率抽选 (Decimation-in-Frequency）基2FFT算法，并简要介 绍其他基的快速傅里叶变换算法。前两种算法特别 适用于N等于2的幂的情况，后一种算法适合于N为复合数的 情况。图3.14 的特性 0-17 3.5.2 时域抽取法基2FFT基本原理 FFT算法基本上分为两大类：时域抽取法F

4、FT(Decimation In Time FFT,简称DIT-FFT)和频域抽取法FFT(Decimation In Frequency FFT,简称DIFFFT)。下面先介绍DITFFT算法。 设序列x(n)的长度为N，且满足为自然数 按n的奇偶把x(n)分解为两个N/2点的子序列8 则x(n)的DFT为由于所以 奇数偶数9 其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，即 (3.2.5) (3.2.6) 由于X1(k)和X2(k)均以N/2为周期，且 ，所以X(k)又可表示为(3.2.7) (3.2.8) 10图4.2.1 蝶形运算符号 11图4.2.2 N点D

5、FT的一次时域抽取分解图(N=8) 12 与第一次分解相同，将x1(r)按奇偶分解成两个N/4长的子序列x3(l)和x4(l)，即那么，X1(k)又可表示为 (3.2.9) 13式中 同理，由X3(k)和X4(k)的周期性和 的对称性 最后得到： (3.2.10) 14 用同样的方法可计算出(3.2.11)其中 15图3.2.3 N点DFT的第二次时域抽取分解图(N=8) 16 图3.2.4 N点DITFFT运算流图(N=8) 同址(原位)计算17图3.19 时间抽选的FFT流程图（N8）18 3.5.3 DITFFT算法与直接计算DFT运算量的比较1.蝶形运算及运算量的比较 每级由N/2个蝶

6、形组成,每一级运算都需要N/2次复数乘和N次复数加(每个蝶形需要两次复数加法)。所以，M级运算总共需要的.复数乘次数为复数加次数为 例如，N=210=1024时P84表格19 2.同址计算(原位) 每一级的蝶形的输入与输出在运算前后可以存储在同一地址的存储单元中,这种同址运算的优点是可以节省存储单元,减低对硬件的要求.20 3. 序列的倒序或变址计算 DITFFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种倒序是很有规律的。由于N= ，所以顺序数可用M位二进制数(nM-1nM-2n1n0)表示。 图4.2.7 形成倒序的树状图(N=23) 21表4.2.1 顺序和倒序二进制数对照

7、表 22 图4.2.8 倒序规律 23 3.5.4 频域抽取法FFT(DIFFFT) 在基2快速算法中，频域抽取法FFT也是一种常用的快速算法，简称DIFFFT。 设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式：24偶数 奇数 将X(k)分解成偶数组与奇数组，当k取偶数(k=2r,r=0,1,N/2-1)时 25当k取奇数(k=2r+1,r=0,1,N/2-1)时(3.2.15) 定义x1(n)和x2(n)分别代入(3.2.14)和(3.2.15)式，可得 (3.2.16) 26图3.2.10 DIFFFT蝶形运算流图符号 27图3.2.11

8、 DIFFFT一次分解运算流图(N=8) 28 由于N是2的整数幂，所以N2仍然是偶数。这样，可将N2点DFT的输出再分为偶数组和奇数组，也就是将N2点的DFT计算进一步分解为两个N4点的DFT计算。29图3.2.12 DIFFFT二次分解运算流图(N=8) 30图3.2.13 DIFFFT运算流图(N=8) 31图3.19 时间抽选的FFT流程图（N8）32显然，一个蝶形的计算包括1次复数乘法和2次复数加法。图3.27所示的整个流程图共有 级迭代运算，每级有 个蝶形。因此，总计算量为：复数乘法次数：复数加法次数：频率抽选的FFT算法的计算量与时间抽选FFT算法的计算量是一样的。33 3.5.

9、6 IDFT的高效算法 上述FFT算法流图也可以用于离散傅里叶逆变换(Inverse Discrete Fourier Transform，简称IDFT)。比较DFT和IDFT的运算公式： 比较式(3.84) 和式(3.85)，可以看出，只要把DFT公式中的系数 ， 并乘以系数 ，就可用FFT算法来计算IDFT，这就得到了IFFT的算法。当把时间抽选FFT算法用于IFFT计算时，由于原来输入的时间序列x(n)现在变为频率序列X(k)，原来是将x(n)偶奇分的，现在变成对X(k)进行偶奇分了，因此这种算法改称为频率抽选IFFT算法。类似的，当把频率抽选FFT算法用于计算IFFT时，应该称为时间抽

10、选IFFT算法。34图4.2.16 DITIFFT运算流图 35 另一种方法, 如果希望直接调用FFT子程序计算IFFT，则可用下面的方法： 由于 对上式两边同时取共轭，得36 3.6 N为合数的FFT算法 上面讨论的以2为基(即N= )的时间抽选和频率抽选FFT算法，由于具有程序简单、计算效率高、对存储量要求不很高等优点，因而在实际中得到最广泛的应用。如N不等于2的幂 ，通常有两种处理办法 (1)用补零的办法将x(n)延长为 。例如N=60，可在序列x(n)末尾填补4个0，即令x(60)=x(61)=x(62)=x(63)=0，使N达到 =64，这样就可使用基2FFT算法。有限长序列补零以后

11、，只是频谱的取样点有所增加而不会影响它的频谱 的形状。37（2）采用以任意数为基数的FFT算法。这种算法举例说明如下：设N 等于两个整数p和q的乘积，即N=p*q，则可将N点DFT分解成p个q点DFT或q个p点DFT来计算。为此，首先将x(n)分成p组，每组长为q，即38然后，将N点DFT也分解为p组来计算，即 （3.87）3940组图3.30 一个 点的DFT分组示意图41在最一般的情况下，设 （3.90）其中 是m个素因子。首先把N分解为两个因子，即 , 其中 ，并用以上讨论的方法将DFT分解为 个 点DFT；然后，将 分解为 ，其中 ，即将每一个 点DFT分解为 个 点DFT；这样，通过

12、m次分解，最后达到 点DFT。这种算法可以使DFT的运算获得最高效率。 42图3.31 一个 的DFT流程图433.7 FFT 的应用 DFT的快速算法FFT的出现， 使DFT在数字通信、 语言信号处理、 图像处理、 功率谱估计、 仿真、 系统分析、 雷达理论、 光学、 医学、 地震以及数值分析等各个领域都得到广泛应用。3.7.1 四种傅立叶变换1、非周期连续时间信号的傅立叶变换（FT）4445 2、周期连续时间信号的傅立叶级数（FS）3、序列傅立叶变换46 4、周期序列的傅立叶级数（DFS）注意周期性、离散性在时域和频域中的对称47图3.7.1 四种傅里叶变换示意图48表3.7.1 四种傅里

13、叶变换的特点49 3.7.2 用FFT对信号进行谱分析 所谓信号的谱分析就是计算信号的傅里叶变换。 连续信号与系统的傅里叶分析显然不便于直接用计算机进行计算，使其应用受到限制，而DFT是一种时域和频域均离散化的变换，适合数值运算，成为分析离散信号和系统的有力工具。 1. 谱分析中参数的选择 工程实际中， 经常遇到的连续信号xa(t)， 其频谱函数Xa(j)也是连续函数。 5051 采样信号的频谱是原模拟信号的频谱沿频率轴，每间隔采样角频率s重复出现一次，或者说采样信号的频谱是原模拟信号的频谱以s为周期，进行周期性延拓而成，乘以系数1/T。 设xa(t)是带限信号，最高截止频率为0，其频谱Xa(

14、j)如图(a)所示。 理想取样52 _53频率归一化现在讨论离散时间信号 的频谱 与取样信号 的频谱 之间的关系。假设离散时间信号 是模拟信号 通过周期性取样得到的，即其傅里叶变换为54将式（2.63）代入上式得 (2.67) 式（2.67）表明，在 的条件下，离散时间信号 的频谱与取样信号频谱相等。由于 （ 为取样频率）是 对 归一化的结果，故可以认为离散时间信号 的频谱是取样信号的频谱经频率归一化后的结果，如图2.25(c)所示。5556与序列傅立叶变换的关系平面57 58 显然，X (jf)仍是f的连续周期函数, xa(nT)和X (jf)如图3.4.5(b)所示。 对 X(jf)在区间

15、0, fs上等间隔采样N点， 采样间隔为F，F也叫频率分辨率 如图3.4.5(c)所示。 参数fs 、 tp、 N和F满足如下关系式： 由于NT= tp ， tp为信号的观测时间，所以 N越大,F越小,分辨率越高59根据取样定理，为了避免混叠失真，要求最小记录长度必须按所需的频率分辨率来选择，即从式（3.91）可看出，当提高信号最高频率 时，必须减少T；在N固定的情况下，由式（3.92）可看出，记录长度 将缩短，这意味着频率分辨率要降低。反之，若要提高频率分辨率，就必须增加 ，在给定N的情况下，这必然导致T的增加，结果使能分析的最高频率降低。在保持分辨率不变的情况下，若希望增加所分析的信号的最

16、高频率，或在保持信号最高频率不变的情况下，提高分辨率，唯一的办法是增加在记录长度内的取样点数N。如果 和F都给定，那么N必须满足条件602.频谱分析的步骤(1).数据准备 已知(2).用FFT计算频谱61 (3). 计算振幅谱 相位谱 功率谱例62 63例 连续信号 ,式中 ， 的最高频率是20Hz, 现对 进行采样，要求信号的频率分辨率为2 Hz，试计算：1、最小采样频率2、信号的观测时间3、信号的采样点数N，如果用FFT计算频谱，N的取值应该是多少？64 3.7.3 用FFT计算线性卷积 如果0kL-1则由时域循环卷积定理有 Y(k)=FFTy(n)=X1(k)X2(k), 0kL-165

17、 由此可见， 循环卷积既可在时域直接计算， 也可以按照图3.4.1所示的计算框图， 在频域计算。 由于DFT有快速算法FFT， 当N很大时， 在频域计算的速度快得多， 因而常用DFT(FFT)计算循环卷积。 图 3.4.1 用FFT计算循环卷积 FFTFFTIFFT66 两个有限长序列的线性卷积可以用循环卷积来代替，而循环卷积可使用FFT来计算。为使循环卷积与线性卷积计算结果相等，必须把x(n)和h(n)都延长到N点 ，延长的部分均为零取样值。这样，y(n)的计算由以下5步来完成。 （1）将x(n)和h(n)都延长到N点， （2）计算x(n)的N点FFT,即X(k)=FFTx(n) （3）计算

18、h(n)的N点FFT,即H(k)=FFTh(n) （4）计算 （5）计算Y(k)的反变换，即y(n)=IFFTX(k)H(k)67图 3.4.3 用FFT计算线性卷积框图 补N- N1个零点N点FFT补N- N2个零点N点FFTN点IFFTy(n)h(n)x(n)68计算量的比较N1=N281610244096r0.570.9429.2100例:首先讨论 与 的长度相等的情况，这时 ，于是 69其次，讨论信号 相对于单位取样响应很长即 的情况，此时 ，于是显然，r值将下降，从而使循环卷积算法的优点不能发挥。下一节将介绍的分段卷积方法，适合 的情况。70 长序列和短序列的线性卷积直接利用DFT计

19、算的缺点：(1) 信号要全部输入后才能进行计算，延迟太多(2) 内存要求大(3) 算法效率不高解决问题方法：采用分段卷积分段卷积可采用重叠相加法 和 重叠保留法71 设序列h(n)长度为M， x(n)为无限长序列。 将x(n)均匀分段， 每段长度取L， 则子段之间不存在重叠。注意，在这种分段方式中，xk(n)的时间起点和x(n)的原点一致.于是， h(n)与x(n)的线性卷积可表示为(3.4.4) 72重叠相加法处理步骤73图 3.4.4 重叠相加法卷积示意图 （a）将x(n)分成长度为L的互不重叠的几段（b）每一段与h(n)卷积结果74 在重叠相加法中，虽然子段xk(n) (长度为L)之间没

20、有重叠，但相邻的yk(n)(长度为N=L+M-1)有重叠，最终结果通过对这些有重叠的部分结果求和得到，这也正是重叠相加法名字的由来。*75 方法：将x(n)长序列分段，长度为L，然后将重叠相加法中分段序列尾部补零的部分改为保留原来输入序列对应位置的值； (2) 各段序列xk(n)分别与 M点短序列h(n)进行长度为N的循环卷积；重叠保留法7677 (3) 从各段循环卷积中提取线性卷积结果。因 yk(n)=xk(n) L h(n)前M-1个点不是线性卷积结果中的点,故x(n)分段时，每段与其前一段有M-1个点重叠。78可见，在重叠保留法中，将原始序列进行有重叠的分段，子段与单位脉冲响应之间进行的

21、是循环卷积，每次保留与线性卷积的结果相同的点，通过上式求和得到最终结果，这也是重叠保留法名字的由来。当然，子段xk(n)和h(n)的循环卷积可以借助FFT在频域实现。7980 -x (n(M1)M-1-N (M1)N-1x 0(n)x 1(n)2N-Mn第一段前需补M-1个零81重叠保留法设序列h(n)的长度为M，则对长序列x(n)的分段方法如下：先在序列x(n)前补M-1个0，然后对补0后的序列进行分段，每段的长度为N=M+L-1,即：对每一段 ，通过循环卷积 ，获得二者的线性卷积 。而输入的每段序列重叠M-1点，故每段的循环卷积输出应去掉前面M-1点只保留后面L点，即 ：8283续图 3.3584 例 已知序列x(n)=n+2,0n12, h(n)=1,2,1试用重叠保留法计算线性卷积, 取L=3, 即N=M+L-1 5 85 解: 重叠保留法y(n)=2, 7, 12, 16, 20, 2