离散傅里叶变换（DFT）是为了能让计算机处理经过采集的离散信号提出的，DFT的推导以及代码实现在下面的两篇文章中也已经介绍了：
DFT推导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−1xne−jN2πnp,p∈{0,1,…,N−1} 为了更清楚地计算DFT的算法复杂度，我们将其表示为矩阵表达：
[
X
0
X
1
⋯
X
N
−
1
]
=
[
W
N
0
W
N
0
W
N
0
⋯
W
N
0
W
N
0
W
N
1
W
N
2
⋯
W
N
N
−
1
⋮
⋮
⋮
⋯
⋮
W
N
0
W
N
N
−
1
W
N
2
(
N
−
1
)
⋯
W
N
(
N
−
1
)
2
]
[
x
0
x
1
⋮
x
N
−
1
]
\begin{bmatrix} X_0 \\ X_1 \\\cdots \\X_{N-1} \end{bmatrix} =\begin{bmatrix} 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}\\ \vdots & \vdots & \vdots & \cdots &\vdots\\ W_N^0 & W_N^{N-1} & W_N^{2(N-1)} & \cdots & W_N^{(N-1)^2} \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\\vdots \\x_{N-1} \end{bmatrix}
X0X1⋯XN−1
=
WN0WN0⋮WN0WN0WN1⋮WNN−1WN0WN2⋮WN2(N−1)⋯⋯⋯⋯WN0WNN−1⋮WN(N−1)2
x0x1⋮xN−1
式中：
W
N
=
e
−
j
2
π
N
W_N=e^{\:-j\frac{2\pi}{N}}
WN=e−jN2π
2. DFT复杂度计算
因为
N
∗
N
N*N
N∗N 的变换矩阵中的每一个元素可以离线计算并存储，这样上述计算过程仅涉及
N
2
N^2
N2 次复数乘法运算和
N
(
N
−
1
)
N(N-1)
N(N−1) 次复数加法运算。
因为变换矩阵的第一行和第一列都是
1
1
1 ，这样就节约了
2
N
−
1
2N - 1
2N−1 次的复数乘法运算，所以实际上，上述计算过程只需要进行
N
2
−
(
2
N
−
1
)
N^2-(2N-1)
N2−(2N−1) 次复数乘法运算，即共需要
(
N
−
1
)
2
(N-1)^2
(N−1)2 次复数乘法运算和
N
(
N
−
1
)
N(N-1)
N(N−1) 次复数加法运算。
因为
1
1
1 次复数乘法运算相当于
4
4
4 次实数乘法运算和
2
2
2 次实数加法运算，
1
1
1 次复数加法运算相当于
2
2
2 次实数加法运算，所以：
一共需要进行
4
(
N
−
1
)
2
4(N-1)^2
4(N−1)2 次实数乘法运算和
2
(
N
−
1
)
2
+
2
N
(
N
−
1
)
=
4
(
N
−
0.5
)
(
N
−
1
)
2(N-1)^2+2N(N-1)=4(N-0.5)(N-1)
2(N−1)2+2N(N−1)=4(N−0.5)(N−1) 次实数加法运算。
所以离散傅里叶变换算法复杂度是
O
(
n
2
)
O(n^2)
O(n2) ，如果要处理的数据点数过多，就会需要大量的计算时间，这会给工程应用带来很大的限制。
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−1xne−jN2πnp,p∈{0,1,…,N−1} 假设
N
N
N 为
2
2
2 的正整数次幂，则把上述公式等号右边分成两部分，分别是
n
n
n 为奇数时的有限项级数求和以及
n
n
n 为偶数时的有限项级数求和：
X
p
=
∑
n
=
0
N
2
−
1
x
2
n
e
−
j
2
π
N
(
2
n
)
p
+
∑
n
=
0
N
2
−
1
x
2
n
+
1
e
−
j
2
π
N
(
2
n
+
1
)
p
=
∑
n
=
0
N
2
−
1
x
2
n
e
−
j
2
π
(
N
/
2
)
n
p
+
e
−
j
2
π
N
p
∑
n
=
0
N
2
−
1
x
2
n
+
1
e
−
j
2
π
(
N
/
2
)
n
p
=
A
p
+
W
N
p
B
p
\begin{aligned} X_{p}&=\sum_{n=0}^{\frac{N}{2}-1} x_{2 n} \mathrm{e}^{-j \frac{2 \pi}{N}(2 n) p}+\sum_{n=0}^{\frac{N}{2}-1} x_{2 n+1} \mathrm{e}^{-j \frac{2 \pi}{N}(2 n+1) p}\\ &=\sum_{n=0}^{\frac{N}{2}-1} x_{2 n} e^{-j \frac{2 \pi}{(N / 2)} n p}+e^{-j \frac{2 \pi}{N} p} \sum_{n=0}^{\frac{N}{2}-1} x_{2 n+1} e^{-j \frac{2 \pi}{(N / 2)} n p}\\ &= A_p + W_N^pB_p \end{aligned}
Xp=n=0∑2N−1x2ne−jN2π(2n)p+n=0∑2N−1x2n+1e−jN2π(2n+1)p=n=0∑2N−1x2ne−j(N/2)2πnp+e−jN2πpn=0∑2N−1x2n+1e−j(N/2)2πnp=Ap+WNpBp
其中：
A
p
=
∑
n
=
0
N
2
−
1
x
2
n
e
−
j
2
π
(
N
/
2
)
n
p
B
p
=
∑
n
=
0
N
2
−
1
x
2
n
+
1
e
−
j
2
π
(
N
/
2
)
n
p
W
N
p
=
e
−
j
2
π
N
p
\begin{aligned} A_p&=\sum_{n=0}^{\frac{N}{2}-1} x_{2 n} e^{-j \frac{2 \pi}{(N / 2)} n p}\\ B_p&=\sum_{n=0}^{\frac{N}{2}-1} x_{2 n+1} e^{-j \frac{2 \pi}{(N / 2)} n p}\\ W_N^p&= e^{-j \frac{2 \pi}{N} p} \end{aligned}
ApBpWNp=n=0∑2N−1x2ne−j(N/2)2πnp=n=0∑2N−1x2n+1e−j(N/2)2πnp=e−jN2πp 注意看，这里
A
p
A_p
Ap 是序列
{
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
B_p
Bp 是序列
{
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
W_N^p
WNp 可以事先计算并存储，所以如果先计算出
A
p
A_p
Ap 和
B
p
B_p
Bp ，再进行
O
(
n
)
O(n)
O(n) 次乘法和加法运算，岂不是就可以算出
X
p
X_p
Xp ? 那么运算量岂不是减少了将近一半？那如果再把
A
p
A_p
Ap 和
B
p
B_p
Bp 分别分成左右两部分（假设
N
N
N 为
2
2
2 的正整数次方），以同样的方式进行计算，如此穷尽下去…等下？先打住、有点过于乐观了…
因为
X
p
X_p
Xp 中，
p
p
p 的取值范围是
p
∈
{
0
,
1
,
…
,
N
−
1
}
p \in\{0,1, \ldots, N-1\}
p∈{0,1,…,N−1} , 而
A
p
A_p
Ap 和
B
p
B_p
Bp 中，
p
p
p 的取值范围却是
p
∈
{
0
,
1
,
…
,
N
/
2
−
1
}
p \in\{0,1, \ldots, N/2-1\}
p∈{0,1,…,N/2−1} ，所以在完成
A
p
A_p
Ap 和
B
p
B_p
Bp 的计算后，我们就只能算出当
p
∈
{
0
,
1
,
…
,
N
/
2
−
1
}
p \in\{0,1, \ldots, N/2-1\}
p∈{0,1,…,N/2−1} 时的
X
p
X_p
Xp ,也就是说傅里叶变换只完成了前一半。那运算量减少将近一半的说法不是笑话吗？还真不是，继续往下看。
因为：
X
p
=
∑
n
=
0
N
2
−
1
x
2
n
e
−
j
2
π
(
N
/
2
)
n
p
+
e
−
j
2
π
N
p
∑
n
=
0
N
2
−
1
x
2
n
+
1
e
−
j
2
π
(
N
/
2
)
n
p
X_{p}=\sum_{n=0}^{\frac{N}{2}-1} x_{2 n} e^{-j \frac{2 \pi}{(N / 2)} n p}+e^{-j \frac{2 \pi}{N} p} \sum_{n=0}^{\frac{N}{2}-1} x_{2 n+1} e^{-j \frac{2 \pi}{(N / 2)} n p}
Xp=n=0∑2N−1x2ne−j(N/2)2πnp+e−jN2πpn=0∑2N−1x2n+1e−j(N/2)2πnp 所以：
X
p
+
N
/
2
=
∑
n
=
0
N
2
−
1
x
2
n
e
−
j
2
π
(
N
/
2
)
n
(
p
+
N
/
2
)
+
e
−
j
2
π
N
(
p
+
N
/
2
)
∑
n
=
0
N
2
−
1
x
2
n
+
1
e
−
j
2
π
(
N
/
2
)
n
(
p
+
N
/
2
)
X_{p+N/2}=\sum_{n=0}^{\frac{N}{2}-1} x_{2 n} 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_{2 n+1} e^{-j \frac{2 \pi}{(N / 2)} n (p+N/2)}
Xp+N/2=n=0∑2N−1x2ne−j(N/2)2πn(p+N/2)+e−jN2π(p+N/2)n=0∑2N−1x2n+1e−j(N/2)2πn(p+N/2) 对等号右边的部分表达进行化简，有：
e
−
j
2
π
(
N
/
2
)
n
(
p
+
N
/
2
)
=
e
−
j
2
π
(
N
/
2
)
n
p
e
−
j
2
π
N
(
p
+
N
/
2
)
=
−
e
−
j
2
π
N
p
\begin{aligned} e^{-j \frac{2 \pi}{(N / 2)} n(p+N / 2)}&=e^{-j \frac{2 \pi}{(N / 2)} n p}\\ e^{-j \frac{2 \pi}{N} (p+N/2)} &= -e^{-j \frac{2 \pi}{N}p} \end{aligned}
e−j(N/2)2πn(p+N/2)e−jN2π(p+N/2)=e−j(N/2)2πnp=−e−jN2πp 将化简结果带入上式，有：
X
p
+
N
/
2
=
∑
n
=
0
N
2
−
1
x
2
n
e
−
j
2
π
(
N
/
2
)
n
p
−
e
−
j
2
π
N
p
∑
n
=
0
N
2
−
1
x
2
n
+
1
e
−
j
2
π
(
N
/
2
)
n
p
=
A
p
−
W
N
p
B
p
\begin{aligned} X_{p+N/2}&=\sum_{n=0}^{\frac{N}{2}-1} x_{2 n} e^{-j \frac{2 \pi}{(N / 2)} n p}-e^{-j \frac{2 \pi}{N}p} \sum_{n=0}^{\frac{N}{2}-1} x_{2 n+1} e^{-j \frac{2 \pi}{(N / 2)} n p}\\ &=A_p - W_N^pB_p \end{aligned}
Xp+N/2=n=0∑2N−1x2ne−j(N/2)2πnp−e−jN2πpn=0∑2N−1x2n+1e−j(N/2)2πnp=Ap−WNpBp
所以有：
{
A
p
+
W
N
p
B
p
=
X
p
A
p
−
W
N
p
B
p
=
X
p
+
N
/
2
\begin{cases} A_p + W_N^pB_p = X_p\\ A_p - W_N^pB_p =X_{p+N/2} \end{cases}
{Ap+WNpBp=XpAp−WNpBp=Xp+N/2 这样，如果先计算出
A
p
A_p
Ap 和
B
p
B_p
Bp ，再进行
O
(
n
)
O(n)
O(n) 次乘法和加法运算，就可以算出当
p
∈
{
0
,
1
,
…
,
N
−
1
}
p \in\{0,1, \ldots, N-1\}
p∈{0,1,…,N−1} 时的所有
X
p
X_p
Xp ，这就是著名的蝶形运算，基本运算单元如下图所示(图中的
W
p
W^p
Wp 和本文的
W
N
p
W_N^p
WNp 是一样的 )：
在
N
N
N 是一个比较大的数时，进行一次这样的处理就能节约将近一半的计算量，下图给出了当
N
=
8
N=8
N=8 时候的第一步的蝶形运算过程：
A
p
A_p
Ap 是序列
{
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
B_p
Bp 是序列
{
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
A_p
Ap 和
B
p
B_p
Bp 采用上面的方法进行计算，当序列长度较大时，则又能节约将近一半的计算量！即
A
p
A_p
Ap 和
B
p
B_p
Bp 采用下面的公式计算：
A
p
=
α
p
+
W
N
/
2
p
β
p
A
p
+
N
/
2
=
α
p
−
W
N
/
2
p
β
p
B
p
=
α
p
′
+
W
N
/
2
p
β
p
′
B
p
+
N
/
2
=
α
p
′
−
W
N
/
2
p
β
p
′
\begin{aligned} A_{p}=\alpha_{p}+W_{N / 2}^{p} \beta_{p} \quad A_{p+N / 2}=\alpha_{p}-W_{N / 2}^{p} \beta_{p}\\ B_{p}=\alpha_{p}^{\prime}+W_{N / 2}^p \beta_{p}^{\prime} \quad B_{p+N / 2}=\alpha_{p}^{\prime}-W_{N / 2}^p \beta_{p}^{\prime} \end{aligned}
Ap=αp+WN/2pβpAp+N/2=αp−WN/2pβpBp=αp′+WN/2pβp′Bp+N/2=αp′−WN/2pβp′ 式中：
W
N
/
2
p
=
e
−
j
2
π
N
/
2
p
=
W
N
2
p
W_{N/2}^p = e^{-j \frac{2 \pi}{N/2} p}=W_{N}^{2p}
WN/2p=e−jN/22πp=WN2p 下图给出了当
N
=
8
N=8
N=8 时候的进一步的蝶形运算过程：
如果递归式地进行这样的处理，就能把算法复杂度从
O
(
n
2
)
O(n^2)
O(n2) 降到
O
(
n
log
n
)
O(n\log n)
O(nlogn) .
本文首发于微信公众号：“振动信号研究所”