快速傅里叶变换 第一部分：库利-图基算法

# 快速傅里叶变换 第一部分：库利-图基算法 来源: https://connorboyle.io/2025/09/11/fft-cooley-tukey.html 2025年9月11日，作者：Connor Boyle 标签: mathematics (https://connorboyle.io/archive.html#mathematics) software (https://connorboyle.io/archive.html#software) ***我计划撰写一系列关于快速傅里叶变换（FFT）算法的文章。第一篇将介绍库利-图基（Cooley-Tukey）算法，这是原始且最著名的 FFT 算法。*** ## 离散傅里叶变换 如果 $x$ 是一个长度为 $|x|$ 且起始索引为 0 的复数序列，那么 $x$ 的离散傅里叶变换（DFT），记作 $\mathcal{F}\{x\}$，定义如下： $$|\mathcal{F}\{x\}| = |x|$$ $$\mathcal{F}\{x\}[k] = \sum_{j=0}^{|x|-1} x[j] \cdot e^{-i 2 \pi jk \frac{1}{|x|}}$$ 由于复数指数运算在傅里叶变换中非常常见，我们定义一个辅助项 $W_N$ 如下： $$W_N \triangleq e^{-i 2 \pi \frac{1}{N}}$$ 即 $W_N = W_N^1$ 是复平面上从 1 开始的 $\frac{1}{N}$ 圈旋转。$W_N^2$ 是复平面上的 $\frac{2}{N}$ 圈旋转，以此类推。代入原始的离散傅里叶变换定义，我们得到： $$\mathcal{F}\{x\}[k] = \sum_{j=0}^{|x|-1} x[j] \cdot W_{|x|}^{jk}$$ *朴素地*（Naïvely）计算这个方程，针对 DFT 的每一个不同的输出频率桶（$k = 0, 1, \ldots, |x| - 2, |x| - 1$），都需要对信号中的 $|x|$ 个样本进行复数乘积求和，因此任何朴素的 DFT 算法的时间复杂度均为 $O(|x|^2)$。 ## 库利-图基算法 如果 $|x|$ 是一个合数，我们可以选择两个自然数 $r$ 和 $d$，使得： $$|x| = r \cdot d$$ 这允许我们将对 $j$ 的单一求和转换为嵌套求和： $$\mathcal{F}\{x\}[k] = \sum_{j_1=0}^{d-1} \sum_{j_0=0}^{r-1} x[j_1 r + j_0] \cdot W_{|x|}^{(j_1 r + j_0)k}$$ $$= \sum_{j_0=0}^{r-1} \sum_{j_1=0}^{d-1} x[j_1 r + j_0] \cdot W_{|x|}^{(j_1 r + j_0)k}$$ 类似地，我们可以定义变量 $k_0$ 和 $k_1$，使得 $k = k_1 d + k_0$。令： $$k_1 \triangleq \lfloor \frac{k}{d} \rfloor$$ $$k_0 \triangleq k - k_1 d$$ 换句话说，$k_1$ 是 $k$ 除以 $d$ 的欧几里得除法 (https://en.wikipedia.org/wiki/Euclidean_division)[^1] 的商，$k_0$ 是余数。 这允许我们再次重新表述离散傅里叶变换： $$\mathcal{F}\{x\}[k] = \sum_{j_0=0}^{r-1} \sum_{j_1=0}^{d-1} x[j_1 r + j_0] \cdot W_{|x|}^{(j_1 r + j_0) (k_1 d + k_0)}$$ $$= \sum_{j_0=0}^{r-1} \sum_{j_1=0}^{d-1} x[j_1 r + j_0] \cdot W_{|x|}^{j_1 r k_1 d + j_1 r k_0 + j_0 (k_1 d + k_0)}$$ $$= \sum_{j_0=0}^{r-1} \sum_{j_1=0}^{d-1} x[j_1 r + j_0] \cdot W_{|x|}^{j_1 k_1 |x| + j_1 r k_0 + j_0 k}$$ $$= \sum_{j_0=0}^{r-1} \sum_{j_1=0}^{d-1} x[j_1 r + j_0] \cdot W_{|x|}^{j_1 k_1 |x|} \cdot W_{|x|}^{j_1 r k_0} \cdot W_{|x|}^{j_0 k}$$ 由于 $W_{|x|}^{j_1 k_1 |x|} = (e^{-i \frac{2 \pi |x|}{|x|}})^{j_1 k_1} = 1^{j_1 k_1} = 1$，因此： $$\mathcal{F}\{x\}[k] = \sum_{j_0=0}^{r-1} \sum_{j_1=0}^{d-1} x[j_1 r + j_0] \cdot W_{|x|}^{j_1 r k_0} \cdot W_{|x|}^{j_0 k}$$ 此时，我们可以将 $x$ 的元素分割成对应于模类的子序列。令 $x_r^{j_0}$ 为一个序列，其元素等于 $x$ 中索引模 $r$ 同余于 $j_0$ 的元素。更正式地说，这些序列（总共有 $r$ 个）可以定义如下： $$x_r^{j_0}[j_1] = x[j_1 r + j_0]$$ $$|x_r^{j_0}| = \frac{|x|}{r} = d$$ 代入这个序列定义，我们得到：[^2] $$\mathcal{F}\{x\}[k] = \sum_{j_0=0}^{r-1} \sum_{j_1=0}^{d-1} x_r^{j_0}[j_1] W_{|x|}^{k_0 j_1 r} W_{|x|}^{j_0 k}$$ $$= \sum_{j_0=0}^{r-1} \sum_{j_1=0}^{d-1} x_r^{j_0}[j_1] W_d^{k_0 j_1} W_{|x|}^{j_0 k}$$ $$= \sum_{j_0=0}^{r-1} \mathcal{F}\{ x_r^{j_0} \}[k_0] W_{|x|}^{j_0 k}$$ 让我们考虑在这种表述下计算离散傅里叶变换需要多长时间： - 首先，我们需要预计算 $r$ 个不同序列 $x_r^{j_0}$ （即 $\mathcal{F}\{ x_r^0 \}, \mathcal{F}\{ x_r^1 \}, \ldots \mathcal{F}\{ x_r^{r-1} \}$）的离散傅里叶变换，每个序列的长度为 $d$。假设每个 DFT 都是朴素执行的，这将需要 $O(r \cdot d^2) = O(|x| \cdot d)$ 次操作。 - 然后，我们需要计算 $\mathcal{F}\{x\}$ 的 $|x|$ 个值，每个值都需要对 $\mathcal{F}\{ x^{(j_0)} \}[k_0]$ 的 $r$ 个不同 $j_0$ 值进行求和，需要 $O(|x| \cdot r)$ 次操作。 将这两个子程序加在一起，需要 $O(|x| \cdot d + |x| \cdot r) = O(|x| \cdot (d + r))$ 次操作，根据 $r$ 和 $d$ 的值，这可能比原始朴素表述的 $O(|x|^2)$ 复杂度有显著改进。更重要的是，这种操作可以递归应用。具体而言，每个长度为 $d$ 的序列 $x_r^{j_0}$ 的 $r$ 个离散傅里叶变换都可以分解为 $r'$ 个长度为 $d'$ 的傅里叶变换，前提是存在两个自然数使得 $d = r' \cdot d'$。[^3] 在最理想[^4] 的情况下，其中 $|x| = 2^n$，$n \in \mathbb{N}$，计算库利-图基算法将需要 $O(|x| \cdot (2 + 2 + \ldots + 2)) = O(|x| \cdot 2 \cdot \log_2(|x|)) = O (|x| \cdot \log(|x|))$ 次操作。 库利-图基算法也可用于计算逆离散傅里叶变换，只需非常轻微的修改。事实上，原始的库利-图基论文（参见“相关阅读”） specifically 描述的是一种计算*逆*离散傅里叶变换的算法，而不是“前向” DFT。我将把库利-图基 iDFT 算法留作读者的练习。 但是，请注意库利-图基算法对于长度为素数的输入序列没有加速效果，并且当输入长度的因子包含大素数时，提供的加速效果也相对有限。为了有效地计算素数甚至非高合成数长度序列的 DFT，我们需要额外的算法。然而，最终，这些其他的 FFT 算法通常依赖于库利-图基算法进行部分计算。我计划在未来的博客文章中至少介绍其中一种技术——Bluestein 算法。 ### 库利-图基交互可视化 此可视化展示了如何使用库利-图基算法计算某个信号 $x$ 的离散傅里叶变换。最底部的黑色方块是输入信号。虽然 DFT 可以应用于复数信号，但为了简单起见，我将输入信号的样本值限制为实数（这模仿了一些现实世界的应用，例如对音频录音执行 DFT）。你可以点击并拖动输入方块来更改它们的值。 灰色圆圈以及其中有时可见的白色“时钟指针”代表某个 $N$（例如 $|x|$、$d$、$d'$ 等）和某个 $x$ 的复数指数 $W_N^x = e^{-2 i \pi \frac{x}{N}}$。这些复数指数（等效于复平面上的旋转）被应用于相关的输入值。与通常的惯例不同，我决定将复平面的实部显示为垂直方向（“向上”为正实部），虚部显示为水平方向（“向右”为正虚部）。你可以将鼠标悬停在某个“旋转”框上，查看输入值是从哪里绘制的。 这些旋转后的输入值相加，以计算离散傅里叶变换的一个元素。在可视化中将鼠标悬停在离散傅里叶变换的白色“输出”框上，以突出显示它求和来源的“旋转”框列。 $|x| =$ 可用因子： ## 关于措辞的一点牢骚 我注意到许多人在谈论离散傅里叶变换时有一种令人恼火且令人困惑的倾向——包括已发表的作者 (https://www.google.com/books/edition/The_Sparse_Fourier_Transform/I4ZTDwAAQBAJ?hl=en&gbpv=1&dq=%22the+FFT+of%22&pg=PT78&printsec=frontcover)——当他们想表达“离散傅里叶变换”（或“DFT”）时，经常使用短语“快速傅里叶变换”（或者更常见的是缩写“FFT”）。我认为这是错误且令人困惑的；要理解为什么，想象你有一个列表： 将以下列表称为前一个列表的“归并排序”（mergesort）有意义吗？ 我认为大多数人会觉得这非常奇怪。我们可以说第二个列表是排序第一个列表的结果，但我们不知道用于排序列表的具体算法是什么。它可能是使用归并排序、快速排序、堆排序、冒泡排序、猴子排序或任何其他排序算法 (https://en.wikipedia.org/wiki/Sorting_algorithm#Popular_sorting_algorithms) 进行排序的（实际上，我刚才只是在脑海中排好了这个列表）。然而，即使我*确实*使用了归并排序，$y$ 仍然不会*是* $x$ 的“归并排序”，它仍然只是“排序 $x$ 的结果”。 同样，FFT 算法的输出不应被称为任何事物的“FFT”。理论上，使用库利-图基算法计算序列的 DFT 与使用朴素实现的 DFT 算法计算该 DFT 给出完全相同的结果。在实践中，库利-图基可能会给出略微更准确的结果，因为总计算量较少，因此浮点舍入错误的机会也较少。 这不仅让我感到恼火；我认为这会在公众中引起混淆。我的一位朋友——一位聪明的数学专业学生和软件工程师——最近问我：“如果我用快速傅里叶变换代替‘普通’傅里叶变换，我会失去什么？我一直只使用普通傅里叶变换”。他可能从许多人随意使用“FFT”一词的方式中推断出，“快速傅里叶变换”是某种近似值或与离散傅里叶变换或一般傅里叶变换相关但不同的概念，而不是计算它们的*算法*。 - An Algorithm for the Machine Calculation of Complex Fourier Series (https://web.stanford.edu/class/cme324/classics/cooley-tukey.pdf) by James W. Cooley & John W. Tukey *感谢我的朋友 Andre Archer 帮助校对本文的早期版本。任何错误均由本人负责。* **脚注：** 标签: *mathematics*-*software* --- [^1]: https://connorboyle.io/2025/09/11/fft-cooley-tukey.html#fn:euclidean-division [^2]: https://connorboyle.io/2025/09/11/fft-cooley-tukey.html#fn:dft-xj0 [^3]: https://connorboyle.io/2025/09/11/fft-cooley-tukey.html#fn:radix-one [^4]: https://connorboyle.io/2025/09/11/fft-cooley-tukey.html#fn:ideal