快速傅里叶变换为什么这么快？

好的，我们来从 傅里叶变换 开始，一步步深入到 快速傅里叶变换（FFT） 。

傅里叶变换 (Fourier Transform, FT)

想象一下，你听到一段复杂的音乐，里面有钢琴声、小提琴声、鼓声等等。虽然这些声音混合在一起，但你的大脑能够分辨出不同乐器发出的不同音高。耳朵能够分辨不同频率，最重要的部分在于内耳耳蜗中的基底膜。这条膜具有从基部到顶部逐渐变化的物理特性（硬度和宽度），使得 不同频率的声音能在基底膜的不同位置引起最大振动：高频声音主要振动靠近基部的窄硬区域，低频声音则主要振动靠近顶部的宽软区域。这种频率在基底膜上的空间定位是核心，它使得特定位置的毛细胞被特定频率的声音所激活，并将这种基于位置的信息转化为神经信号传递给大脑，大脑根据信号来源的位置最终解读出我们所感知的不同音高。人耳对声音的感知涉及到心理声学的内容，如果你对此感兴趣，可以阅读我们的其他文章（放上文章链接，如果有的话）。

傅里叶变换在数学和信号处理中的作用与此类似：它是一种将信号（可以是时间相关的声音、图像中的亮度变化，或其他任何类型的函数）分解成其组成频率（或称为频谱分量）的技术。任何复杂的周期性信号都可以表示为一系列不同频率、不同振幅和不同相位的正弦波（或余弦波）的叠加。傅里叶变换将这个思想推广到了非周期信号，它可以告诉我们一个信号中包含哪些频率成分，以及这些频率成分的强度（振幅）和相位。正向傅里叶变换: 将信号从时域（或空域）转换到频域。

使用MathJax渲染数学公式 \[ X(f) = \int_{-\infty}^{\infty} x(t) e^{-j2\pi ft} dt \]

其中：x(t) 是时域信号；X(f) 是频域表示，它是一个复数函数，包含了频率 f 的振幅和相位信息；f 是给定的频率；t 是时间；e-j2πft 是复平面上的单位周期函数，通常称之为基函数。事实上 e-j2πft 就类似于人耳基底膜中的一个特定位置，它只对特定的频率敏感（实际上它只对周期和它一样为1/(2πf)的成分敏感，和其他频率成分是正交的）。可以通过以下公式做一个不严谨的验证，假设 x(t)=ej2πgt ，其中g≠f，则：

使用MathJax渲染数学公式 \[ \begin{aligned} X(f) &= \int_{-\infty}^{\infty} x(t) e^{-j2\pi ft} dt \\ &= \int_{-\infty}^{\infty} e^{j2\pi gt} e^{-j2\pi ft} dt \\ &= \int_{-\infty}^{\infty} e^{j2\pi (g-f)t} dt \\ &= \delta(g-f)\\ &=0\\ \end{aligned} \]

只有当上式中 g=f 时，积分才不为 0。知道了一个函数的频率成分和对应的相位，也就是 X(f)，以此推导出时域上的函数是很自然的，只需要将这些函数的周期和相位作为初始条件在时域上反演即可。想象你已经知道一个三角函数的周期和初始时刻的相位，其在整个时域的值也就确定了，而 X(f) 不过是众多三角函数和对应相位的集合罢了。可以用以下函数求对应 X(f) 的时域表示：

使用MathJax渲染数学公式 \[X(t) = \int_{-\infty}^{\infty} X(f) e^{j2\pi ft} df \]

当将连续信号在时间上进行采样，得到一个离散的时间序列时，就不能直接使用傅里叶变换。这时，需要离散傅里叶变换（DFT），DFT 是傅里叶变换在离散时域和离散频域上的对应形式。将傅里叶变换推广到离散信号也是很自然的（至少形式上是这样）：

使用MathJax渲染数学公式 \[X(f) = \int_{-\infty}^{\infty} x(t) e^{-j2\pi ft} dt \quad \Rightarrow \quad X[k] = \sum_{n=0}^{N-1} x[n] e^{-j2\pi kn/N} \quad \text{for } k=0,1,\ldots,N-1 \]

相比于连续函数傅里叶变换，x(t) 变为了 DFT 中 x[n] 的离散采样信号，其有 N 个采样点，积分变为了求和，f 中连续的周期变为了离散周期k。DFT 假设输入的离散序列 x[n] 是一个无限周期序列的一个主值周期（这可以通过对信号在 [-∞,∞] 周期延拓来做验证，延拓后的信号 DFT 和初始信号的 DFT 结果相同）。由于 x[n] 确实不一定是一个无限周期序列中的一个周期，因此 DFT 会引入一些误差，通常称为频谱泄露，可以通过窗函数进行抑制，本文不过多表述。 和连续傅里叶变换一样，e-j2πkn/N 也被称为核函数，或者“旋转因子”。设 WN=e-j2π/N，则核函数为 WNkn，离散傅里叶变换为：

使用MathJax渲染数学公式 \[X[k] = \sum_{n=0}^{N-1} x[n] W_N^{kn} \]

WN 的一个很有意思的性质是 WNN/2=0：

使用MathJax渲染数学公式 \[ \begin{aligned} W_N^{N/2} &= \left(e^{-j2\pi/N}\right)^{N/2} \\ &= e^{-j2\pi (N/2)/N} \\ &= e^{-j\pi} \\ &= \cos(-\pi) + j\sin(-\pi) \\ &= -1 \end{aligned} \]

上式使用了欧拉公式 ejθ=cos(θ)+jsin(θ) 。DFT 使得我们能够在计算机上对采样后的信号进行频域分析。它是数字信号处理的基石之一。X[k] 中的每个元素都是一个复数，表示在频率 fk=k⋅(f_s/N) 处的信号分量的振幅和相位（其中 fs 是采样频率）。

计算量问题

仔细观察 DFT 的公式，对于每一个 X[k] 的计算，需要 N 次复数乘法（x[n] 乘以 e-j2πkn/N)， N-1 次复数加法。由于我们需要计算 N 个不同的 X[k]（从 k=0 到 N-1），所以总的计算量大约是 N⋅N+N⋅(N-1)≈2N2 次。因此，DFT 的计算复杂度是 O(N2) 。当 N 很大时，举个例子，对于一个4秒钟的信号，采样率为44100/s，N=44100×4=176400，大约需要 2×176002≈622 亿次计算，这显然不能接受。

快速傅里叶变换 (Fast Fourier Transform, FFT)

FFT 不是一种新的变换，而是计算 DFT 的一种或一类高效算法。它能够显著减少计算 DFT 所需的乘法和加法次数，从而大大提高计算速度。

FFT 算法的核心思想是分治法 (有些人也称为分而治之)。它巧妙地利用了 DFT 计算中旋转因子 WNkn 的对称性和周期性，将一个大规模的 DFT 分解成若干个较小规模的 DFT，然后将这些小规模 DFT 的结果组合起来得到最终的 DFT 结果。

以最著名的 Cooley-Tukey 算法为例解释其原理。其基本步骤如下，假设 N 是2的整数次幂，即 N=2M：

  • 分解 : 将长度为 N 的输入序列 x[n] 分成两个长度为 N/2 的子序列：
    · 偶数索引的序列： g[m]=x[2m]
    奇数索引的序列： h[m]=x[2m+1]
  • 变换 : 分别计算这两个子序列的 N/2 点DFT：
    G[k] ( g[m] 的 N/2 点DFT)
    H[k] ( h[m] 的 N/2 点DFT)
  • 合并 : 利用旋转因子W_N^kn将 G[k] 和 H[k] 合并得到原始序列的 N 点DFT X[k]。注意DFT的定义是 X[k]=∑_n=0N-1x [n] WNkn。可以将其拆分为偶数项和奇数项：
使用MathJax渲染数学公式 \[ 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)} \]

也就是：

使用MathJax渲染数学公式 \[ X[k] = \sum_{m=0}^{N/2-1} g[m] \left(W_N^2\right)^{km} + W_N^k \sum_{m=0}^{N/2-1} h[m] \left(W_N^2\right)^{km} \]

由于 WN2=e-j2π/(N/2) =WN/2 ，上式变为：

使用MathJax渲染数学公式 \[ X[k] = G[k’] + W_N^k H[k’] \]

这里 G[k’] 和 H[k’] 是 N/2 点的DFT。当 k 从 0 到 N-1 时， k’ 在 N/2 点DFT中会周期性重复。

具体来说，对于 k=0,1,…,N/2-1：

使用MathJax渲染数学公式 \[ X[k] = G[k’] + W_N^k H[k’] \]

对于 k=N/2,…,N-1，可以令 k’=k-N/2。利用 WNk+N/2=WNk WNN/2=-WNk （使用式 2）和 G[k+N/2]=G[k] (因为 G 是 N/2 点 DFT，周期为 N/2) 以及 H[k+N/2]=H[k]:

使用MathJax渲染数学公式 \[ X\left[k+\frac{N}{2}\right] = G[k] – W_N^k H[k] \]

公式(3)(4)构成了所谓的“蝶形运算”，是 FFT 算法的基本计算单元，从这里就可以看出计算 X[k] 只需要计算两个一半信号的 G[k] 和 WNk H[k] ，并将他们组合即可。

  • 递归：当在整个周期 X[k] 时，实际上只需要计算一半点的 G[k] 和一半点的 WNk H[k] 。上述分解和合并的过程可以递归地应用于子序列的DFT计算，每次问题规模减半，直到DFT的规模减小到非常小（例如只剩1点或2点），这时问题将很容易解决。如果 N=2M，由于每次问题规模减半，那么这个递归过程会进行 M=log2 N 次。
计算量的降低：
  1. 一个 N 点的 DFT 被分解为两个 N/2 点的 DFT（计算G[k]和H[k]），加上 N/2 次复数乘法 (计算 WNk H[k] 和 N 次复数加/减法，整个操作的复杂度是O(N)。
  2. 设 C(N) 是计算 N 点 FFT 的运算量。则有递归关系： C(N)=2C(N/2)+O(N) (其中 O(N) 代表上面蝶形运算的复杂度)。
  3. 解这个递归关系，可以假设一个解，或者使用算法分析中的主定理，得到的总计算复杂度大约是 O(NlogN)，而 DFT 是 O(N2 )。

对比一下 O(N2 ) 和 O(NlogN)：当 N 更大时，例如 N=220 (约一百万点，常见于高分辨率图像或较长音频)：

  • DFT: (220 )2=240 次运算。
  • FFT: 220×20 次运算，虽然依旧很大，但比 240 小了非常非常多。速度提升了约 220/20≈52428 倍。

除了使用分而治之外，还有一些其他的方法加速运算，如当递归时问题分割到一定小的程度（此时函数调用的开销已经大于直接运算小规模问题的 DFT），递归停止并直接进行朴素的 DFT 计算，这种思想在很多算法中都可以得到体现，比如归并排序，小规模问题上的直接运算也更适合运用 SIMD 和多线程等技术。此外FFT也不止本文所述的 Cooley-Tukey 算法，还有一些基于其的变体如分裂基法、素因子法和用于处理素数长度的雷德法等。常见的快速傅里叶变换库如 FFTW 都集成了上述算法及优化，并且可以根据特定的输入规模、数据类型以及硬件环境进行测试并动态地选择不同的算法组合。因此不建议在通用计算机中使用自己实现的 FFT，而应该使用现有的FFT库，除非你需要在高度特化的计算环境中进行指向性的优化。