数字信号处理原理与实践 第4章

可以按上述方法继续分解，将N/2点的DFT 按奇偶部分分解为两个N/4点的DFT，进一 步减少复数乘法的运算量 。

以N＝8为例，将一个8点的DFT分解为两个 4点的DFT，其分解过程如图所示

同样，再把每个4点的DFT分解为两个点的 DFT，对于两点的DFT来说，已经不需要再 分解。因此得到以下结果：






例 4-3 fft2的应用。 X=[2 3 7 5;2 4 6 5;2 4 5 6;1 2 3 4]; Y=fft2(X); 运行结果： Y= 61.0000 -14.0000 + 7.0000i -5.0000 14.0000 - 7.0000i 0 - 7.0000i -3.0000 + 2.0000i 4.0000 - 1.0000i 1.0000 + 2.0000i 7.0000 -2.0000 + 1.0000i 1.0000 -2.0000 1.0000i 0 + 7.0000i -1.0000 - 2.0000i 4.0000 + 1.0000i 3.0000 - 2.0000i
运算可以用蝶形信号流图表示。

求N点的值，需要N2次复数乘法和N(N-1)次 复数加法。对于N/2点DFT有(N/2)2=N2/4次 复数乘法和(N/2)( N/2-1)次复数加法，因此， 两个N/2点DFT有N2/2次复数乘法和N( N/21)次复数加法。

把两个N/2的DFT合成一个N点DFT时，需要 N/2个蝶型运算，包括N/2次复数乘法和N 次复数加法。由于N>>1，总共的复数乘法 为，复数加法次数为，通过这样的分解， 运算量几乎可以减少到一半。





1．fft和ifft函数 调用方式 (1)Y＝fft(x) 参数说明 如果x是向量，则采用傅里叶变换来求解x的离散傅里叶变换； 如果x是矩阵，则计算该矩阵每一列的离散傅立叶变换； (2)Y=fft(x，N) 参数说明 N是进行离散傅立叶变换的x的数据长度，可以通过对x进行补零或截取来获得。 (3)Y=fft(x，[ ]，dim)或Y＝fft(x，N，dim) 参数说明 在参数dim指定的维上进行离散傅立叶变换； 当x为矩阵时，dim用来指定变换的实施方向:dim=1，表明变换按列进行， dim=2表明变换按行进行。 ifft的参数应用与fft完全相同。

IFFT的流图如图
4.6 分裂基FFT算法
分裂基FFT算法又称基2/4算法，或混合基算法，它和前面 讲的基2算法有关，也和基4算法有关。基4算法比基2FFT 算法更快，从理论上说，用较大的基数可以进一步减少运 算次数，但程序过于复杂，所以一般不取基数大于8。 1984年，法国的杜梅尔和霍尔曼首先提出了“分裂基” 算法，把基2算法和基4算法结合起来，其基本思路是对偶 数序号输出采用基2算法，对奇数数序号输出采用基4算法。 由于分裂基算法是目前已知的对于N＝2M算法中具有最少 的加法次数和乘法次数，并且具有非常好的结构，因此被 认为是最好的FFT算法。后来的研究表明，该算法最接近 理论上的所需乘法次数的最小值。首先介绍一下基4FFT算 法。
4.5 IDFT的快速计算方法IFFT

IDFT的快速计算方法是快速傅里叶反变换， 用英文缩写IFFT来表示，无论时选法还是频 选法的FFT算法均可用于IDFT运算。对于N 点序列，其DFT运算公式为

其IDFT运算公式为
只要将DFT运算中的改为，在将结果乘以1/N，就可以 用FFT的计算方法来实现IFFT。频域序列为输入，时域序列 为输出，所以时选FFT对应时选的IFFT，频选FFT对应频选 IFFT。对于系数1/N要做如下处理，把1/N分解为（1/2） M，在M的每一级运算都要乘以1/2，这样做的目的是为了 避免在运算过程中出现溢出。其中N为IFFT变换的点数，N ＝2M，且M为整数。

以N＝8为例，其自然序号是0，1，2，3， 4，5，6，7。第一次按奇偶分开，得到两 组N/2点的DFT，的序列号是0，2，4，6和 1，3，5，7。对每一组再做奇偶分开，这 时抽取后得到四组，每组的序号是0，4和2， 6和1，5和3，7。

这种从十进制看来很乱的输入顺序,实际上 是按二进制“倒序位”排列的。仍然以N=8 为例,x(0),x(1),…,x(7)对应是 x(000),x(001),x(010),x(011),x(100), x(101),x(110),x(111),将二进制码翻转得到 x(000),x(100),x(010),x(110), x(001),x(101),x(011),x(111),它们对应的十进 制序号分别是 x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)。
4.1 离散傅里叶变换存在的问题
对于N点序列，其DFT为
X (k )

n 0
N 1
kn x(n)W N
kn x ( n) 和 WN 一般情况下， 都是复数，每计算一个值，需要作N 次复数乘法和(N-1)次复数加法，求N点的值，需要N2次复 数乘法，以及N(N-1)次复数加法。由于实现一次复数乘法 需要四次实数乘法和两次实数加法，一次复数加法需要两 次实数加法。
以N=16为例,X(k)的第一级的分解可以得到4 个分裂基,X1(k)的第二级分解可以得到2个分 裂基,一个基-4的4点DFT和2个基-2的2点DFT, 而X2(k)和X3(k)的第二级的分解分别是基-4的 4点DFT,对N =16 的分裂基FFT 的示意图如 图所示。
4.7快速傅里叶变换的MATLAB程序实现
4.4 运算量进一步减少的方法
虽然前面讨论的FFT算法，已经大大减 少了运算量，但是还可以通过一些措施， 进一步减少运算量。 对于N＝2M，一共需要M级运算，每级 N/2个蝶型单元，每个蝶型单元进行一次复 数乘法，则总共需要MN/2次复数乘法，如 果＝1，－1，j，-j，则与这些因子相乘时不 作实际的复数乘法运算。

MATLAB为计算快速傅里叶变换，提供了一 系列丰富的数学函数，主要有fft，Ifft，fft2， ifft2等。当所处理的数据的长度为2的幂时， 采用基—2算法进行计算，计算速度会显著 增加。所以，要尽可能使所要处理的数据 长度为2的整数次幂或者用添零方式来添补 数据使其长度成为2的整数次幂。






例4-1 fft(x)函数的应用。 X=[1 4 2 3] Y=fft(x) 运行结果： Y = 10.0000 -1.0000 - 1.0000i -4.0000 1.0000i 例 4-2 ifft(x)函数的应用。 Y = [10.0000 -1.0000 - 1.0000i -4.0000 1.0000i] X=ifft(Y) 运行结果： X=1 4 2 3

当N＝8时，将一个N/2点的DFT分解成两个 N/4个点的DFT的分解过程 如图所示。

这样一个8点的DFT就分解成为四个2点的 DFT 。

如果N＝16，32，或者2的更高次幂，则可 以按上述方法继续分解下去，直到只有两 点DFT为止。
4.2.2 算法的讨论

用上面所述的方法，任何一个N＝2M点的 DFT，都可以经过M次分解，最后分成两个 点的DFT运算，如图所示，每分一次，称为 一“级”运算，所以N点的DFT可以分为M 级。
则
k X (k ) X 1 (k ) W N X 2 (k )
k=0，1，…N/2-1
N / 2k 后半部分的 X ( N / 2 k ) X1 ( N / 2 k ) WN X 2 (N / 2 k)
k X1(k ) WN X 2 (k )

这样，N点 X (k ) 的就可以用X1(k ) 和 X 2 (k ) 来表示， 分别为N/2个点的DFT。
-1.0000 + -1.0000 +

2．fft2和ifft2函数 调用方式 (1)Y＝fft2(x) 参数说明 如果x是向量，则与一维快速傅里叶变换fft相同。 如果x是矩阵，则计算该矩阵的二维快速傅里叶变 换，数据二维傅里叶变换fft2(x)相当于fft(fft(x))， 即先对x的列做一维傅里叶变换，然后对变换结果 的行做傅里叶变换。

DFT运算分为对应的两组，于是有
N 1 2 r 0 N 1 2 r 0
X (k )
x(n)W
n 0
N 1 2 r 0
N 1
nk N

x(2r )W
N 1 2 r 0
2 rk N

x(2r 1)W
( 2 r 1) k N
rk k rk x(2r )WN W x ( 2 r 1 ) W /2 N N /2
4.2.1 算法的推导

设N＝2M，M为正整数，N是2的M次方，如 果不满足这个条件，可以人为的加上若干 个零值点，达到这个要求。这种N为2的整 数次幂的FFT称为基2 FFT。

将 x(n) 按照奇数和偶数分解为两个序列，将N 点的DFT运 算分解为两个N/2点的DFT运算，即令n=2r，n=2r+1，r=0， 1，2，….N/2-1

如果N＝16，通过上述推导可以把16点的 DFT分成4个4点的DFT，其信号第一级流图 如图所示。
可以用同样的方法对X1(k)、X2(k)、X3(k)做第 三级分解,把r 为偶序号的x1(r)作N/4点的DFT, 把r 为奇序号的x1(r)作N8点的DFT。同样的, 对X2(k)、X3(k)做同样的处理。
第四章 快速傅里叶变换
数字信号处理原理与实践 (修订版)@NEU
离散傅里叶变换（DFT）在数字信号 处理中是最常用的运算之一，由于DFT计 算比较繁琐，在相当长的时间内并没有得 到真正的应用。直到1965年库利 (J.W.Cooley)和图基(J.W.Tukey)提出了DFT的 一种快速算法，又有桑德(G.Sande)和图基 的快速算法相继出现，之后又出现了各种 各样的DFT算法，这些方法统称为快速傅 里叶变换（FFT）。


由于每一级都含有N/2个蝶型单元，每一个 蝶型单元又只需要一次复数乘法，两次复 数加法，则完成一次时选FFT的总计算量为： 复数乘法次数：Mf = (N/2)log2N=MN/2 复数加法次数：Af = Nlog2N=MN
当全部FFT完成后，变换后输出的序列依照 正序排列，即存储单元中顺序放着 X(0),X(1),…,X(N-1)。但输入却不是原来的自 然顺序，而是x(0),x(4),x(2),…,x(N-1) 是由于作奇偶分开所产生的。

式中将X(k)分解为偶数组和奇数组,分别令 k=2r,k=2r+1,r=0,1,2,…,N/2-1,则
令
得到两个N/2点DFT
Βιβλιοθήκη Baidu

这样，一个N点的DFT被分解为两个N/2点 的DFT。当N＝8时，上述过程如图所示：

将一个N点DFT分解为4个N/4点DFT的过程 如图所示。
是一个完整的N＝8，按频率抽选的FFT结构如 图所示。
令 X 1 (k )
x (r )W
1 r 0
N 1 2
rk N /2

x(2r )W
r 0
N 1 2
rk N /2
k=0，1，…N/2-1
X 2 (k )
x
r 0
N 1 2
rk WN 2 (r ) /2

x(2r 1)W
r 0
N 1 2
rk N /2
k=0，1，…N/2-1
4.3 按频率抽取基2 FFT算法
与按时间抽取相对应，常用的FFT算法还有 基2频率抽取（decimation-in-frequency,缩写 为DIF）FFT算法，它不是将按时序的奇偶分 解，而是将频域的序号按奇偶分开，下面 说明频率抽取的基2FFT算法。 令序列点数N＝2M，M为整数，先将N点的 DFT分成上下两个部分，即
在实际运算过程中，运算量并没有如此之大， 因为在DFT运算中包含着大量的重复运算， 同时一些系数是不需要乘法运算的。
4.2 按时间抽取的基2 FFT算法

nk W 利用 N 的周期性，对称性，把长度为N点的
DFT运算逐次分解为较短序列的DFT运算。 这种方法叫时间抽取法(decimation-in-time, 缩写为DIT)，简称时选法。