详解FFT(快速傅里叶变换FFT.

第四章 快速傅里叶变换

有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长 序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换 (FFT). 1965 年，Cooley 和 Tukey 提出了计算离散傅里叶变换（DFT）的快 速算法，将 DFT 的运算量减少了几个数量级。从此，对快速傅里叶变换（FFT） 算法的研究便不断深入，数字信号处理这门新兴学科也随 FFT 的出现和发 展而迅速发展。根据对序列分解与选取方法的不同而产生了 FFT 的多种算 法，基本算法是基２DIT 和基２DIF。FFT 在离散傅里叶反变换、线性卷积 和线性相关等方面也有重要应用。

快速傅里叶变换（FFT）是计算离散傅里叶变换（DFT）的快速算法。 DFT 的定义式为

N −1

kn

X (k ) = ∑ x(n)WN

n =0

RN (k )

在所有复指数值 W Nkn 的值全部已算好的情况下，要计算一个 X (k ) 需要 N 次复数乘法和 N－1 次复数加法。算出全部 N 点 X (k ) 共需 N 次复数乘法

2

和 N ( N − 1) 次复数加法。即计算量是与 N 成正比的。

2

FFT 的基本思想：将大点数的 DFT 分解为若干个小点数 DFT 的组合， 从而减少运算量。

WN 因子具有以下两个特性，可使 DFT 运算量尽量分解为小点数的 DFT

运算：

= W N = W N

（1） 周期性： W ( k + N ) n

kn

N

( n + N ) k

（2） 对称性：W

( k + N / 2 )

= −W k

N

N

利用这两个性质，可以使 DFT 运算中有些项合并，以减少乘法次数。例子： 求当 N＝4 时，X(2)的值

3

∑

n =0

4

4

4

4

4

X (2) = x(n)W 2 n = x(0)W 0 + x(1)W 2 + x(2)W 4 + x(3)W 6

= [ x(0) + x(2)]W 0 + [ x(1) + x(3)]W 2

4

0＝[ x(0) + x(2)]－[ x(1) + x(3)]W 4

4

(周期性)

(对称性)

通过合并，使乘法次数由 4 次减少到 1 次，运算量减少。

FFT 的算法形式有很多种，但基本上可以分为两大类：按时间抽取 （DIT）和按频率抽取（DIF）。

4.1 按时间抽取（DIT）的 FTT

为了将大点数的 DFT 分解为小点数的 DFT 运算，要求序列的长度 N 为 复合数，最常用的是 N = 2 的情况（M 为正整数）。该情况下的变换称为

M

基 2FFT。下面讨论基 2 情况的算法。 先将序列 x(n)按奇偶项分解为两组

⎧x(2r ) = x1 (r ) ⎨

1) = x2 (r ) ⎩x(2r +

将 DFT 运算也相应分为两组

N −1

N

r = 0,1,L, − 1

2

X (k ) = DFT [ x(n)] = ∑ x(n)W kn

N

n =0

N −1

N

N −1

N

= ∑ x(n)W kn＋ ∑ x(n)W kn

n=0 n为偶数

n =0 n 为奇数

N / 2−1

∑

(2 )

2 rk

+

N / 2−1

∑ (2

+ 1)

r W

( 2 r +1) k

x r

W＝

x N

N

r =0

r =0

N / 2−

1

2 rk

∑+

N / 2−1 k

∑

2 rk

＝

x1 (r )WN

WN x2 (r )WN

r =0

r =0

N / 2−1

k

rk + WN / 2−1

＝ ∑ N / 2

N

x1 (r )W

∑

rk

x2

(r )WN / 2

r =0

r =0

= X 1 (k ) + Wk

N X 2 (k )

W

2 rk N

= W rk

N / 2 ）

（因为

其中 X 1 (k ) 、 X 2 (k ) 分别是 x1 (n)、x2 (n) 的 N/2 点的 DFT

N / 2−1

rk

N / 2−1

rk

N

X 1 (k ) ＝ ∑ x1 (r )WN / 2 = ∑ x(2r )WN / 2 ,0 ≤ k ≤ − 1

r =0

r =0

2

rk

N / 2 −1

rk

N / 2 −1

N

X 2 (k ) ＝ ∑ x 2 (r )W N / 2 = ∑ x(2r + 1)W N / 2 ,0 ≤ k ≤ −1

r =0

r =0

2

至此，一个 N 点 DFT 被分解为两个 N/2 点的 DFT。

上面是否将全部 N 点的 X (k ) 求解出来了？

分析： X 1 (k ) 和

N − 1 ），则 由

X 2 (k ) 只有 N/2 个点（ k = 0,1,L,

2

k

X (k ) = X (k ) + W X (k ) 只能求出 X (k ) 的前 N/2 个点的 DFT，要求出

1 N 2

全部 N 点的 X (k ) ，需要找出 X 1 (k ) 、 X 2 (k ) 和 X (k + N / 2) 的关系，其

N

中 k = 0,1,L, − 1。由式子 X (k ) = X 1 2

k

(k ) + W NX 2 (k ) 可得

k + N / 2X (k + N / 2) = X (k + N / 2) + W X 2 (k + N / 2) 化简得

1 N

k

N

X (k + N / 2) ＝ = X 1 (k ) − WN X 2 (k ) ， k = 0,1,L, − 1 2

这样 N 点 DFT 可全部由下式确定出来：

X (k ) = X (k ) + W X (k ) ⎪ 1 N 2

⎨ N k 2

X (k + N / 2) = X 1 (k ) − W X (k )

算。

k

k = 0,1,L , N − 1 （＊）

2

上式可用一个专用的碟形符号来表示，这个符号对应一次复乘和两次复加运

a

a + W k b

N

b

W

k

a − WN k b

-1

N

图 蝶形运算符号

N 2 N 2

次复数乘 通过这样的分解以后，每一个 N／2 点的 DFT 只需要 ( ) = 2

4

N 2 N 2

次复乘，再加上将两个 N／2 点 法，两个 N/2 点的 DFT 需要 2( ) = 2

2

DFT 合并成 为 N 点 DFT 时有 N ／ 2 次与 W 因 子相 乘，一 共需 要

N 2 2

N

+ 2

N 2

≈ 次复乘。可见，通过这样的分解，运算量节省了近一半。 2

因为 N = 2 M ，N/2 仍然是偶数，因此可以对两个 N/2 点的 DFT 再 分别作进一步的分解，将两个 N/2 点的 DFT 分解成两个 N/4 点的 DFT。

例如对 x1 (r) ，可以在按其偶数部分及奇数部分进行分解：

⎧x1 (2l ) = x3 (l ) ⎨

1) = x4 (l ) ⎩x1 (2l +

则的运算可相应分为两组：

N / 4−1

N

l = 0,1,L, − 1

4

2lk N / 2

N / 4−1

( 2l +1) k N / 2

X 1 (k ) ＝ ∑ x1 (2l )W

l =0

+ ∑ x1 (2l + 1)W

l =0

N / 4−1

∑＝

3

( )

lk

/ 4

+

k N / 2

N / 4−1

∑

l =0

( )

4

lk N / 4

x l WN

l =0

W

x l W

N

k = 0,1,L,

2 k

k

= X 3 (k ) + WN / 2 X 4 (k )

4

− 1

k将系数统一为以Ｎ为周期，即W N / 2 = W N ，可得

) X 1 (k ) = X ) + W 3 (k N X 4 (k ⎪

2 k

⎨

2 kk = 0,1,L , N − 1

X 1 (k + N / 4) = X 3 (k ) − WN X 4 (k )

4

同样，对 X 2 (k ) 也可进行类似的分解。一直分解下去，最后是２点的

DFT，２点 DFT 的运算也可用碟形符号来表示。这样，对于一个 N = 2 = 8

3

的 DFT 运算，其按时间抽取的分解过程及完整流图如下图所示。

这种方法，由于每一步分解都是按输入序列在时域上的次序是属于偶数 还是奇数来抽取的，故称为“时间抽取法”。

分析上面的流图， N = 2 M ，一共要进行 M 次分解，构成了从 x(n)到 X(k)的 M 级运算过程。每一级运算都是由 N/2 个蝶形运算构成，因此每一 级运算都需要 N/2 次复乘和 N 次复加，则按时间抽取的 M 级运算后总共需 要

复数乘法次数： mF = N

N ⋅ M = log N

2

2

2

复数加法次数： aF = N ⋅ M = N log 2 N

根据上面的流图，分析ＦＦＴ算法的两个特点，它们对ＦＦＴ的软硬件 构成产生很大的影响。

（１） 原位运算 也称为同址运算，当数据输入到存储器中以后，每一级运算的结果仍然存储 在原来的存储器中，直到最后输出，中间无需其它的存储器。根据运算流图 分析原位运算是如何进行的。原位运算的结构可以节省存储单元，降低设备 成本。

（２） 变址 分析运算流图中的输入输出序列的顺序，输出按顺序，输入是“码位倒置” 的顺序。见图。

自然顺序 0 1 2 3 4 5 6 7 0

二进制表示 000 001 010 011 100 101 110 111 码位倒置顺序 1 2 3 码位倒置 000 100 010 110 001 101 011 111 4 5 码位倒置顺序 0 4 2 6 1 5 3 7 6 7

X(0) X(4) X(2) X(6) X(1) X(5) X(3) X(7) 码位倒置的变址处理

在实际运算中，直接将输入数据 x(n)按码位倒置的顺序排好输入很不

方便，一般总是先按自然顺序输入存储单元，然后通过变址运算将自然顺序

的存储换成码位倒置顺序的存储，这样就可以进行 FFT 的原位运算。变质 的功能如图所示。用软件实现是通用采用雷德（Rader）算法，算出 I 的倒 序Ｊ以后立即将输入数据 X(I)和 X(J)对换。尽管变址运算所占运算量的比例 很小，但对某些高要求的应用（尤其在实时信号处理中），也可设法用适当 的电路结构直接实现变址。例如单片数字信号处理器 TMS320C25 就有专用 于 FFT 的二进制码变址模式。

4.２ 按频率抽取（DIF）的 FTT

除时间抽取法外，另外一种普遍使用的 FFT 结构是频率抽取法。频率 抽取法将输入序列不是按奇、偶分组，而是将Ｎ点 DFT 写成前后两部分：

N −1

X (k ) = DFT [ x(n)] = ∑ x(n)W kn

N

n =0

( N / 2)−1

N −1

= ∑ x(n)W kn＋ ∑ x(n)W kn

n =0

N

n= N / 2

N

N / 2−1

∑ ( )＝

n =0

nk

N / 2−1

+∑ ( +n=0

/ 2)

( n + N / 2) k

x n WN

x n N

WN

nk

N / 2−1

( N / 2) k

＝ ∑[ x(n) + WN

n =0

x(n + N / 2)]WN

因为 W

N / 2 N

= −1,W ( N / 2) k = (−1) k ， k 为 偶数时 (−1) k = 1 ， k 为奇数时

N

(−1) k = −1，由此可将 X(k)分解为偶数组和奇数组：

N / 2−1

∑

N / 2−1

N

X (k )＝ n =0 [ x(n) + (−1) k x(n + N / 2)]W nk

∑

N / 2−1

N

X (2r )＝ n=0 [ x(n) + x(n + N / 2)]W 2 nr

∑

=

n =0

N / 2

[ x(n) + x(n + N / 2)]W nr

N / 2−1

∑

X (2r + 1)＝

N / 2−1

n =0

N

[ x(n) − x(n + N / 2)]W ( 2 r +1) n

∑

=

n =0

N

N / 2

[ x(n) − x(n + N / 2)]W nW nr

⎧x1 (n) = x(n) + x(n + N / 2)

令 ⎨

N

n = 0,1,L, N / 2 − 1

⎩x2 (n) = [ x(n) − x(n + N / 2)]W n

这两个序列都是 N/2 点的序列，对应的是两个 N/2 点的 DFT 运算：

N / 2−1

∑

X (2r )＝

n =0

N / 2

1

x (n)]W nr

N / 2−1

∑

X (2r + 1)＝

n =0

N / 2

2

x (n)W rn

这样，同样是将一个 N 点的 DFT 分解为两个 N/2 点的 DFT 了。频率抽选法 对应的碟形运算关系图如下：

a

a + b

b

-1

(a − b)W n

N

WN

n

对于 N=8 时频率抽取法的 FFT 流图如下：

这种分组的办法由于每次都是按输出 X(k)在频域的顺序上是属于偶数还是 奇数来分组的，称为频率抽取法。与前面按时间抽取的方法相比，相同点 问题：如何利用快速算法计算 IDFT？ 分析 IDFT 的公式：

−nk x(n) = IDFT [ X (k )] = 1 ∑ X (k )W , n = 0,1,L, N − N

N 1k =0

N −1

比较 DFT 的公式：

N −1

X (k ) = DFT [ x(n)] = ∑ x(n)W nk , k = 0,1,L, N − 1

N

n =0

得知可用两种方法来实现 IDFT 的快速算法：（１）只要把 DFT 运算中的每

一个系数 W nk 该为 W − nk ，并且最后再乘以常数 N

N

1

，就可以用时间抽取法

N

或频率抽取的 FFT 算法来直接计算 IDFT。这种方法需要对 FFT 的程序和参 数 稍 加

N −1 1 1 x(n) = {DFT [ X * (k )]}, n = 0,1,L, N − 1 ，也就 * Nnk * [∑ X (k )W ]= N k =0 N

改 动 才 能 实 现 。 （ ２ ） 因 为

是说，可先将 X(k)取共轭变换，即将 X(k)的虚部乘以－１，就可直接调用 FFT 的程序，最后再对运算结果取一次共轭变换并乘以常数 1/N 即可得到 x(n)的值。这种方法中，FFT 运算和 IFFT 运算都可以共用一个子程序块， 在使用通用计算机或用硬件实现时比较方便。

4.1.3 混合基 FFT 算法

以上讨论的是基 2 的 FFT 算法，即 N = 2 的情况，这种情况实际

M

上使用得最多，这种 FFT 运算，程序简单，效率很高，用起来很方便。另 外，在实际应用时，有限长序列的长度 N 到底是多少在很大程度时是由人 为因素确定的，因此，大多数场合人们可以将 N 选定为 N = 2 ，从而可

M

以直接调用以２为基数的 FFT 运算程序。

如果长度 N 不能认为确定，而 N 的数值又不是以 2 为基数的整数次 方，一般可有以下两种处理方法：

（１）

将 x(n)用补零的方法延长，使 N 增长到最邻近的一个

N = 2 M 数值。例如，N=30，可以在序列 x(n)中补进

x(30)＝x(31)＝０两个零值点，使 N=32。如果计算 FFT 的目的是为了了解整个频谱，而不是特定频率点，则此 法可行。因为有限长序列补零以后并不影响其频谱

X (e jw ) ，只是频谱的采样点数增加而已。

（２） 如果要求特定频率点的频谱，则 N 不能改变。如果 N 为复合数，则可以用以任意数为基数的 FFT 算法来计 算。快速傅里叶变换的基本思想就是要将 DFT 的运算

尽量分小。例如，N=6 时，可以按照 N=3×２分解， 将６点 DFT 分解为３组２点 DFT。

举例：N=9 时的快速算法。

4.2 快速傅里叶变换的应用

凡是可以利用傅里叶变换来进行分析、综合、变换的地方，都可以利

用 FFT 算法及运用数字计算技术加以实现。FFT 在数字通信、语音信号处 理、图像处理、匹配滤波以及功率谱估计、仿真、系统分析等各个领域都得 到了广泛的应用。但不管 FFT 在哪里应用，一般都以卷积积分或相关积分 的具体处理为依据，或者以用 FFT 作为连续傅里叶变换的近似为基础。

4.2.1 利用 FFT 求线性卷积—快速卷积

在实际中常常遇到要求两个序列的线性卷积。如一个信号序列 x(n)

通过 FIR 滤波器时，其输出 y(n)应是 x(n)与 h(n)的卷积：

y(n) = x(n) ∗ h(n) =

∞

x(m)h(n − m) ∑

m =−∞

有限长序列 x(n)与 h(n)的卷积的结果 y(n)也是一个有限长序列。假设 x(n)与 h(n)的长度分别为 N1 和 N2，则 y(n)的长度为 N=N1+N2-1。若通过补零使 x(n)与 h(n)都加长到 N 点，就可以用圆周卷积来计算线性卷积。这样得到用 FFT 运算来求 y(n)值（快速卷积）的步骤如下：

（１） 对 序 列 x(n) 与 h(n) 补 零至长为 N ，使 N ≥ N1+N2-1 ，并且

N = 2 M (M为整数) ，即

1 ⎧x(n), n = 0,1,L, N1 −

⎨ x(n) = ⎩0, n = N1, N1 + 1,L, N − 1

1 ⎧h(n), n = 0,1,L, N 2 − h(n) = ⎨

⎩0, n = N 2, N 2 + 1,L, N − 1

（２） 用 FFT 计算 x(n)与 h(n)的离散傅里叶变换

4.2 快速傅里叶变换的应用

h(n) ⇔ ( FFT ) ⇔ H (k )

x(n) ⇔ ( FFT ) ⇔ X (k )

（N 点）

（N 点）

（３） 计算 Y(k)=X(k)H(k)

（４） 用 IFFT 计算 Y(k)的离散傅里叶反变换得：y(n)=IFFT[Y(k)] （N 点）

4.2.2 利用 FFT 求相关—快速相关

互相关及自相关的运算已广泛的应用于信号分析与统计分析，应用于连 续时间系统也用于离散事件系统。

用 FFT 计算相关函数称为快速相关，它与快速卷积完全类似，不同的 是一个应用离散相关定理，另一个应用离散卷积定理。同样都要注意到离散 傅里叶变换固有的周期性，也同样用补零的方法来绕过这个障碍。

设两个离散时间信号 x(n)与 y(n)为已知，离散互相关函数记作 Rxy (n) ， 定义为

∞

Rxy (n) =

x(m) y(n + m) ∑

m =−∞

如果 x(n)与 y(n)的序列长度分别为 N1 和 N2，则用 FFT 求相关的计算步骤 如下：

（ １ ）对序 列 x(n) 与 y(n) 补 零至长为 N ，使 N ≥ N1+N2-1 ，并且

N = 2 M (M为整数) ，即

1 ⎧x(n), n = 0,1,L, N1 −

x(n) = ⎨

⎩0, n = N1, N1 + 1,L, N − 1

y(n), n = 0,1,L, N 2 − 1 ⎧ ⎨ y(n) = ⎩0, n = N 2, N 2 + 1,L, N − 1

（２）用 FFT 计算 x(n)与 y(n)的离散傅里叶变换

x(n) ⇔ ( FFT ) ⇔ X (k )

（N 点）

y(n) ⇔ ( FFT ) ⇔ Y (k )

（N 点）

（３）将 X(k)的虚部 Im[X(k)]改变符号，求得其共轭 X*(k)

4.2.2 利用 FFT 求相关—快速相关

（５） 用 IFFT 求得相关序列 Rxy (n)

（４）计算 Rxy (k ) =X*(k)Y(k)

Rxy (n) ＝IFFT[ Rxy (k ) ]

（N 点）

如果 x(n)＝y(n)，则求得的是自相关序列 Rxx (n)

4.2.3 Chirp－Z 变换

采用 FFT 算法可以很快的计算出全部 DFT 值，也即 Z 变换在单位圆 上的全部等间隔采样值。但是，很多场合下，并非整个单位圆上的频谱都是 很有意义的，例如对于窄带信号过程，往往只需要对信号所在的一段频带进 行分析，这是，希望采样能密集在这段频带内，而对频带以外的部分，则可 以完全不管。另外，有时也希望采样能不局限于单位圆上，例如，语音信号 处理中，往往需要知道其 Z 变换的极点所在频率，如果极点位置离单位圆 较远，，则其单位圆上的频谱就很平滑，如图所示，这是很难从中识别出极 点所在的频率。要是采样不是沿单位圆而是沿一条接近这些极点的弧线进 行，则所得的结果将会在极点所在频率上出现明显的尖峰，从而可以较准确 的测定极点频率。

螺线采样就是一种适用于这种需要的变换，并且可以采用 FFT 来快 速计算，这种变换也称为 Chirp-Z 变换，它是沿 Z 平面上的一段螺线作等分 角的采样，这些采样点可表达为

− k

z = AW , k = 0,1,L, M − 1 k

其中M为采样点的总数，A为起始点位置，这个位置可以进一步用它的半径 A0及相角θ 0 来表示

A = A e jθ 0

0

参数 W 可表示为

W = W 0 e jΦ0

其中 W0 为螺线的伸展率， W0 >１时螺线内缩（反时针方向）； W0 <１时螺

线外伸。Φ 螺线采样点在 Z 平面上的分布 0 为螺线上采样点之间的等分角。可表示为下图。

下面分析这些点上采样值计算的特点。假定 x(n)是长度为 N 的有限 长信号序列，则其 Z 变换在采样点上的值为