try ai
科普
编辑
分享
反馈
  • 快速傅里叶变换 (FFT)

快速傅里叶变换 (FFT)

SciencePedia玻尔百科
核心要点
  • 快速傅里叶变换 (FFT) 将离散傅里叶变换的计算复杂度从难以处理的 O(N2)O(N^2)O(N2) 降低到了高效的 O(Nlog⁡N)O(N \log N)O(NlogN)。
  • 它采用“分治”策略,递归地将一个大变换分解为多个小变换,并通过蝶形运算合并结果。
  • 通过卷积定理,FFT 极大地加速了卷积运算。卷积是信号处理、地球物理学乃至大数乘法中的核心操作。
  • FFT 是谱方法背后的引擎,谱方法通过在频域中将微分运算转化为简单的乘法来求解复杂的微分方程。

引言

离散傅里叶变换 (DFT) 是一种基础的数学工具,它使我们能够将一个复杂的信号(如声波或金融数据)分解为其组成频率。虽然功能强大,但直接应用 DFT 有一个致命的局限:其 O(N2)O(N^2)O(N2) 的计算复杂度使其在处理现代科学和工程中常见的大型数据集时变得异常缓慢。这种“平方的暴政”制造了一道计算壁垒,阻碍了无数领域的分析工作。

本文介绍的快速傅里叶变换 (FFT) 是一种革命性的算法,它能计算出与 DFT 完全相同的结果,但效率却高得惊人。通过巧妙地利用数学对称性,FFT 克服了计算瓶颈,为以前不可能进行的分析打开了大门。在接下来的章节中,您将学习这个优雅算法的工作原理,并见证其带来的变革性影响。“原理与机制”一章将剖析为 FFT 带来惊人的 O(Nlog⁡N)O(N \log N)O(NlogN) 速度和卓越数值精度的“分治”策略。随后的“应用与跨学科联系”一章将展示 FFT 的广泛用途,揭示一个单一算法如何驱动从音频滤波、地震成像到求解物理学和金融学基本方程等各种应用。

原理与机制

想象一下,你正在聆听一个复杂的和弦。你的耳朵,这个神奇的生物工程杰作,会立刻将这团混乱的声波分解为其组成音符:或许是一个 C、一个 E 和一个 G。​​离散傅里叶变换 (DFT)​​ 就是能让我们对任何信号——无论是声波、无线电传输,还是股票价格的波动——做同样事情的数学工具。它接收一个随时间记录的信号(即“和弦”),然后告诉我们它是由哪些“音符”(纯频率)组成的,以及每个音符的音量有多大。

DFT 的定义非常直观,甚至可以说简单得有些迷惑性。对于一个有 NNN 个采样点的信号,我们称之为 x[n]x[n]x[n],我们可以通过对所有采样点求和来找出任何给定频率分量 X[k]X[k]X[k] 的强度,每个采样点都乘以一个旋转的复数,即“旋转因子” WN=exp⁡(−2πi/N)W_N = \exp(-2\pi i/N)WN​=exp(−2πi/N)。公式如下:

X[k]=∑n=0N−1x[n]WNnkX[k] = \sum_{n=0}^{N-1} x[n] W_N^{nk}X[k]=n=0∑N−1​x[n]WNnk​

为了得到完整的频谱,我们必须对 NNN 个可能的频率分量都进行一次这样的计算。这个方法很简单,但其中却隐藏着一个可怕的瓶颈。

平方的暴政

我们来思考一下这其中涉及的计算量。对于 NNN 个输出频率中的每一个,我们都必须遍历所有 NNN 个输入采样点。这意味着我们大约需要执行 N×N=N2N \times N = N^2N×N=N2 次复数乘法和加法。我们说其复杂度是 N2N^2N2 阶,记为 O(N2)O(N^2)O(N2)。

对于少量采样点来说,这不是问题。但在现实世界中,信号可能非常庞大。一秒钟的 CD 音质音频就有 44,100 个采样点。一个 O(N2)O(N^2)O(N2) 算法仅处理一秒钟的声音就需要大约 44,1002≈2044,100^2 \approx 2044,1002≈20 亿次运算!如果我们分析一张高分辨率图像或一段长长的科学数据,NNN 可能达到数百万,而 N2N^2N2 则会成为一个天文数字。一个用更好的算法只需一秒钟就能完成的计算,用直接的 DFT 方法可能需要数年。

这不仅仅是理论上的担忧。在一个 N=4096N = 4096N=4096(在现代计算中这是一个非常适中的尺寸)的网格上进行的简单模拟显示,直接 DFT 方法的速度可能比一种更巧妙的方法慢近 70 倍。这种 O(N2)O(N^2)O(N2) 的行为是一道计算壁垒,是“平方的暴政”,几十年来一直限制着我们分析复杂信号的能力。我们需要一种更快的方法。不是近似方法,而是一种数学上完全相同,但效率却高得多的方法来计算完全相同的变换。

伟大的分治策略

这个突破以​​快速傅里叶变换 (FFT)​​ 的形式出现,该算法是“分治”策略的一个绝佳范例。FFT 的天才之处(最著名的贡献者是 J.W. Cooley 和 John Tukey)在于意识到一个大规模的 DFT 可以被分解,并由小规模的 DFT 构建而成。

让我们看看这个魔法是如何实现的。我们从一个长度为 NNN 的序列的 DFT 求和开始(为简便起见,我们暂且假设 NNN 是 2 的幂):

X[k]=∑j=0N−1x[j]WNjkX[k] = \sum_{j=0}^{N-1} x[j] W_N^{jk}X[k]=j=0∑N−1​x[j]WNjk​

现在,我们做一个非常简单的操作:将这个和分成两部分,一部分是对偶数索引的采样点(j=2mj=2mj=2m)求和,另一部分是对奇数索引的采样点(j=2m+1j=2m+1j=2m+1)求和。

X[k]=∑m=0N/2−1x[2m]WNk(2m)+∑m=0N/2−1x[2m+1]WNk(2m+1)X[k] = \sum_{m=0}^{N/2-1} x[2m] W_N^{k(2m)} + \sum_{m=0}^{N/2-1} x[2m+1] W_N^{k(2m+1)}X[k]=m=0∑N/2−1​x[2m]WNk(2m)​+m=0∑N/2−1​x[2m+1]WNk(2m+1)​

这看起来有点乱,但一个奇妙的简化即将出现。看一看旋转因子。WN2=(exp⁡(−2πi/N))2=exp⁡(−2πi/(N/2))=WN/2W_N^2 = (\exp(-2\pi i/N))^2 = \exp(-2\pi i/(N/2)) = W_{N/2}WN2​=(exp(−2πi/N))2=exp(−2πi/(N/2))=WN/2​。复指数的这个美妙性质是关键所在!利用这一点,我们可以重写方程:

X[k]=∑m=0N/2−1x[2m]WN/2km⏟DFT of even part+WNk⋅∑m=0N/2−1x[2m+1]WN/2km⏟DFT of odd partX[k] = \underbrace{\sum_{m=0}^{N/2-1} x[2m] W_{N/2}^{km}}_{\text{DFT of even part}} + W_N^k \cdot \underbrace{\sum_{m=0}^{N/2-1} x[2m+1] W_{N/2}^{km}}_{\text{DFT of odd part}}X[k]=DFT of even partm=0∑N/2−1​x[2m]WN/2km​​​+WNk​⋅DFT of odd partm=0∑N/2−1​x[2m+1]WN/2km​​​

看看发生了什么!一个大小为 NNN 的 DFT 现在被表示为两个大小为 N/2N/2N/2 的 DFT:偶数索引采样点的 DFT(我们称其输出为 E[k]E[k]E[k])和奇数索引采样点的 DFT(称之为 O[k]O[k]O[k])。

公式变得异常简单:

X[k]=E[k]+WNkO[k]X[k] = E[k] + W_N^k O[k]X[k]=E[k]+WNk​O[k]

这告诉我们如何合并那些较小的结果。这个操作被称为​​蝶形运算​​,因其在图示中的形状而得名。为了让这个想法更具体,考虑一个简单但高度振荡的输入信号 x[n]=(−1)nx[n] = (-1)^nx[n]=(−1)n。当我们将其分解为偶数和奇数部分时,偶数采样点为 x[2n]=(−1)2n=1x[2n] = (-1)^{2n} = 1x[2n]=(−1)2n=1,奇数采样点为 x[2n+1]=(−1)2n+1=−1x[2n+1] = (-1)^{2n+1} = -1x[2n+1]=(−1)2n+1=−1。这个看似复杂的信号分解成了两个平凡的常数序列!。这就是抽取的威力:它能将一个复杂问题转化为多个简单得多的子问题。

我们还没完。这个公式给了我们输出频谱的前半部分(对于从 000 到 N/2−1N/2-1N/2−1 的 kkk)。那后半部分呢?在这里,另一个对称性来拯救我们。事实证明,对于后半部分,组合方式几乎完全相同:

X[k+N/2]=E[k]−WNkO[k]X[k+N/2] = E[k] - W_N^k O[k]X[k+N/2]=E[k]−WNk​O[k]

所以,通过两个大小为 N/2N/2N/2 的 DFT 和一组简单的蝶形加法与乘法,我们就能构建出大小为 NNN 的完整 DFT。关键在于:我们可以将同样的逻辑应用于那些更小的 N/2N/2N/2 大小的 DFT,将它们分解成 N/4N/4N/4 大小的 DFT,以此类推,直到最后只剩下大小为 1 的平凡 DFT。

对数奇迹

为什么这种递归分解如此之快?让我们来计算一下工作量。为了计算一个大小为 NNN 的 DFT,我们进行了两次大小为 N/2N/2N/2 的 DFT,然后在蝶形运算阶段执行了大约 NNN 次操作。如果我们用 T(N)T(N)T(N) 表示计算一个大小为 NNN 的 DFT 所需的时间,我们的递推关系大致是:

T(N)=2T(N/2)+cNT(N) = 2T(N/2) + cNT(N)=2T(N/2)+cN

其中 ccc 是某个常数。如果你展开这个递推关系,你会发现在递归的每一“层”,总工作量都大约是 cNcNcN。那么有多少层呢?由于我们在每一步都将问题规模减半,所以层数是 log⁡2N\log_2 Nlog2​N。因此,总工作量与 Nlog⁡2NN \log_2 NNlog2​N 成正比。

N2N^2N2 和 Nlog⁡2NN \log_2 NNlog2​N 之间的差异不仅仅是小小的改进,它是一种相变,是不可能与常规之间的区别。对于 N=1,048,576N=1,048,576N=1,048,576(一个 1024×10241024 \times 10241024×1024 的图像),N2N^2N2 超过一万亿。但 Nlog⁡2NN \log_2 NNlog2​N 仅约 2000 万。FFT 将一个需要超级计算机耗时数月的任务,变成你的笔记本电脑在几分之一秒内就能完成的事情。这个“对数奇迹”为现代数字信号处理打开了大门。

从优雅的递归到实用的速度

递归的“分治”思想很优美,但在实践中,进行数百万次函数调用可能会很慢。高性能 FFT 库使用的是等效的​​迭代算法​​。要理解它,我们需要问一个问题:如果我们不断将输入数组分解为偶数和奇数部分,那么原始的采样点最终会到哪里去?

答案很有趣。一个最初在索引 jjj 处的输入采样点,最终会到达由 jjj 的二进制表示的​​位倒序​​所给出的位置。例如,在一个大小为 N=8N=8N=8 的变换中,索引为 6(二进制 110)的采样点将与索引为 3(二进制 011)的采样点交换位置。

这导出了一个非常高效的迭代算法:

  1. 首先,根据这个位倒序置换对整个输入数组进行重新排序。这是将所有数据放到正确位置以便开始计算的“洗牌”过程。
  2. 然后,自底向上进行。在第一遍中,合并相邻的采样点对,以计算大小为 2 的 DFT。
  3. 在下一遍中,合并相邻的大小为 2 的 DFT 对,以形成大小为 4 的 DFT。
  4. 按此方式逐层进行,共 log⁡2N\log_2 Nlog2​N 个阶段,在每个阶段将正在构建的 DFT 的大小加倍,直到得到最终的、完整的大小为 NNN 的 DFT。

这种自底向上的迭代方法避免了递归的开销,并允许进行高度优化、内存高效的实现。它正是驱动我们手机、电脑和科学仪器中 FFT 的引擎。

对称之美

傅里叶变换的结构富含对称性,利用这些对称性可以带来更高的优雅度和效率。

考虑一下当我们的输入信号是实值时会发生什么,这几乎是所有来自物理世界的信号的情况。其产生的频谱具有一种称为​​共轭对称​​的特殊性质:频率分量 kkk 是分量 N−kN-kN−k 的复共轭。也就是说,X[k]=X[N−k]‾X[k] = \overline{X[N-k]}X[k]=X[N−k]​。这源于旋转因子自身的一个基本对称性:WNnk‾=WNN−nk\overline{W_N^{nk}} = W_N^{N-nk}WNnk​​=WNN−nk​。这种对称性意味着大约一半的输出是冗余的!如果你知道频谱的前半部分,你就能立即知道后半部分。针对实值输入的巧妙算法利用这一点,将计算量减少了近一半,这是数学提供的又一顿“免费午餐”。

当我们考虑​​逆变换​​——即从频谱重构原始时域信号的过程时,一个更惊人的统一性得以揭示。逆离散傅里叶变换 (IDFT) 的公式与正向 DFT 几乎完全相同,只是指数的符号改变了,并多了一个 1/N1/N1/N 的缩放因子。人们可能期望需要一个全新的算法来进行逆变换。但值得注意的是,我们并不需要。正向 FFT 算法只需稍作修改即可用于计算逆 FFT:你只需取输入的复共轭,运行正向 FFT,然后取结果的复共轭(并乘以 1/N1/N1/N)即可。完全相同的计算机制可以通过共轭操作“反向”运行。这是关于时域和频域深刻对偶性质的有力陈述。

意外的馈赠:精度的礼物

到目前为止,我们一直在称赞 FFT 惊人的速度。但它还有另一个同样重要的优点:它在数值上也更精确。

计算机以有限精度进行算术运算,这意味着每次计算都可能引入微小的​​舍入误差​​。在一个有数十亿次操作的算法中,这些微小的误差会累积起来,破坏最终结果。正是在这一点上,FFT 的效率带来了第二个红利。

直接 DFT 有 O(N2)O(N^2)O(N2) 次操作,为这些误差的累积提供了更多的机会。而 FFT 只有 O(Nlog⁡N)O(N \log N)O(NlogN) 次操作,机会要少得多。因此,FFT 不仅能更快地产生答案,而且还能产生一个更干净的答案。数值实验一致表明,来自 FFT 的误差明显小于来自直接 DFT 的误差,尤其是在 NNN 很大时。

这不仅仅是一个经验观察。它已被严格证明。FFT 计算中误差的范数增长仅与 log⁡N\log NlogN 成正比,而对于直接 DFT,它可能增长得快得多。FFT 相对误差的最终理论界限非常简洁:它与 u⋅log⁡2Nu \cdot \log_2 Nu⋅log2​N 成正比,其中 uuu 是机器的单位舍入精度。这种对数增长异常缓慢,使得 FFT 成为一个非常稳定的算法。

快速傅里叶变换不仅仅是一个巧妙的算法。它证明了数学结构的力量。通过识别和利用隐藏在傅里叶变换定义中的对称性,它将一个计算上难以处理的问题转变为现代科学技术的基石。它不仅给了我们速度,还给了我们精度,这是来自纯数学领域一份稀有而美丽的礼物。

应用与跨学科联系

既然我们已经剖析了快速傅里叶变换,并看到了巧妙的分治策略如何使其速度快得惊人,我们就可以问一个最重要的问题:它有什么用?你会欣喜地发现,答案是:几乎无所不能。FFT 不仅仅是一个快速算法;它是一种审视世界的新视角,通过改变我们的观点,它解决了科学、工程乃至金融领域一系列惊人的问题。它是一个绝佳的例子,展示了一段优美的数学如何能够向外扩散,改变那些表面上彼此毫无关联的领域。

我们的应用之旅将是一次探索这种相互联系的旅程。我们将看到,用于过滤音频信号的相同思想可以用来乘以巨大的数字,而用于模拟遥远恒星湍流的工具,也同样是设计新材料和在华尔街为股票定价的核心。

问题的核心:快速卷积

FFT 最直接、最基础的应用源于​​卷积定理​​。正如我们所见,这个定理是一段数学魔法:在时域或空域中复杂、混乱的卷积过程,在频域中变成了简单的、逐元素的乘法。许多物理相互作用和数学运算,其核心都是卷积。

想一个​​数字信号处理​​中的简单过程。一个音频信号,它只是一串表示压力随时间变化的数字,通过一个滤波器——也许是为了增强低音或消除嘶嘶声。输出信号是原始信号与滤波器“冲激响应”的卷积。直接进行此操作,需要为输出中的每一个点进行大量的乘法和加法运算。但有了 FFT,我们可以变换信号和滤波器的响应,在频域中将它们逐点相乘,然后再变换回来。结果是相同的,但完成的工作量只是原来的很小一部分。为了使其完美工作,我们必须小心地用零填充我们的信号到合适的长度,以确保 DFT 的循环卷积能给我们实际想要的线性卷积,但为如此巨大的速度提升付出这个小代价是值得的。

同样的原理让我们能够深入探视地球内部。在​​地球物理学​​中,地震道被用来绘制地下结构图。一个震源,比如一次小型爆炸或一辆巨大的震源车,向地下发送一个声波(通常称为“子波”)。当这个子波传播时,它会从不同的岩石层反射回来。地表麦克风记录到的信号是一系列的回声。这个记录道可以很好地建模为源子波与地球“反射系数序列”的卷积——这是一个表示不同深度反射强度的数字列表。通过使用 FFT 执行这种卷积,地震学家可以创建合成地震道来检验他们的假设,甚至可以反向工作——一个称为反褶积的过程——试图从记录的数据中揭示隐藏的反射系数序列,从而描绘出我们脚下世界的图景。

从信号到符号:计算领域的一场革命

到目前为止,我们谈论的“信号”都是存在于时间或空间中的事物。但如果这一串数字代表的是更抽象的东西呢?如果它们是多项式的系数呢?例如,序列 [a_0, a_1, a_2] 代表多项式 a0+a1x+a2x2a_0 + a_1x + a_2x^2a0​+a1​x+a2​x2。

这里有一个惊人的发现:两个多项式相乘的行为与对其系数列表进行卷积完全相同。这是同一个数学运算!突然之间,我们的信号处理工具 FFT,变成了一种闪电般快速的代数运算方式。要乘以两个次数约为 nnn 的多项式,传统方法大约需要 n2n^2n2 次运算。而使用 FFT,我们变换系数列表,将它们相乘,然后变换回来,所有这些都可以在大约 O(nlog⁡n)O(n \log n)O(nlogn) 的时间内完成。

我们可以更进一步。什么是大整数?像 528 这样的数字只是多项式 5x2+2x1+8x05x^2 + 2x^1 + 8x^05x2+2x1+8x0 在 x=10x=10x=10 时求值的简写。这意味着我们可以通过将两个巨大整数的数字看作是两个多项式的系数来将它们相乘。我们使用 FFT 对数字序列进行卷积,得到乘积多项式的系数。最后一步是“处理进位”,这只是在这个新的多项式上以 x=10x=10x=10 求值。这种方法,被称为 Schönhage–Strassen 算法(及其后续算法),是计算机科学领域的一场革命。一个源于分析波动的思想,从根本上改变了我们执行最基本算术运算的方式。

频率的微积分:求解宇宙的方程

FFT 还有另一个魔术,能将它带入微积分和自然基本定律的领域。许多这些定律都以偏微分方程 (PDE) 的形式表达,描述了量在空间和时间中的变化。求解这些方程是出了名的困难。

魔术在于:在傅里叶域中,微分运算变成了简单的乘法。对一个函数求导 ddx\frac{d}{dx}dxd​ 等价于将其傅里叶变换乘以 ikikik,其中 iii 是虚数单位,kkk 是波数(频率)。这是因为傅里叶级数的基函数,即正弦和余弦函数(或复指数),是微分算子的特征函数。

这一性质是求解 PDE 的​​谱方法​​的基础。为了计算一个复杂方程中的空间导数,我们可以对我们系统的状态(比如,网格上所有点的流体速度)进行 FFT,在频域中乘以适当的 ikikik 因子,然后使用逆 FFT 返回到物理空间,我们的导数就已经被精确地计算出来了。

这种方法在​​流体动力学​​等领域带来了变革。模拟湍流——流动气体或液体中涡旋的混沌之舞——是计算物理学的重大挑战之一。借助由 FFT 驱动的谱方法,研究人员可以在拥有数十亿个点的网格上进行直接数值模拟 (DNS)。速度的提升不仅仅是微小的改进,而是令人难以置信的。对于一个在 512×512×512512 \times 512 \times 512512×512×512 网格上的三维模拟,使用 FFT 比直接计算快了近​​五百万倍​​。没有 FFT,这些模拟将完全不可能实现。

同样的想法也驱动着原子尺度的模拟。在​​分子动力学​​中,计算成千上万或数百万带电粒子之间的长程静电力是一个计算上令人望而却步的 O(M2)O(M^2)O(M2) 问题。质点-网格方法,如质点-网格-Ewald (PME) 或质点谱 Ewald (PSE) 技术,通过将粒子电荷分布到一个均匀的网格上,使用 FFT 求解泊松方程得到电势,然后将力插值回粒子位置来克服这个问题。这将成本降低到可控的 O(Mlog⁡M)O(M \log M)O(MlogM),从而使那些对于药物发现、材料科学以及理解生命基本机制至关重要的模拟成为可能。

超越物理学:一种通用语言

FFT 的影响并不止于物理科学。它能快速计算大量结果的能力,使其在速度至关重要的领域中具有无价的价值。在​​计算金融​​中,一个金融期权(如在未来以特定价格购买股票的权利)的价值可以用 Black-Scholes 偏微分方程来描述。虽然像 Crank-Nicolson 格式这样的传统方法可以为单个行权价求解这个方程,但基于傅里叶的方法具有显著的优势。通过处理资产价格分布的特征函数,单个 FFT 就能同时为一系列广泛的行权价的期权定价。在高速的量化交易世界中,这种效率改变了游戏规则。

回到实验室,FFT 不仅作为计算引擎,还是一种复杂的分析工具。在​​材料科学​​中,透射电子显微镜 (TEM) 可以产生分辨率极高的图像,以至于晶体中单个原子列都清晰可见。如果我们拍摄这样一张图像并计算其二维 FFT,我们会得到一个反映原子周期性排列的亮点图案——这在数学上等同于一个衍射图样。

有趣的是,将这个计算出的 FFT 与来自同一晶体的实际物理选区电子衍射 (SAED) 图样进行比较。亮点的位置会匹配,证实了晶体结构。但亮点的强度以及精细的细节会有所不同。这些差异不是错误,它们是数据!它们告诉我们简单的数学 FFT 所忽略的物理学:显微镜透镜扭曲图像的方式(衬度传递函数)以及电子在较厚晶体内部复杂的多次散射(动力学)相互作用,这些相互作用产生了像菊池线这样美丽的图案。FFT 提供了一个理想的基准,我们可以借此理解更丰富的现实物理。

深入探究:代数之美

我们已经看到 FFT 做了这么多不同的事情。是否存在一个统一的原则,一个其非凡力量背后的更深层原因?答案是肯定的,它就在代数结构之中。

考虑一种特殊类型的矩阵,称为​​循环矩阵​​,其中每一行都是其上一行的循环移位。如果你写下表示卷积运算的矩阵,你会发现它就是一个循环矩阵。现在,秘密在于:所有循环矩阵都可以被傅里叶变换矩阵对角化。

这是一个深刻的陈述。它意味着如果你将基底变换为傅里叶基(这正是 FFT 所做的),任何以及所有的卷积运算都会变成简单的对角矩阵。乘以一个对角矩阵就是逐元素乘法。这就是卷积定理背后的终极“原因”。FFT 不仅仅是一个巧妙的算法;它是解开与循环对称性相关的一整类问题的自然、对角基底的钥匙。

从这个角度看,FFT 不再仅仅是处理信号或图像的工具。它是理解和操纵具有潜在周期性或循环性结构的系统的基本钥匙,无论这些系统是岩石层、多项式的系数,还是晶体中的原子。它是数学统一性及其描述我们世界的惊人、深远力量的美丽证明。