快速傅里叶变换

傅里叶变换、离散傅里叶变换与快速傅里叶变换

The Fourier Transform is a tool that breaks a waveform (a function or signal) into an alternate representation, characterized by the sine and cosine functions of varying frequencies. The Fourier Transform shows that any waveform can be re-written as the sum of sinusoidals.

傅里叶变换(Fourier transform，FT)将函数(在应用上称：信号，波形）分解为不同特征的正弦函数的和，从而将其由时域转换到频域。逆傅里叶变换则反向操作，因其基本思想首先由法国学者约瑟夫·傅里叶系统地提出，所以以其名字来命名以示纪念。

一般情况下，若“傅里叶变换”一词不加任何限定语，则指的是“连续傅里叶变换”（连续函数的傅里叶变换）。在工程应用中，得益于数字技术的应用，绝大多数傅里叶变换的应用都是采用离散傅里叶变换（DFT），更确切的说，是它的快速算法快速傅里叶变换（FFT，Fast Fourier Transform)。离散傅里叶变换就是对连续傅里叶变换离散化, 等价于对信号进行采样，将其变为有限序列的傅里叶变换。

The Discrete Fourier Transform (DFT) is the equivalent of the continuous Fourier Transform for signals known only at N instants separated by sample times T (i.e. a finite sequence of data) ^4

快速傅里叶变换（Fast Fourier Transform，FFT）是一种可在 O(nlogn) 时间内完成的离散傅里叶变换（Discrete Fourier transform，DFT）算法。而从 computer science 的视角，FFT 主要是被用来将 O(n^2)复杂度多项式乘法加速为 O(nlogn) ，从而快速计算多项式乘法。

FFT

naive polynomial multiplication

考虑朴素的多项式乘积算法，考虑到两个多项式 A(x), B(x) 的乘积 C(x)，假设 A(x)的项数为 n, 其系数构成 n 维向量 $(a_0, a_1,\cdots,a_n)$，B(x) 的项数为 m，对应 m 维向量 $(b_0, b_1,\cdots,b_m)$。则 C(x) 的系数构成 n+m-1 维的向量，其算法如下：

1
2
3

for i in range(len(a)):
    for j in range(len(b)):
        c[i + j] += a[i] * b[j]

多项式乘法本质上是加法卷积，而 FFT 可以加速卷积操作。卷积是信号处理重要操作，公式如下：

\[ (f*g)[k] = \sum_{i=\infty}^{\infty} f[k]g[k-i] \]

则多项式乘法可以表示为下面公式。也就是说要计算出 C(x)的 $x^k$ 项的系数，需要 A(x) 的 $x^i$ 项的系数和 B(x)的 $x^{k-i}$ 项的系数乘积。

\[ C[k] = \sum_{k=\infty}^{\infty} A[k]B[k-i] \]

对应代码：

def bound(idx, arr): return arr[idx] if idx in range(len(arr)) else 0
for k in range(len(a)+len(b)-1):
    for i in range(k+1):
        c[k] += bound(i, a) * bound(k-i, b)

注意补 0 操作：对于索引超出 [0,m] 的 f(x) 的多项式系数以及索引超出 [0，n] 的 g(x) 多项式系数设置为 0

也可以写成：$C[k] = \sum_{i+j=k} A[i]B[j]$

形如 $C[k] = \sum_{i \bigoplus j=k} A[i]B[j]$ 的式子称为卷积，其中 $\bigoplus$ 为某种运算。

polynomial representations

对于 degree = n-1 多项式 A(x)，

Coefficient representation: 所有项的系数组成的系数向量$(a_0, a_1,a_2,...,a_{n-1})$ 唯一确定多项式:$A(x) = \sum_{i=0}^{n-1}a_ix^i$
Point-value representation: 一组互不相同的$(x_0, x_1,x_2, \cdots,x_{n})$ 代入 A(x),得到 n 个取值 $(y_0,y_1, \cdots, y_n)$, 其中 $y_i = \sum_{j=0}^{n}a_j x_i$，可以通过拉格朗日插值方法唯一确定多项式。

定理: 一个 n-1 次多项式在 n 个不同点的取值唯一确定了该多项式。

已知一组插值点 $(x_0, x_1,x_2, \cdots,x_{n})$中 A(x) 与 B(x)（假设多项式阶数相同，没有的视为 0，这也就是为什么补 0 的原因）对应的点值向量分别为 $(y_{a0},y_{a1}, \cdots, y_{an})$ 和 $(y_{b0},y_{b1}, \cdots, y_{bn})$, 则 C(x)=A(x)B(x)的点值表达式可以在 O(n) 的时间求出来，为 $(y_{a0}y_{b0},y_{a1}y_{b1}, \cdots, y_{an}y_{bn})$. 因为 C(x) 的项数为二者之和，而 A(x), B(x) 分别有 n,m 项，所以我们代入的差值节点至少有 n+m 个。

在多项式的点值表示法下多项式乘法可以在 O(n) 时间内求出，远快于系数表示法的 O(n^2)暴力卷积。但是如果按照定义根据系数表示求一个多项式的点值表示，时间复杂度为 O(n^2). 已知多项式的点值表示，求其系数表示，使用朴素的插值算法时间复杂度为 O(n^2)。因此，FFT 通过以下思路优化多项式乘法:

polynomial multiplication under Point-value representation
choose root of unity to accelerate the convertion between Point-value representation and class representation.

当然，如果从矩阵角度理解的话，FFT 通过以下思路优化多项式乘法: factor the matrix of polynomial multiplication with the property of root of unity. [^7]

FFT

FFT 算法流程：把系数表达转换为点值表达 -> 点值表达相乘 ->把点值表达转换为系数表达。

“把系数表达转换为点值表达”的算法叫做 DFT， “把点值表达转换为系数表达”的算法叫做 IDFT(DFT 的逆运算)

其实插值点的选择没有规定，只要选择的插值点不同，个数足够就可以进行 DFT。但是二者转换的复杂度就是 O(n^2) 了，和多项式乘法一样没有必要。而傅立叶，一个很懂复数的人，把选择了单位圆的单位根 root of unity 代入了多项式,求点值表达, 利用单位根的特性就实现了 O(nlogn) 的效率。

复数与欧拉函数

Euler's Formula on unit of circle: for any real x , $e^{ix} =cosx+isinx$, where i is the imaginary unit, $i^2= −1$.

根据复数的知识，向量的计算，乘除法比较方便，直接将两个相量的模相乘/除、辐角相加/减即可；而加减法就需要做反复的直角坐标-极坐标的相互转换，计算非常麻烦。

具有实部和虚部的复数做乘法等运算实际上很不方便，但是借助于欧拉公式，虽然加减法计算依然很麻烦，但是复数的乘除法变得简单起来: 复数相乘时,模长相乘,幅角相加：

设单位圆上复数点 $x_1,x_2$ 相乘，

\[ (cosx_1+isinx_1)(cosx_2+isinx_2) = e^{ix_1}*e^{ix_2} = e^{i(x_1+x_2)}= cos(x_1+x_2)+i sin(x_1+x_2) \]

将单位圆的角度均分为 n 份,并选择对应的从原点向外辐射的半径做向量，其中幅角为正且最小的向量称为 n 次单位向量，记为 $\omega_n^1$。 $w_n^k = e^{\frac{2 \pi k}{n} i}$.

因此，单位向量有以下运算性质:

$\omega_i = \omega^i$, where $\omega$ is a generator.
$\omega_n^k = \omega_n^{k-1}*w^1_n \quad (2 \le k \le n)$
$\omega_{2n}^{2k} = w^k_n$，从几何理解分 2n 份和分 n 份，然后比例一样，辐角一样。
$w^{k+\frac{n}{2}} = -w^k_n$, 多转了半圈，方向相反。

注意: root of unity 的阶选择 2^m，2 的次幂。这与子群的阶和群的阶的关系有关

Discrete Fourier Transform

FFT 利用单位根的性质，将多项式分成奇偶，采用分治方式操作。

DFT: 考虑 n 项（deg=n-1）多项式 A(x), 其系数向量为 $(a_0, a_1, \cdots, a_{n-1})$, 将 n 次单位根的 0~n-1 次幂分别代入 A(x) 得到其点值向量 $(A(w_n^0), A(w_n^1), \cdots, A(w_n^{n-1}))$. 注意: 这里假设 n 是 2 的幂次。

对于 $A(x) = a_0 + a_1x^1 + \cdots + a_{n-1}x^{n-1}$，将其奇偶分开，得到:

$A(x) = a_0 + a_2x^2 + a_4x^4 + \cdots + a_{n-2}x^{n-2} + x(a_1 + a_3x^2 + a_5x^4 + \cdots + a_{n-1}x^{n-2})$

$AL(x) = a_0 + a_2x^2 + a_4x^4 + \cdots + a_{n-2}x^{\frac{n-2}{2}}$

$AR(x) = a_1 + a_3x^2 + a_5x^2 + \cdots + a_{n-1}x^{\frac{n-2}{2}})$

$A(x) = AL(x^2) + xAR(x^2)$

然后将单位根分别带入 AL(x^2) 和 AR(x^2)

设 $0 \le k \le \frac{n}{2}-1 \quad k \in Z$ , 根据前面提到的单位根的特性，$A(w_n^k) = AL(w_n^{2k}) + w_n^{k} AR(x^{2k}) = AL(w_{n/2}^{k}) + w_n^{k} AR({n/2}^{k})$

当 $ k+ n-1$,

\[ A(w_n^{k+\frac{n}{2}}) = AL(w_{n}^{2k+n}) + w_n^{k+\frac{n}{2}} AR(w_{n}^{2k+n}) = AL(w_{n/2}^{k}) - w_n^{k} AR(w_{n/2}^{k}) \]

注意 k 与 k+n/2 取遍了[0，n-1]中的 n 个整数，保证了可以由这 n 个点值反推解出系数.

因此，只需要知道$AL(x)$ 和 $AR(x)$ 分别在 $w_{n/2}^0,w_{n/2}^1,\cdots,w_{n/2}^{n/2-1}$ 的取值，就可以在 O(n) 时间得出 A(x) 的取值，而二者是 A(x)的一般规模，可转化为子问题递归求解，时间复杂度: T(n) = 2T(n/2)+O(n) = O(nlogn).

Inverse Discrete Fourier Transform

使用快速傅里叶变换将点值表示的多项式转化为系数表示，这个过程叫做离散傅里叶反变换（Inverse Discrete Fourier Transform）。

References

推荐 Introduction to Algorithms，讲解思路很好，从 evaluation on point-value polynomial 到 FFT,通过 Vandermonde matrix 求逆来讲解 inverse FFT.