1965年：FFT让大规模傅里叶分析走向实用

詹姆斯·库利与约翰·图基于1965年4月发表的五页论文，推动了一项繁重计算成为常用方法：找出序列中的频率成分。他们的快速傅里叶变换把计算重新组织为更小、可以重复利用的部分。当序列长度为2的幂时，运算量大致按点数与其对数的乘积增长，而不是按点数的平方增长。更早的数学家已发现过相关方法。真正留下持久影响的是这一方法在电子计算中的有效实现和传播。速度提升让更大规模的计算变得可行，却不会改善原始测量，也不会消除数值计算与采样的限制。

封面：由AI生成的概念插图，表现合成信号及其频率成分。它既不描绘当年的计算机，也不展示实验数据。

听出信号内部的成分

想象一下，你在录制一个音乐和弦。麦克风给出的是一串振幅值，而你想知道的是其中有哪些音。傅里叶分析在这两种描述之间建立数学联系。研究振动或图像中的图案时，也会遇到类似的问题。

离散傅里叶变换（DFT）把有限序列表示为离散的频率成分。快速傅里叶变换（FFT）则是高效计算同一个变换的方法，并不会产生一种新的频谱。NumPy的技术文档把输入描述为时域信号，把输出描述为它在频域中的表示。DFT文档

实际困难在于重复计算。直接计算会把每一个输入点与每一个输出频率结合起来。点数翻倍，主要的运算量就增加到四倍。因此，当科学问题需要处理数千个点时，一条数学上简单的公式也可能成为计算瓶颈。

库利与图基提出了一般性的因式分解，以及适用于长度为2的幂的便捷实现，其中还包括重复使用输入数组存储空间的方法。论文报告了IBM 7094上的程序及运行时间，而不只是一个抽象提案。收稿日期为1964年8月17日，论文刊登在Mathematics of Computation的1965年4月号。这里尚未确认确切的出版日。原论文，大学提供的副本

一项有着漫长前史的突破

这个想法并非在1965年凭空出现。历史研究者Michael Heideman、Don Johnson和Sidney Burrus考察了更早的文献，在高斯的工作中找到了等价的分解。他们推断其写作时间约为1805年；原稿本身没有明确注明日期，直到高斯去世后的1866年才发表。他们的研究还讨论了后来的方法，包括Danielson与Lanczos在1942年的工作，以及Good的素因子方法。Good方法的约束与一般性的库利–图基因式分解不同。历史考察及可访问的作者副本

在1993年的回忆中，库利与图基描述了Richard Garwin在连接相关工作、推动其发展方面的作用。他们也承认更早的发现，并解释为什么大规模电子计算使这一方法具有价值。这是参与者的回顾性叙述，需要考虑记忆本身的局限。库利与图基的回忆

持久的转折在于，一种可重复利用的计算方法进入了实践。普林斯顿大学关于IEEE里程碑的介绍，把它与科学计算、信号处理、医学成像和数据传输联系起来。这些应用也依赖仪器、物理模型及其他算法；FFT提供了一个强有力的计算模块。后来的纪念与认可支持其历史重要性，但不能把近期纪念活动的日期当作发现日期。普林斯顿大学与IEEE里程碑

四个点：先看如何复用，再看一般公式

使用无量纲序列**[1, 2, 3, 4]。这个刻意缩小的例子是我们的教学计算，并非1965年论文中的数据。用i表示虚数单位，满足i² = −1**。不进行缩放、指数取负号的DFT给出四个输出：

频率索引	输出
0	10
1	−2 + 2i
2	−2
3	−2 − 2i

第一个输出把四个输入相加。其余输出则组合正、负及虚数权重。一个复数输出记录两个分量；它的模和相位共同描述一个频率成分。

现在，把偶数索引的输入[1, 3]与奇数索引的输入[2, 4]分开。每一组只需要求和与求差：

偶数索引组：E = [4, −2]。
奇数索引组：O = [6, −2]。

再把这些较短的结果组合起来。第一个蝶形运算得到4 + 6 = 10和4 − 6 = −2。第二个先用−i旋转奇数索引组的结果，得到**(−i)(−2) = 2i**，随后产生−2 + 2i与−2 − 2i两个输出。

“蝶形”是这种成对求和与求差运算的名称。示意图中的交叉线用于整理计算流程，不是波之间的物理连接。两个输出复用同样的中间结果。规模增大时，每个较短的变换还可以继续拆分：复用由一次性的捷径变为层次结构。

从具体例子到运算量的增长

对于N个输入值，我们采用如下约定：

$X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i nk/N}.$

这里，n是输入索引，k是输出索引，二者都从0取到N − 1。这些因子表示模为1的复数旋转。我们的正向变换不含归一化因子；逆变换使用正指数及1/N因子。还存在其他约定。原论文从正指数的傅里叶和出发，因此必须明确比较符号和归一化方式。现代约定

当N为偶数时，把求和按输入索引的奇偶拆开。将两个长度为N/2的变换写成Eₖ与Oₖ，并令W_N = exp(−2πi/N)。对于0 ≤ k < N/2，重构公式为：

$X_k = E_k + W_N^k O_k, \qquad X_{k+N/2} = E_k - W_N^k O_k.$

第二个等式利用了旋转半圈时因子变号这一点。于是，我们计算两个较短的变换，再进行与N成正比的组合操作。当N = 2ᵐ时重复拆分，会得到m = log₂N层。

每层的工作量与序列长度成正比。因此，总工作量按N log₂N增长。O(N log N)这一记号描述的是增长规律，而不是精确的处理器指令数或耗时秒数。这里用于教学的基2构造仅适用于2的幂；一般性的库利–图基因式分解也可用于其他合数长度。

真实计算机的速度还取决于数据布局、内存传输与实现方式。仅凭算术运算次数减少，不能确定你的设备上实际运行时间会缩短多少。

更快的计算无法解决什么

如果等间隔测量的间隔为Δt秒，那么采样率为1/Δt赫兹，频率格点间隔为1/(NΔt)赫兹。解释较大的输出索引时，需要考虑所采用的正、负频率排列顺序。频率约定

有限采样仍会限制能作出的推断。这里，t表示以秒为单位的时间，f表示以赫兹为单位的信号频率，fₛ = 1/Δt表示以赫兹为单位的采样率。例如，在t = n/fₛ处采样cos(2πft)，把f换成f + fₛ不会改变样本，因为增加的相位为2πn。任何更快的算法，都无法仅凭这些样本区分这两个信号。噪声和不充分的测量模型同样仍是测量问题。

数值运算还有另一种限制。复数旋转因子通常需要数值近似，改变运算顺序也可能改变舍入结果。FFTW关于精度的讨论强调了旋转因子精度的重要性，并记录了不同实现与环境之间的差异。数学上等价，并不保证浮点数的位模式相同。精度讨论

检查教学计算

我们的检查使用CPython 3.12.10，且只使用标准库。四点变换的根[1, −i, −1, i]能够精确表示，因此，对这些小规模整数输入采用完全相等的判据是合适的。不能把这一判据机械地用于更长的浮点变换。

在留有记录的检查中，执行了以下核心计算：

PYTHON

ROOTS = (1, -1j, -1, 1j)

def direct(x, roots=ROOTS):
    return [
        sum(x[n] * roots[(n*k) % 4] for n in range(4))
        for k in range(4)
    ]

def butterfly(x):
    even = (x[0] + x[2], x[0] - x[2])
    odd = (x[1] + x[3], x[1] - x[3])
    return [
        even[0] + odd[0], even[1] - 1j*odd[1],
        even[0] - odd[0], even[1] + 1j*odd[1],
    ]

在**[1, 2, 3, 4]、脉冲[1, 0, 0, 0]、常数序列[1, 1, 1, 1]**、全零序列，以及[1+i, −2, 3−i, 2i]上比较两个函数。五种情况的结果全部完全一致。第一种情况也与上表吻合。把根的虚部符号反过来，会改变非对称情况的频谱；检查因此拒绝了它，因为它采用的是不同约定。

这验证的是固定的教学案例，既不证明每一种实现都正确，也不构成对历史性能的独立复现。原论文没有提供完成这种复现所需的完整程序、基准测试输入和完整的计时条件。

本文由AI系统根据所引资料撰写并翻译，并按News编辑方法进行自动检查。教学计算已按文中说明执行；这不等于同行评审、人工科学验证或对历史程序的复现。