FFT为什么这么快？

离散傅里叶变换（DFT）是为了能让计算机处理经过采集的离散信号提出的，DFT的推导以及代码实现在下面的两篇文章中也已经介绍了：

1. DFT推导
2. DFT代码实现

个人认为，DFT是数字信号处理领域最重要的数学思想，而快速傅里叶变换（FFT）是该领域最重要的算法。FFT不是一种新的变换，他是DFT的快速实现。正是因为FFT的提出，才让傅里叶变换算法能够在工程界大显身手。

1. 离散傅里叶变换的矩阵表达

DFT变换的公式为：
X p = ∑ n = 0 N − 1 x n e   − j 2 π N n p , p ∈ { 0 , 1 , … , N − 1 } X_{p}=\sum_{n=0}^{N-1} x_{n} \mathrm{e}^{\:-j \frac{2 \pi}{N} n p}, \quad p \in\{0,1, \ldots, N-1\} \\ Xp​=n=0∑N−1​xn​e−jN2π​np,p∈{0,1,…,N−1}
为了更清楚地计算DFT的算法复杂度，我们将其表示为矩阵表达：
$$\left[\begin{array}{c}{X}_0\\ X_1\\ \cdots \\ X_{N-1}\end{array}\right]=\left[\begin{array}{ccccc}{W}_N^0& W_N^0& W_N^0& \cdots & W_N^0\\ W_N^0& W_N^1& W_N^2& \cdots & W_N^{N-1}\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ W_N^0& W_N^{N-1}& W_N^{2(N-1)}& \cdots & W_N^{(N-1{)}^2}\end{array}\right]\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{N-1}\end{array}\right]$$
式中：

$$W_N=e^{-j\frac{2\pi }{N}}$$

2. DFT复杂度计算

因为 $N\ast N$ 的变换矩阵中的每一个元素可以离线计算并存储，这样上述计算过程仅涉及 $N^2$ 次复数乘法运算和 $N(N-1)$ 次复数加法运算。

因为变换矩阵的第一行和第一列都是 $1$ ，这样就节约了 $2N-1$ 次的复数乘法运算，所以实际上，上述计算过程只需要进行 $N^2-(2N-1)$ 次复数乘法运算，即共需要 $(N-1{)}^2$ 次复数乘法运算和 $N(N-1)$ 次复数加法运算。

因为 $1$ 次复数乘法运算相当于 $4$ 次实数乘法运算和 $2$ 次实数加法运算，$1$ 次复数加法运算相当于 $2$ 次实数加法运算，所以：

一共需要进行 $4(N-1{)}^2$ 次实数乘法运算和 $2(N-1{)}^2+2N(N-1)=4(N-0.5)(N-1)$ 次实数加法运算。

所以离散傅里叶变换算法复杂度是 $O(n^2)$ ，如果要处理的数据点数过多，就会需要大量的计算时间，这会给工程应用带来很大的限制。

3. FFT算法

Cooley和Tukey在1965年宣布，他们“第一次”发现了FFT算法。而实际上，他们发现的FFT算法的核心原理，早在150多年前（1806年），就被高斯记在了自己的小本本上。那么这个算法究竟用了什么奇技淫巧呢，本文将一探究竟。

DFT变换的基础公式为：
X p = ∑ n = 0 N − 1 x n e   − j 2 π N n p , p ∈ { 0 , 1 , … , N − 1 } X_{p}=\sum_{n=0}^{N-1} x_{n} \mathrm{e}^{\:-j \frac{2 \pi}{N} n p}, \quad p \in\{0,1, \ldots, N-1\} \\ Xp​=n=0∑N−1​xn​e−jN2π​np,p∈{0,1,…,N−1}
假设 $N$ 为 $2$ 的正整数次幂，则把上述公式等号右边分成两部分，分别是 $n$ 为奇数时的有限项级数求和以及 $n$ 为偶数时的有限项级数求和：
$$\begin{array}{rl}{X}_p& =\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n}{\mathrm{e}}^{-j\frac{2\pi }{N}(2n)p}+\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n+1}{\mathrm{e}}^{-j\frac{2\pi }{N}(2n+1)p}\\ & =\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n}{e}^{-j\frac{2\pi }{(N/2)}np}+e^{-j\frac{2\pi }{N}p}\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n+1}{e}^{-j\frac{2\pi }{(N/2)}np}\\ & =A_p+W_N^p{B}_p\end{array}$$

其中：
$$\begin{array}{rl}{A}_p& =\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n}{e}^{-j\frac{2\pi }{(N/2)}np}\\ B_p& =\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n+1}{e}^{-j\frac{2\pi }{(N/2)}np}\\ W_N^p& =e^{-j\frac{2\pi }{N}p}\end{array}$$
注意看，这里 $A_p$ 是序列 { x 2 n } = { x 0 , x 2 , … x N − 4 , x N − 2 } \left\{x_{2 n}\right\}=\left\{x_{0}, x_{2}, \ldots x_{N-4}, x_{N-2}\right\} {x2n​}={x0​,x2​,…xN−4​,xN−2​} 的 DFT 变换结果，而 $B_p$ 是序列 { x 2 n + 1 } = { x 1 , x 3 , … x N − 3 , x N − 1 } \left\{x_{2 n+1}\right\}=\left\{x_{1}, x_{3}, \ldots x_{N-3}, x_{N-1}\right\} {x2n+1​}={x1​,x3​,…xN−3​,xN−1​} 的 DFT 结果。因为 $W_N^p$ 可以事先计算并存储，所以如果先计算出 $A_p$ 和 $B_p$ ，再进行 $O(n)$ 次乘法和加法运算，岂不是就可以算出 $X_p$ ? 那么运算量岂不是减少了将近一半？那如果再把 $A_p$ 和 $B_p$ 分别分成左右两部分（假设 $N$ 为 $2$ 的正整数次方），以同样的方式进行计算，如此穷尽下去…等下？先打住、有点过于乐观了…

因为 $X_p$ 中，$p$ 的取值范围是 p ∈ { 0 , 1 , … , N − 1 } p \in\{0,1, \ldots, N-1\} p∈{0,1,…,N−1} , 而 $A_p$ 和 $B_p$ 中，$p$ 的取值范围却是 p ∈ { 0 , 1 , … , N / 2 − 1 } p \in\{0,1, \ldots, N/2-1\} p∈{0,1,…,N/2−1} ，所以在完成 $A_p$ 和 $B_p$ 的计算后，我们就只能算出当 p ∈ { 0 , 1 , … , N / 2 − 1 } p \in\{0,1, \ldots, N/2-1\} p∈{0,1,…,N/2−1} 时的 $X_p$ ,也就是说傅里叶变换只完成了前一半。那运算量减少将近一半的说法不是笑话吗？还真不是，继续往下看。

因为：
$$X_p=\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n}{e}^{-j\frac{2\pi }{(N/2)}np}+e^{-j\frac{2\pi }{N}p}\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n+1}{e}^{-j\frac{2\pi }{(N/2)}np}$$
所以：
$$X_{p+N/2}=\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n}{e}^{-j\frac{2\pi }{(N/2)}n(p+N/2)}+e^{-j\frac{2\pi }{N}(p+N/2)}\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n+1}{e}^{-j\frac{2\pi }{(N/2)}n(p+N/2)}$$
对等号右边的部分表达进行化简，有：
$$\begin{array}{rl}{e}^{-j\frac{2\pi }{(N/2)}n(p+N/2)}& =e^{-j\frac{2\pi }{(N/2)}np}\\ e^{-j\frac{2\pi }{N}(p+N/2)}& =-e^{-j\frac{2\pi }{N}p}\end{array}$$
将化简结果带入上式，有：
$$\begin{array}{rl}{X}_{p+N/2}& =\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n}{e}^{-j\frac{2\pi }{(N/2)}np}-e^{-j\frac{2\pi }{N}p}\sum _{n=0}^{\frac{N}{2}-1}{x}_{2n+1}{e}^{-j\frac{2\pi }{(N/2)}np}\\ & =A_p-W_N^p{B}_p\end{array}$$

所以有：
$$\{\begin{array}{l}{A}_p+W_N^p{B}_p=X_p\\ A_p-W_N^p{B}_p=X_{p+N/2}\end{array}$$
这样，如果先计算出 $A_p$ 和 $B_p$ ，再进行 $O(n)$ 次乘法和加法运算，就可以算出当 p ∈ { 0 , 1 , … , N − 1 } p \in\{0,1, \ldots, N-1\} p∈{0,1,…,N−1} 时的所有 $X_p$ ，这就是著名的蝶形运算，基本运算单元如下图所示(图中的 $W^p$ 和本文的 $W_N^p$ 是一样的 )：

在 $N$ 是一个比较大的数时，进行一次这样的处理就能节约将近一半的计算量，下图给出了当 $N=8$ 时候的第一步的蝶形运算过程：

$A_p$ 是序列 { x 2 n } = { x 0 , x 2 , … x N − 4 , x N − 2 } \left\{x_{2 n}\right\}=\left\{x_{0}, x_{2}, \ldots x_{N-4}, x_{N-2}\right\} {x2n​}={x0​,x2​,…xN−4​,xN−2​} 的 DFT 变换结果， $B_p$ 是序列 { x 2 n + 1 } = { x 1 , x 3 , … x N − 3 , x N − 1 } \left\{x_{2 n+1}\right\}=\left\{x_{1}, x_{3}, \ldots x_{N-3}, x_{N-1}\right\} {x2n+1​}={x1​,x3​,…xN−3​,xN−1​} 的 DFT 变换结果。那么如果再分别对 $A_p$ 和 $B_p$ 采用上面的方法进行计算，当序列长度较大时，则又能节约将近一半的计算量！即$A_p$ 和 $B_p$ 采用下面的公式计算：
$$\begin{array}{r}{A}_p={\alpha }_p+W_{N/2}^p{\beta }_p{A}_{p+N/2}={\alpha }_p-W_{N/2}^p{\beta }_p\\ B_p={\alpha }_p^{\prime }+W_{N/2}^p{\beta }_p^{\prime }{B}_{p+N/2}={\alpha }_p^{\prime }-W_{N/2}^p{\beta }_p^{\prime }\end{array}$$
式中：
$$W_{N/2}^p=e^{-j\frac{2\pi }{N/2}p}=W_N^{2p}$$
下图给出了当 $N=8$ 时候的进一步的蝶形运算过程：

如果递归式地进行这样的处理，就能把算法复杂度从 $O(n^2)$ 降到 $O(n\log n)$ .

本文首发于微信公众号：“振动信号研究所”