Featured image of post FFT学习笔记
学习

FFT学习笔记

这是一个副标题

引入

假设我们有两个多项式：

\( A(x) = a_0 + a_1x + a_2x^2 + ... + a_{n-1}x^{n-1} \)

\( B(x) = b_0 + b_1x + b_2x^2 + ... + b_{m-1}x^{m-1} \)

我们希望计算它们的乘积：

\( C(x) = A(x) \times B(x) = c_0 + c_1x + c_2x^2 + ... + c_{(n+m-2)}x^{(n+m-2)} \)

其中，每个系数 \( c_k = \sum_{i=0}^{k} a_i b_{k-i} \)。

如果直接计算（卷积），我们需要计算大约 \( n \times m \) 次乘法和加法。如果 \( n \) 和 \( m \) 都很大，这个 \( O(n^2) \) 的复杂度会非常慢。

而快速傅里叶变换（FFT）提供了一种将复杂度降至 \( O(n \log n) \) 的方法。

变换

多项式的标准表示法是系数表示法（Coefficient Representation），即我们上面用的 \( (a_0, a_1, ..., a_{n-1}) \)。还有一种等价的表示法，叫做点值表示法（Point-Value Representation）。

众所周知，两点唯一确定一条直线。那么三个点，四个点，甚至更多点呢？已有证明，一个 \( n-1 \) 次多项式可以由其在 \( n \) 个不同点 \( (x_0, x_1, ..., x_{n-1}) \) 上的取值 \( (A(x_0), A(x_1), ..., A(x_{n-1})) \) 唯一确定，咱知道就好。

点值表示法下的多项式乘法变得极其简单： 假设我们有多项式 \( A(x) \) 和 \( B(x) \)。

  1. 为 \( A(x) \) 取 \( N \) 个点：\( (x_k, A(x_k)) \)，\( k=0,1,...,N-1 \)
  2. 为 \( B(x) \) 取同样的 \( N \) 个点：\( (x_k, B(x_k)) \)，\( k=0,1,...,N-1 \) 注意，这里 \( N \) 必须至少是 \(deg(A) + deg(B) - 1 \)，否则乘积多项式 \( C(x) \) 无法被唯一确定。
  3. 那么乘积 \( C(x) = A(x) \times B(x) \) 在这 \( N \) 个点上的值就是： \( (x_k, C(x_k)) = (x_k, A(x_k) \times B(x_k)) \)，\( k=0,1,...,N-1 \)

此时点值相乘只需要 \( O(N) \) 的时间！这可比 \( O(N^2) \) 好多了。

于是，整个多项式乘法的策略就变成了：系数表示 → 点值表示 → 点值相乘 → 系数表示

  1. 求值（Evaluation）：将 \( A(x) \) 和 \( B(x) \) 从系数形式转换为点值形式。
  2. 点乘（Pointwise Multiplication）：计算 \( C(x_k) = A(x_k) \times B(x_k) \)。
  3. 插值（Interpolation）：将 \( C(x) \) 从点值形式转换回系数形式。

现在的核心问题是：如何快速地在系数表示和点值表示之间进行转换？ 如果随便选 \( N \) 个点，第一步（求值）和第三步（插值）的复杂度仍然是 \( O(N^2) \)（因为计算每个 \( A(x_k) \) 需要 \( O(N) \) 时间）。

FFT 的精妙就在于它精心选择了一组特殊的点，并利用分治策略，将求值和插值的过程都优化到了 \( O(N \log N) \)。

DFT 与单位根

FFT 选择的这组特殊的点，是单位根（Roots of Unity）。

离散傅里叶变换（DFT） 就是将一个多项式 \( A(x) \) 在 \( N \) 个 \( N \)次单位根 \( \omega_N^0, \omega_N^1, ..., \omega_N^{N-1} \) 上进行求值，得到的结果 \( (A(\omega_N^0), A(\omega_N^1), ..., A(\omega_N^{N-1})) \) 就是其DFT系数。

\( N \)次单位根是方程 \( z^N = 1 \) 的 \( N \) 个复数解。它们可以被表示为： \( \omega_N^k = e^{2\pi i k / N} = \cos(\frac{2\pi k}{N}) + i \sin(\frac{2\pi k}{N}) \)， 其中 \( k = 0, 1, ..., N-1 \)，\( i \) 是虚数单位。

单位根有几个极其重要的性质，是FFT分治的基础：

  1. 消去引理: \( \omega_{2N}^{2k} = \omega_N^k \)
  2. 折半引理: 如果 \( N \) 是偶数，那么 \( \omega_N^{k + N/2} = -\omega_N^k \)
  3. 求和引理：稍后在逆变换中会用到。

分治的魔法

FFT 是快速计算 DFT 的算法。我们假设 \( N \) 是 2 的幂（如果不是，可以补零到最近的 2 的幂）。

我们将多项式 \( A(x) \) 按奇偶项分成两个新的、规模减半的多项式：

\( A^{[0]}(x) = a_0 + a_2x + a_4x^2 + ... + a_{N-2}x^{N/2 - 1} \) （偶数索引项）

\( A^{[1]}(x) = a_1 + a_3x + a_5x^2 + ... + a_{N-1}x^{N/2 - 1} \) （奇数索引项）

很容易发现，原多项式可以表示为： \( A(x) = A^{[0]}(x^2) + x A^{[1]}(x^2) \)

现在，我们的问题从“求 \( A(x) \) 在 \( N \) 个点 \( (\omega_N^0, \omega_N^1, ..., \omega_N^{N-1}) \) 上的值”，变成了“求 \( A^{[0]}(x) \) 和 \( A^{[1]}(x) \) 在 \( N/2 \) 个点 \( ((\omega_N^0)^2, (\omega_N^1)^2, ..., (\omega_N^{N-1})^2) \) 上的值”。

根据消去引理：\( (\omega_N^k)^2 = \omega_{N/2}^k \) 并且，\( (\omega_N^{k + N/2})^2 = (\omega_N^k \omega_N^{N/2})^2 = (\omega_N^k \cdot -1)^2 = (\omega_N^k)^2 = \omega_{N/2}^k \)

这意味着，我们需要求值的 \( N \) 个平方点，实际上只有 \( N/2 \) 个不同的值：\( \omega_{N/2}^0, \omega_{N/2}^1, ..., \omega_{N/2}^{N/2 - 1} \)！

于是，我们可以进行分治递归：

设 \( y_k^{[0]} = A^{[0]}(\omega_{N/2}^k) \) 设 \( y_k^{[1]} = A^{[1]}(\omega_{N/2}^k) \) 这些是规模为 \( N/2 \) 的DFT问题。

那么，原问题 \( A(\omega_N^k) \) 的解可以通过合并子问题的解得到：

  1. 对于 \( k = 0, 1, ..., N/2 - 1 \):

    \( A(\omega_N^k) = A^{[0]}(\omega_{N/2}^k) + \omega_N^k \cdot A^{[1]}(\omega_{N/2}^k) = y_k^{[0]} + \omega_N^k y_k^{[1]} \)

  2. 对于 \( k = N/2, ..., N-1 \)（令 \( k’ = k - N/2 \)）:

    \( \omega_N^{k} = \omega_N^{k’ + N/2} = -\omega_N^{k’} \) \( A(\omega_N^{k}) = A(\omega_N^{k’ + N/2}) = A^{[0]}(\omega_{N/2}^{k’}) + \omega_N^{k’ + N/2} \cdot A^{[1]}(\omega_{N/2}^{k’}) = y_{k’}^{[0]} - \omega_N^{k’} y_{k’}^{[1]} \)

这个过程就是FFT的核心！

复杂度分析：

  • 每层递归：我们将一个规模为 \( N \) 的问题分解为 2 个规模为 \( N/2 \) 的问题，并进行 \( O(N) \) 的合并操作（乘以 \( \omega_N^k \) 并相加）。
  • 递归树的高度是 \( \log_2 N \)。
  • 总复杂度：\( T(N) = 2T(N/2) + O(N) = O(N \log N) \)。

IDFT

现在我们有了 \( C(x) \) 的点值表示 \( ((\omega_N^0, C(\omega_N^0)), (\omega_N^1, C(\omega_N^1)), ..., (\omega_N^{N-1}, C(\omega_N^{N-1}))) \)，如何变回系数 \( (c_0, c_1, ..., c_{N-1}) \)？这个过程称为逆离散傅里叶变换（IDFT）。

神奇的是，IDFT的过程和DFT几乎一模一样。 DFT的矩阵形式是：

$$ \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_N & \omega_N^2 & \cdots & \omega_N^{N-1} \\ 1 & \omega_N^2 & \omega_N^4 & \cdots & \omega_N^{2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_N^{N-1} & \omega_N^{2(N-1)} & \cdots & \omega_N^{(N-1)(N-1)} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{N-1} \end{bmatrix} $$
这个范德蒙德矩阵 \( V \) 几乎是可逆的。其逆矩阵 \( V^{-1} \) 与 \( V \) 非常相似，只是把每个元素 \( \omega_N^k \) 替换为 \( \omega_N^{-k} \)，并且整体乘以系数 \( 1/N \)。

也就是说，计算IDFT的过程就是：

  1. 将点值向量 \( \vec{y} \) 作为输入。
  2. 将单位根 \( \omega_N^k \) 全部替换为其共轭 \( \omega_N^{-k} \)。
  3. 执行一遍和FFT完全相同的算法！
  4. 最后将得到的结果除以 \( N \)。

所以，IFFT 的复杂度也是 \( O(N \log N) \)。

小结

假设我们要计算 \( A(x) \times B(x) \)，它们的次数界分别为 \( n \) 和 \( m \)。

  1. 补零（Zero-padding）：

    • 确定次数界 \( N \geq n + m - 1 \)，并且 \( N \) 是 2 的幂。
    • 将 \( A(x) \) 和 \( B(x) \) 的系数向量长度都补到 \( N \)，后面补零。

  2. 求值（FFT）：

    • 对补零后的 \( A(x) \) 的系数向量执行FFT，得到其点值表示 \( \vec{A} = (A(\omega_N^0), A(\omega_N^1), ..., A(\omega_N^{N-1})) \)。
    • 对补零后的 \( B(x) \) 的系数向量执行FFT，得到其点值表示 \( \vec{B} = (B(\omega_N^0), B(\omega_N^1), ..., B(\omega_N^{N-1})) \)。

  3. 点乘（Pointwise Multiply）：

    • 计算 \( \vec{C} = \vec{A} \circ \vec{B} \)，即 \( C(\omega_N^k) = A(\omega_N^k) \times B(\omega_N^k) \) for \( k=0,1,...,N-1 \)。

  4. 插值（IFFT）：

    • 对点值向量 \( \vec{C} \) 执行逆FFT（即将单位根取共轭，做FFT，再除以N），得到的结果就是乘积多项式 \( C(x) \) 的系数向量 \( (c_0, c_1, ..., c_{N-1}) \)。

最终复杂度为 \( O(N \log N) \)，远优于直接的 \( O(N^2) \)。

这就是FFT为何在信号处理、图像处理、数值计算等众多领域成为基石算法的原因。它通过巧妙的数学变换和分治策略，极大地加速了卷积（多项式乘法）运算。

© 2025 ❤️M.H. · 静心写作，分享所知