「离散傅里叶变换」和「离散傅里叶反变换」

数学

离散傅里叶变换(Discrete Fourier Transform, DFT)和离散傅里叶反变换(Inverse Discrete Fourier Transform, IDFT) 是鼎鼎大名的快速傅里叶变换(Fast Fourier Transform, FFT)的前置知识.

其中 FFT 用于加速两个多项式 $A(x)$、$B(x)$ 的乘积 $C(x)$ 的计算,DFT 和 IDFT 是 FFT 的两个中间步骤.

需要明确,给定的是两个多项式:

其中 $A(x)$ 是一个 $n$ 项多项式,$B(x)$ 是一个 $m$ 项多项式.

我们要求的多项式函数 $C(x)$ 为一个 $n+m-1$ 项的多项式(考虑最低项 $x$ 的幂次为 $1$,最高项为 $n-1+m-1=n+m-2$,共 $1+(n+m-2)=n+m-1$ 项).

要了解 DFT 和 IDFT,还需要一些前置知识.


多项式的系数表示法

我们平时用的多项式系数 $a_0,a_1,\dots,a_{n-1}$ 就属于多项式的系数表示法.

我们称 $n$ 维系数向量 $(a_0,a_1,\dots,a_{n-1})$ 为 $A(x)$ 的系数表示,其中:

对于 $A(x)$ 的系数表示 $(a_0,a_1,\dots,a_{n-1})$ 和 $B(x)$ 的系数表示 $(b_0,b_1,\dots,b_{m-1})$,要求 $C(x)$ 的系数表示,显然可以通过以下代码实现:

1
2
3
for (int i = 0; i < n; i++)
for (int j = 0; j < m; j++)
c[i + j] += a[i] * b[j];

得到的 $C(x)$ 为 $n+m-1$ 项多项式.

这个朴素算法的时间复杂度是 $O(n^2)$ 的.


多项式的点值表示法

点值表示法是 DFT 算法的核心.

我们在 $\mathbb{C}$ 上任意取一组不同的 $x_i$,共 $n$ 个,记为向量 $(x_0,x_1,\dots,x_{n-1})$,这里的 $x_i$ 称为插值节点.

将每个 $x_i$ 分别代入 $A(x)$,得到一组共 $n$ 个运算结果 $y_i$,记为向量 $(y_0,y_1,\dots,y_{n-1})$.

也就是说,求

向量 $(y_0,y_1,\dots,y_{n-1})$ 称为多项式 $A(x)$ 的点值表示.


下面我们证明,一个 $n$ 项多项式($n-1$ 次)在 $n$ 个不同点的取值唯一确定该多项式.

考虑反证法.

假设两个不同的多项式 $A(x)$ 和 $B(x)$,它们满足,对 $\forall i\in[0,n-1]$,总有 $A(x_i)=B(x_i)$.

令 $T(x)=A(x)-B(x)$,则对 $\forall i\in[0,n-1]$,都有 $T(x)=0$.

也就是说,$n-1$ 次多项式 $T(x)$ 有 $n$ 个根,和代数基本定理(一个 $n-1$ 次多项式在复数域上有且仅有 $n-1$ 个根)矛盾,假设不成立,证毕.


通过系数表示求点值表示,复杂度为 $O(n^2)$.

通过点值表示求系数表示,可以使用插值算法,其朴素形式复杂度为 $O(n^2)$.

若已知 $A(x)$ 和 $B(x)$ 的点值表示 $(a_0,a_1,\dots,a_{n-1})$ 和 $(b_0,b_1,\dots,b_{n-1})$,那么多项式乘积 $C(x)=A(x)\cdot B(x)$ 可以简单地表示为

设 $A(x)$ 的点值表示有 $r$ 项,$B(x)$ 的有 $s$ 项,为了保证所得能确定唯一的 $C(x)$,我们需要在 $A$、$B$ 中至少取 $r+s$ 个插值节点.

倘若能快速地转换点值表示法系数表示法,求 $C(x)$ 的过程就会被大大加速.


确定了方向,让我们来捋一捋 DFT 求算 $C(x)$ 的思路,以便于思考.

  1. 输入 $A(x)$ 和 $B(x)$ 的系数表示;
  2. 将其分别转换为点值表达式;
  3. 通过这两个点值表达式计算 $C(x)$ 的点值表达式;
  4. 将 $C(x)$ 的点值表示转化为系数表示,即为结果.

单位复数根模型

互联网上有很多相关的资料,这里只提一下后面要用的引理.

  1. $\omega_n^k=\cos(2\pi\cdot\frac{k}{n})+i\sin(2\pi\cdot\frac{k}{n})$;
  2. $\omega_n^k=\omega_n^{k\mod n}$;
  3. (折半引理)$\omega^{2k}_{2n}=\omega^k_n$;
  4. (消去引理)$\omega^{k+n/2}_n=-\omega^k_n$.

确保理解前置知识过后,我们进入正题.


离散傅里叶变换

考虑一个 $n$ 项多项式 $A(x)$,其中 $n=2^x$,$x\in\mathbb{N}^*$,不足 $2^x$ 的补系数为 $0$ 的项.

倘若取 $n$ 次单位根的 $0\sim n-1$ 次幂,分别代入 $A(x)$,得到一个点值向量 $(A(\omega_n^0),A(\omega_n^1),\dots,A(\omega_n^{n-1}))$.

对给定多项式进行这样操作的过程,称为离散傅里叶变换(Discrete Fourier Transform, DFT).

显然,如果直接代入,DFT 的朴素实现仍然很低效.

考虑对其优化.


对于任意多项式 $A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}$,将每一项按照幂次的奇偶分组:

那么,有

对于 $\forall k\in[0,n-1]$,$k\in\mathbb{R}$,

考虑如何简化右式.


info

注意,请牢记 $n$ 是 $2$ 的整数次幂,即 $n=2^x$.

对于 $k\in[0,\frac{n}{2}-1]$ 的部分($\frac{n}{2}$ 是整数):

对于 $k+\frac{n}{2}\in[\frac{n}{2},n-1]$ 的部分:

注意到:

则有:

$k$ 与 $k+\frac{n}{2}$ 取遍了 $0\sim n-1$ 这 $n$ 个整数,因此能通过这 $n$ 个点值反解系数.


把这两个式子写在一起:

只要知道 $\omega_{n/2}^0,\omega_{n/2}^1,\dots,\omega_{n/2}^{n/2-1}$,转化为子问题,就可以在 $O(n\log n)$ 的时间内求得 $A(x)$.

这个过程就是 DFT.


离散傅里叶反变换

将用单位根求出的点值表示的多项式转化为系数表示,这个过程即为离散傅里叶反变换(Inverse Discrete Fourier Transform).

也就是说,要通过插值为单位根的 $n$ 维点值向量 $(A(\omega_n^0),A(\omega_n^1),\dots,A(\omega_n^{n-1}))$ 推出 $n$ 维系数向量 $(a_0,a_1,\dots,a_{n-1})$.

可以把离散傅里叶反变换理解为离散傅里叶变换的逆操作.

假设 $(a_0,a_1,\dots,a_{n-1})$ 进行一次 DFT 的结果是 $(d_0,d_1,\dots,d_{n-1})$.

构造辅助多项式:

在 $F(x)$ 中,$(d_0,d_1,\dots,d_{n-1})$ 是系数向量.


取插值节点 $x_k=\omega_n^{-k}$($k\in[0,n-1]$),将其代入 $F(x)$,得到另一组点值向量 $(c_0,c_1,\dots,c_{n-1})$.

也就是说:

将 $d_i$ 还原,并化简 $c_k$:


考虑这个式子的右半部分,令:

则有 $c_k=\sum\limits_{j=0}^{n-1}a_j\cdot S(j-k)$,

注意到 $\omega_n^\delta\neq1$ 时,这是一个等比数列.

分类讨论 $\delta$.


如果 $\omega_n^\delta=1$,有 $\delta=0$,则:

$j=k=0$,即 $j=k$ 时,$S(\delta)=n$.


如果 $\omega_n^\delta\neq1$,有 $\delta\neq0$,则 $S(\delta)$ 是公比为 $\omega_n^\delta$ 的等比数列.

根据等比数列求和公式:

也就是说,$j\neq k$ 时,$S(\delta)=0$.


故 $S(j-k)$ 是这样一个分段函数:

将其代入原式:

得到:

这里的 $a_k$ 就是原式的系数.

(忍不住说一声,这个式子太简洁了,向研究傅里叶变换的前辈们致以崇高的敬意!)


由此我们可以得到以下步骤:

  1. 对多项式 $A(x)$,将其用插值节点 $(\omega_n^0,\omega_n^1,\dots,\omega_n^{n-1})$ 进行一次 DFT,得到 $(d_0,d_1,\dots,d_{n-1})$.

  2. 我们再将 $(\omega_n^0,\omega_n^1,\dots,\omega_n^{n-1})$ 作为系数,用插值节点 $(\omega_n^{-0},\omega_n^{-1},\dots,\omega_n^{-(n-1)})$ 再进行一次 DFT,得到 $(c_0,c_1,\dots,c_{n-1})$.

  3. 将 $(c_0,c_1,\dots,c_{n-1})$ 每一项除以 $n$,得到的向量就是我们要求的 $(a_0,a_1,\dots,a_{n-1})$.

这个过程就是 IDFT.


快速傅里叶变换

略讲一下 FFT 的原理.

  1. 给定多项式 $A(x)$ 和 $B(x)$,得到其系数表示.
  2. 对 $A(x)$ 和 $B(x)$ 分别进行 DFT,得到两个点值向量.
  3. 对这两个点值向量每一项相乘,得到一个点值向量.
  4. 对这个点值向量进行 IDFT,即得到所求 $C(x)$ 的系数表示.

以后可能会在另一篇文章补上 FFT 的代码实现,敬请期待.


Reference

本文作者:Xecades

本文链接:https://blog.xecades.xyz/articles/DFT-IDFT/

文章默认使用 CC BY-NC-SA 4.0 协议进行许可,使用时请注意遵守协议。

评论