傅里叶变换：FFT 与 DFT 原理及算法

图 TC.3.1 按时间抽取算法的第一步。

通过计算 $N/4$ 点 DFT，我们可以从以下关系式得到 $N/2$ 点 DFT $F_1(k)$ 和 $F_2(k)$

$$\begin{align*} F_1(k) &= \mathrm{F}{f_1(2n)} + W_{N/2}^k \mathrm{F}{f_1(2n+1)}, \ F_1!\left(k+\frac{N}{4}\right) &= \mathrm{F}{f_1(2n)} - W_{N/2}^k \mathrm{F}{f_1(2n+1)}, \ F_2(k) &= \mathrm{F}{f_2(2n)} + W_{N/2}^k \mathrm{F}{f_2(2n+1)}, \ F_2!\left(k+\frac{N}{4}\right) &= \mathrm{F}{f_2(2n)} - W_{N/2}^k \mathrm{F}{f_2(2n+1)} \end{align*} \quad (k = 0,1,\dots,\frac{N}{4}-1;\ n = 0,1,\dots,\frac{N}{4}-1)$$

数据序列的抽取可以反复进行，直到结果序列缩减为单点序列。对于 $N=2^v$，这种抽取可以进行 $v = \log_2 N$ 次。因此，复数乘法的总次数减少到 $(N/2)\log_2 N$。复数加法的次数为 $N\log_2 N$。

为便于说明，图 TC.3.2 描绘了 $N=8$ 点 DFT 的计算。我们观察到计算分三个阶段进行：首先计算四个两点 DFT，然后计算两个四点 DFT，最后计算一个八点 DFT。图 TC.3.3 说明了 $N=8$ 时较小 DFT 组合形成较大 DFT 的过程。

图 TC.3.2 计算 $N=8$ 点 DFT 的三个阶段。

图 TC.3.3 八点按时间抽取 FFT 算法。

图 TC.3.4 按时间抽取 FFT 算法中的基本蝶形运算。

一个重要的观察点是输入数据序列经过 $(v-1)$ 次抽取后的顺序。例如，若考虑 $N=8$ 的情况，我们知道第一次抽取产生序列 $x(0), x(2), x(4), x(6), x(1), x(3), x(5), x(7)$，第二次抽取产生序列 $x(0), x(4), x(2), x(6), x(1), x(5), x(3), x(7)$。这种输入数据序列的混洗（shuffling）具有明确的顺序，可从图 TC.3.5 中看出，该图说明了八点序列的抽取过程。

图 TC.3.5 数据的混洗与位反转。

另一种重要的基 -2 FFT 算法称为按频率抽取算法（decimation-in-frequency algorithm），通过分治法得到。为推导该算法，我们首先将 DFT 公式分成两个求和式，其中一个涉及前 $N/2$ 个数据点的求和，第二个涉及后 $N/2$ 个数据点的求和。因此我们得到

$$\begin{align*} X(k) &= \sum_{n=0}^{N/2-1} x(n) W_N^{nk} + \sum_{n=N/2}^{N-1} x(n) W_N^{nk} \ &= \sum_{n=0}^{N/2-1} \left[ x(n) + x!\left(n+\frac{N}{2}\right) W_N^{Nk/2} \right] W_N^{nk} \end{align*}$$

现在，让我们将 $X(k)$ 分成偶数编号和奇数编号的样本。因此我们得到

$$\begin{align*} X(2k) &= \sum_{n=0}^{N/2-1} \left[ x(n) + x!\left(n+\frac{N}{2}\right) \right] W_{N/2}^{,nk}, && k = 0, 1, \ldots, \frac{N}{2}-1, \ X(2k+1) &= \sum_{n=0}^{N/2-1} \left[ x(n) - x!\left(n+\frac{N}{2}\right) \right] W_N^{,n} W_{N/2}^{,nk}, && k = 0, 1, \ldots, \frac{N}{2}-1. \end{align*}$$

其中我们利用了 $W_N^2 = W_{N/2}$ 这一事实。

上述计算过程可以通过对 $N/2$ 点 DFT $X(2k)$ 和 $X(2k+1)$ 的抽取来重复。整个过程涉及 $v = \log_2 N$ 个抽取阶段，每个阶段涉及 $N/2$ 个如图 TC.3.7 所示类型的蝶形运算。因此，通过按频率抽取 FFT 计算 $N$ 点 DFT 需要 $(N/2)\log_2 N$ 次复数乘法和 $N\log_2 N$ 次复数加法，与按时间抽取算法相同。为便于说明，八点按频率抽取算法示于图 TC.3.8。

图 TC.3.6 按频率抽取 FFT 算法的第一阶段。

图 TC.3.7 按频率抽取中的基本蝶形运算。

图 TC.3.8 $N=8$ 点按频率抽取 FFT 算法。

从图 TC.3.8 我们观察到，输入数据 $x(n)$ 按自然顺序出现，但输出 DFT 按位反转顺序出现。我们还注意到计算是原位进行的（in place）。然而，可以重新配置按频率抽取算法，使输入序列按位反转顺序出现，而输出 DFT 按正常顺序出现。此外，如果我们放弃计算必须原位进行的要求，也可以使输入数据和输出 DFT 都按正常顺序出现。

基 -4 FFT 算法

当 DFT 中的数据点数 $N$ 是 4 的幂（即 $N=4^v$）时，我们当然总是可以使用基 -2 算法进行计算。然而，在这种情况下，采用基 -4 FFT 算法在计算上更为高效。

让我们首先简要描述基 -4 按时间抽取 FFT 算法。我们将 $N$ 点输入序列分成四个子序列：$x(4n), x(4n+1), x(4n+2), x(4n+3)$，其中 $n = 0, 1, \ldots, N/4-1$。

$$\begin{align*} X(k) &= \sum_{n=0}^{N/4-1} x(4n) W_N^{,4nk} + \sum_{n=0}^{N/4-1} x(4n+1) W_N^{,(4n+1)k} \ &\quad + \sum_{n=0}^{N/4-1} x(4n+2) W_N^{,(4n+2)k} + \sum_{n=0}^{N/4-1} x(4n+3) W_N^{,(4n+3)k} \ &= \sum_{l=0}^{3} W_N^{,lk} \sum_{n=0}^{N/4-1} x(4n+l) W_{N/4}^{,nk}. \end{align*}$$

因此，从上述方程得到的四个 $N/4$ 点 DFT $F(l,q)$ 被组合起来产生 $N$ 点 DFT。组合 $N/4$ 点 DFT 的表达式定义了一个基 -4 按时间抽取蝶形运算，可以用矩阵形式表示为

$$\begin {bmatrix} X(k) \ X(k+N/4) \ X(k+N/2) \ X(k+3N/4) \end {bmatrix} = \begin {bmatrix} 1 & 1 & 1 & 1 \ 1 & -j & -1 & j \ 1 & -1 & 1 & -1 \ 1 & j & -1 & -j \end {bmatrix} \begin {bmatrix} W_N^0 F(0,k) \ W_N^k F(1,k) \ W_N^{2k} F(2,k) \ W_N^{3k} F(3,k) \end {bmatrix}$$

基 -4 蝶形运算示于图 TC.3.9a，更紧凑的形式示于图 TC.3.9b。注意，每个蝶形运算涉及三次复数乘法（因为 $W_N^0=1$）和 12 次复数加法。

图 TC.3.9 基 -4 FFT 算法中的基本蝶形运算。

通过分两步执行加法，可以将每个蝶形运算的加法次数从 12 减少到 8。这可以通过将前述线性变换的矩阵表示为两个矩阵的乘积来实现，如下所示：

$$\begin{bmatrix} X(0,q) \ X(1,q) \ X(2,q) \ X(3,q) \end{bmatrix} = \begin{bmatrix} 1 & 0 & 1 & 0 \ 0 & 1 & 0 & -j \ 1 & 0 & -1 & 0 \ 0 & 1 & 0 & j \end{bmatrix} \begin{bmatrix} 1 & 0 & 1 & 0 \ 1 & 0 & -1 & 0 \ 0 & 1 & 0 & 1 \ 0 & 1 & 0 & -1 \end{bmatrix} \begin{bmatrix} W_N^0 F(0,q) \ W_N^1 F(1,q) \ W_N^2 F(2,q) \ W_N^3 F(3,q) \end{bmatrix}$$

图 TC.3.10 输入为正常顺序、输出为数字反转顺序的 16 点基 -4 按时间抽取算法

图 TC.3.11 示出了一种 16 点基 -4 按频率抽取 FFT 算法。其输入为正常顺序，输出为数字反转顺序。其计算复杂度与基 -4 按时间抽取 FFT 算法完全相同。

图 TC.3.11 输入为正常顺序、输出为数字反转顺序的 16 点基 -4 按频率抽取算法。

为便于说明，让我们重新推导基 4 频率抽取算法，通过将 $N$ 点 DFT 公式分解为四个较小的 DFT。我们有

$$\begin {aligned} X(k) &= \sum_{n=0}^{N-1} x(n) W_N^{kn} \ &= \sum_{n=0}^{N/4-1} x(n) W_N^{kn} + \sum_{n=N/4}^{N/2-1} x(n) W_N^{kn} + \sum_{n=N/2}^{3N/4-1} x(n) W_N^{kn} + \sum_{n=3N/4}^{N-1} x(n) W_N^{kn} \ &= \sum_{n=0}^{N/4-1} x(n) W_N^{kn} + W_N^{Nk/4} \sum_{n=0}^{N/4-1} x\left (n+\frac {N}{4}\right) W_N^{kn} + \ &\quad W_N^{Nk/2} \sum_{n=0}^{N/4-1} x\left (n+\frac {N}{2}\right) W_N^{kn} + W_N^{3Nk/4} \sum_{n=0}^{N/4-1} x\left (n+\frac {3N}{4}\right) W_N^{kn} \end {aligned}$$

根据旋转因子的定义，我们有 $W_N^{kN/4} = (-j)^k, \quad W_N^{kN/2} = (-1)^k, \quad W_N^{3kN/4} = (j)^k$。

因此

$$X(k) = \sum_{n=0}^{N/4-1} \left [ x(n) + (-j)^k x\left (n+\frac {N}{4}\right) + (-1)^k x\left (n+\frac {N}{2}\right) + (j)^k x\left (n+\frac {3N}{4}\right) \right] W_N^{nk}$$

该关系式不是 $N/4$ 点 DFT，因为旋转因子依赖于 $N$ 而非 $N/4$。为将其转换为 $N/4$ 点 DFT，我们将 DFT 序列分成四个 $N/4$ 点子序列：$X(4k), X(4k+1), X(4k+2)$ 和 $X(4k+3)$，其中 $k = 0, 1, \ldots, N/4$。因此我们得到基 -4 按频率抽取 DFT 为

$$\begin{aligned} X(4k) &= \sum_{n=0}^{N/4-1} \left[ x(n) + x\left(n+\frac{N}{4}\right) + x\left(n+\frac{N}{2}\right) + x\left(n+\frac{3N}{4}\right) \right] W_{N/4}^{nk} \ X(4k+1) &= \sum_{n=0}^{N/4-1} \left[ x(n) - jx\left(n+\frac{N}{4}\right) - x\left(n+\frac{N}{2}\right) + jx\left(n+\frac{3N}{4}\right) \right] W_N^n W_{N/4}^{nk} \ X(4k+2) &= \sum_{n=0}^{N/4-1} \left[ x(n) - x\left(n+\frac{N}{4}\right) + x\left(n+\frac{N}{2}\right) - x\left(n+\frac{3N}{4}\right) \right] W_N^{2n} W_{N/4}^{nk} \ X(4k+3) &= \sum_{n=0}^{N/4-1} \left[ x(n) + jx\left(n+\frac{N}{4}\right) - x\left(n+\frac{N}{2}\right) - jx\left(n+\frac{3N}{4}\right) \right] W_N^{3n} W_{N/4}^{nk} \end{aligned}$$

其中我们利用了性质 $W_N^{4kn} = W_{N/4}^{kn}$。注意，每个 $N/4$ 点 DFT 的输入是四个信号样本经旋转因子加权后的线性组合。此过程重复 $v$ 次，其中 $v = \log_4 N$。

分裂基 FFT 算法

观察图 TC.3.8 所示的基 -2 按频率抽取流程图可知，DFT 的偶数编号点可以独立于奇数编号点进行计算。这提示我们可以在算法的独立部分使用不同的计算方法，以减少计算次数。分裂基 FFT（SRFFT）算法利用这一思想，在同一个 FFT 算法中同时使用基 -2 和基 -4 分解。

首先，我们回顾在基 -2 按频率抽取 FFT 算法中，$N$ 点 DFT 的偶数编号样本由下式给出

$$X(2k) = \sum_{n=0}^{N/2-1} \left [x(n) + x\left (n+\frac {N}{2}\right)\right] W_{N/2}^{nk}, \quad k = 0, 1, \ldots, \frac {N}{2}-1$$

此计算使用基 -2 已足够。

DFT 的奇数编号样本 $ ext{ extbraceleft} X(2k+1) ext{ extbraceright}$ 需要先将输入序列与旋转因子 $W_N^n$ 相乘。对于这些样本，基 -4 分解能产生一定的计算效率，因为四点 DFT 具有最大的无乘法蝶形运算。事实上，可以证明使用大于 4 的基不会显著降低计算复杂度。

如果我们对 $N$ 点 DFT 的奇数编号样本使用基 -4 按频率抽取 FFT 算法，我们得到以下 $N/4$ 点 DFT：

$$\begin{aligned} X(4k+1) &= \sum_{n=0}^{N/4-1} \left{\left [x(n) - x\left (n+\frac {N}{2}\right)\right] - j\left [x\left (n+\frac {N}{4}\right) - x\left (n+\frac {3N}{4}\right)\right]\right} W_N^n W_{N/4}^{nk} \ X(4k+3) &= \sum_{n=0}^{N/4-1} \left{\left [x(n) - x\left (n+\frac {N}{2}\right)\right] + j\left [x\left (n+\frac {N}{4}\right) - x\left (n+\frac {3N}{4}\right)\right]\right} W_N^{3n} W_{N/4}^{nk} \end{aligned}$$

图 TC.3.12 示出了原位 32 点按频率抽取 SFFT 算法的流程图。

图 TC.3.12 Duhamel（1986）论文中的长度 32 分裂基 FFT 算法；经 IEEE 许可转载

图 TC.3.13 SRFFT 算法的蝶形运算。

表 TC.3.1 计算 $N$ 点复数 DFT 所需的非平凡实数乘法和加法次数

DFT 定义与实例分析

1. DFT 的定义

连续时间傅里叶变换的定义为：

$$X(f) = \int_{-\infty}^{+\infty}x(t), e^{-j2\pi f t},dt$$

离散傅里叶变换（DFT）由连续傅里叶变换离散化得到：

$$X(m) = \sum_{n=0}^{N-1} x(n), e^{-j\frac{2\pi m n}{N}}$$

式中：

• $N$：时域离散信号的点数
• $n$：时域离散信号序号，满足 $n \in {0, 1, \dots, N-1}$
• $m$：频域信号序号，满足 $m \in {0, 1, \dots, N-1}$

DFT 输入为 $N$ 个时域离散样本，输出为 $N$ 个复数值频域样本。

由欧拉公式：$e^{-j\omega} = \cos\omega - j\sin\omega$。

DFT 可展开为实部虚部形式：

$$X(m) = \sum_{n=0}^{N-1} x(n)\left[ \cos\left(\frac{2\pi m n}{N}\right)- j\sin\left(\frac{2\pi m n}{N}\right) \right]$$

2. DFT 的实例分析

设时域连续信号为：

$$x(t) = \sin(2\pi\cdot 1000\cdot t)+ 0.5\sin\left(2\pi\cdot 2000\cdot t + \frac{3\pi}{4}\right)$$

以采样频率 $f_s=8000\ \text{Hz}$ 采样，采样周期 $T_s=1/8000\ \text{s}$，得到离散信号：

$$\begin{aligned} x(n) &= \sin(2\pi\cdot 1000\cdot n T_s)+ 0.5\sin\left(2\pi\cdot 2000\cdot n T_s + \frac{3\pi}{4}\right) \ &= \sin\left(\frac{2\pi n}{8}\right)+ 0.5\sin\left(\frac{4\pi n}{8}+\frac{3\pi}{4}\right) \end{aligned}$$

取 $N=8$，即 $n = 0, 1, \dots, 7$，得到离散样本：

$$\begin{aligned} &x(0)= 0.3535,\quad x(1)= 0.3535,\quad x(2)= 0.6464,\quad x(3)= 1.0607,\ &x(4)= 0.3535,\quad x(5)=-1.0607,\quad x(6)=-1.3535,\quad x(7)=-0.3535 \end{aligned}$$

Fig.1 输入信号

视角一：以 $n$ 为自变量计算 DFT

对固定序号 $m=i$：

• $X(i)$ 的实部：$x(n)$ 与 $\cos\left(\frac{2\pi i n}{8}\right)$ 做点积
• $X(i)$ 的虚部：$x(n)$ 与 $-\sin\left(\frac{2\pi i n}{8}\right)$ 做点积

本例 DFT 计算结果：

$$\begin{aligned} &X(0)= 0.0-j0.0,\quad X(1)= 0.0-j4.0,\quad X(2)= 1.414+j1.414,\quad X(3)= 0.0+j0.0,\ &X(4)= 0.0-j0.0,\quad X(5)= 0.0-j0.0,\quad X(6)= 1.414-j1.414,\quad X(7)= 0.0+j4.0 \end{aligned}$$

Fig. 2 $X(0)$ 的计算方法

Fig. 3 $X(1)$ 的计算方法

Fig. 4 $X(2)$ 的计算方法

Fig. 5 $X(3)$ 的计算方法

Fig. 6 $X(4)$ 的计算方法

Fig. 7 $X(5)$ 的计算方法

Fig. 8 $X(6)$ 的计算方法

Fig. 9 $X(7)$ 的计算方法

视角二：以频率索引 $k$ 为自变量观察 DFT

在频域视角下，离散傅里叶变换（DFT）从频率维度拆解时域信号的构成，频率索引 $k$ 作为自变量可揭示频域信号的构成规律：对每个 $k$ 计算 DFT，本质是将时域信号 $x(n)$ 与复指数基函数 $e^{-j2\pi kn/N}$ 做内积，即对 $N$ 个等距频率分量的加权求和。

1. 数学推导基础

依据 DFT 的标准定义：

$$X(k)=\sum_{n=0}^{N-1} x(n),e^{-j\frac{2\pi kn}{N}}\qquad k=0,1,\dots ,N-1$$

利用欧拉公式 $e^{-j\theta}= \cos\theta - j\sin\theta$ 将复指数展开为正、负正弦形式，可得：

$$\begin{aligned} X(k) &= \sum_{n=0}^{N-1} x(n)\bigl[\cos!\bigl(\tfrac{2\pi kn}{N}\bigr)-j\sin!\bigl(\tfrac{2\pi kn}{N}\bigr)\bigr] \ &= \underbrace{\sum_{n=0}^{N-1} x(n)\cos!\bigl(\tfrac{2\pi kn}{N}\bigr)}{\displaystyle \operatorname{Re}{X(k)}} ;-;j\underbrace{\sum{n=0}^{N-1} x(n)\sin!\bigl(\tfrac{2\pi kn}{N}\bigr)}_{\displaystyle \operatorname{Im}{X(k)}} . \end{aligned}$$

2. 实部与虚部的规范表述

3. 要素说明

每个频域点 $X(k)$ 由 $N$ 组不同频率的正余弦信号叠加构造，关键要素如下：

注意：$k$ 与 $n$ 取值范围相同，但物理意义截然不同；$k$ 固定、$n$ 遍历时，求得单一频域点 $X(k)$ 的值；$k$ 取遍 $0, 1, \dots ,N-1$ 时，可得到完整的频谱 $ ext{ extbraceleft} X(0), X(1), \dots ,X(N-1) ext{ extbraceright}$。

4. 完整 DFT 表达式

$$X(k)=\sum_{n=0}^{N-1} x(n),e^{-j\frac{2\pi kn}{N}} =\sum_{n=0}^{N-1} x(n)\Bigl[\cos!\bigl(\tfrac{2\pi kn}{N}\bigr)-j\sin!\bigl(\tfrac{2\pi kn}{N}\bigr)\Bigr], \qquad k=0,1,\dots ,N-1$$

该形式清晰展示了频率索引 $k$ 作为自变量、时域采样 $x(n)$ 与基函数相位 $\frac{2\pi kn}{N}$ 交互产生的加权求和过程，便于在频域分析中系统刻画信号的频率分量。

DFT 结果与视角一一致：

$$\begin{aligned} &X(0)= 0.0-j0.0,\quad X(1)= 0.0-j4.0,\quad X(2)= 1.414+j1.414,\quad X(3)= 0.0+j0.0,\ &X(4)= 0.0-j0.0,\quad X(5)= 0.0-j0.0,\quad X(6)= 1.414-j1.414,\quad X(7)= 0.0+j4.0 \end{aligned}$$

Fig. 10 DFT 计算的第二视角

3. DFT 结果

本例 DFT 输出：

$$\begin{aligned} &X(0)= 0.0-j0.0,\quad X(1)= 0.0-j4.0,\quad X(2)= 1.414+j1.414,\quad X(3)= 0.0+j0.0,\ &X(4)= 0.0-j0.0,\quad X(5)= 0.0-j0.0,\quad X(6)= 1.414-j1.414,\quad X(7)= 0.0+j4.0 \end{aligned}$$

分别绘制实部、虚部、幅值、相位：

Fig. 11 DFT 计算结果

4. DFT 的共轭对称性

由计算结果可观察到：

$$\begin{aligned} |X(5)| &= |X(3)|, \quad & |X(6)| &= |X(2)|, \quad & |X(7)| &= |X(1)| \ \operatorname{Re}{X(5)} &= \operatorname{Re}{X(3)}, \quad & \operatorname{Re}{X(6)} &= \operatorname{Re}{X(2)}, \quad & \operatorname{Re}{X(7)} &= \operatorname{Re}{X(1)} \ \operatorname{Im}{X(5)} &= -\operatorname{Im}{X(3)}, \quad & \operatorname{Im}{X(6)} &= -\operatorname{Im}{X(2)}, \quad & \operatorname{Im}{X(7)} &= -\operatorname{Im}{X(1)} \end{aligned}$$

上述关系对应 DFT 的共轭对称性：

$$X(m) = X^*(N-m)$$

式中 $(\cdot)^*$ 表示复共轭。

受共轭对称性约束，DFT 输出中后 $N/2-1$ 个样本存在信息冗余，有效输出可视为 $N/2+1$ 个独立复数。

5. 序号 $m$ 的物理含义

在 $N$ 点 DFT 中，序号 $m$ 对应信号在 $N$ 个点内包含 $m$ 个周期。

对采样频率 $f_s$，第 $m$ 根谱线对应的频率为：

$$f_m = m\cdot\frac{f_s}{N}$$

本例中：

• $m=1 \Rightarrow f_1 = 1000\ \text{Hz}$
• $m=2 \Rightarrow f_2 = 2000\ \text{Hz}$

6. DFT 幅值的含义

原始信号中：

• 1000 Hz 分量幅值为 1
• 2000 Hz 分量幅值为 0.5

而 DFT 幅值分别为 4 和 2。原因在于：傅里叶变换输出为频谱密度，与实际频谱之间相差系数 $N/2$。本例 $N=8$，对应比例为 4。

7. IDFT 的定义

离散傅里叶反变换（IDFT）定义为：

$$x(n) = \frac{1}{N}\sum_{m=0}^{N-1} X(m), e^{j\frac{2\pi m n}{N}}$$

与 DFT 相比有两处区别：

• 多系数 $1/N$
• 指数项为 $+j\frac{2\pi m n}{N}$

利用欧拉公式展开：

$$x(n) = \frac{1}{N}\sum_{m=0}^{N-1} X(m)\left[ \cos\left(\frac{2\pi m n}{N}\right)+ j\sin\left(\frac{2\pi m n}{N}\right) \right]$$

IDFT 以 $N$ 个复频域样本为输入，恢复出 $N$ 个时域离散样本。