总觉得不会写 FFT 有愧于天地。

# 傅里叶变换 —— 连续到离散

高等数学的知识告诉我们，周期信号可以被分解成一系列正弦信号之和（当然要满足一定的条件，不过工科一般也碰不到吧），也就是所谓的傅里叶级数：

$$f(t)=\frac{A_0}{2}+\sum_{n=1}^{\infty}A_n\cos(n\Omega t+\varphi_n)$$

这个式子也可以写成复指数的形式：

$$f(t)=\sum_{n=-\infty}^{\infty}F_n e^{jn\Omega t}$$

对于非周期的信号，可以将其视作周期无穷大的周期信号，在这种情况下，傅里叶级数中$\Omega$ 会趋近于无穷小，离散的频谱概念则不再适用。因此需要引入频谱密度的概念：

$$F(j\Omega)=\lim_{T\to\infty}\frac{F_n}{1/T}=\lim_{T\to\infty}F_nT$$

这样，有：

$$F_nT=\int_{-T/2}^{T/2}f(t)e^{-jn\Omega t}dt$$

即

$$f(t)=\sum_{n=-\infty}^{\infty}F_ne^{jn\Omega t}$$

当$T\to\infty$，$\Omega$ 趋近无穷小，记为$d\omega$，则有：

$$F(j\omega)=\lim_{T\to\infty}F_nT=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt$$

反之有：

$$f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(j\omega)e^{j\omega t}d\omega$$

这是傅里叶变换的定义。

不管怎么样，傅里叶变换看起来都还是不错的。不过对于计算机来说却不是 —— 傅里叶变换处理连续的函数，得到连续的频谱。这显然难以通过计算机实现，因此，我们需要将其离散化。

首先，我们将时间轴离散化，将时间轴上的点看作是一个个的采样点，这样，我们就得到了一个离散的序列$f_N(k)$。采用和傅里叶变换类似的定义，有：

$$F(e^{j\theta})=\sum_{k=-\infty}^{\infty}f_N(k)e^{-jk\theta}$$

这也就是离散时间傅里叶变换（DTFT）。然而，这里的$\theta$ 是一个连续的值：$N\to\infty$，$n\frac{2\pi}{N}\to\theta$，这使得$\theta$ 成为了一个连续的变量（也就是频谱依旧是一个连续量）。DTFT 的逆变换为：

$$f(k)=\frac{1}{2\pi}F(e^{j\theta})e^{jk\theta}d\theta$$

这里仍然需要对一段连续的区间进行积分，对于计算机来说依旧是难以计算的。

如果考虑一个周期序列，周期序列对应的频谱是一个离散谱，那么如果能够将需要计算的序列$f_N(k)$ 转换为一个周期序列，就能够得到便于计算机处理的离散谱。

考虑一个有限长序列$f(k)$，其在区间$[0,N-1]$ 中有取值，其余区段取值为 0：

$$f(k)=\begin{cases} f(k),&0\leq k\leq N-1\\ 0,&\text{otherwise} \end{cases}$$

将其延拓到整个时间轴上，有：

$$f_N(k)=\sum_{l=-\infty}^{\infty}f(k+lN),\ l\in\mathbb{Z}$$

衍伸傅里叶级数的定义，得到：

$$F(n)=\sum_{k=0}^{N-1}f_N(k)e^{-j\frac{2\pi}{N}nk}\\ f(k)=\frac{1}{N}\sum_{n=0}^{N-1}F(n)e^{j\frac{2\pi}{N}nk}$$

这便是离散傅里叶变换（DFT）和离散傅里叶逆变换（IDFT）。

我们利用 DFT 的公式可以比较简单的实现离散傅里叶变换：

def dft(x):

    N = len(x)

    n = np.arange(N)

    return np.array([np.sum(x * np.exp(-1j * 2 * np.pi * k * n / N)) for k in range(N)])

def idft(X):

    N = len(X)

    n = np.arange(N)

    return np.array(

        [np.sum(X * np.exp(1j * 2 * np.pi * k * n / N)) / N for k in range(N)]

    )

显然，这里的时间复杂度为$\mathcal{O}(n^2)$，明显是不够好的。我们需要寻找更好的算法。

# 快速傅里叶变换 —— 从多项式乘法开始

在引入快速傅里叶变换的概念之前，先来考虑一个问题 —— 如何快速地计算多项式乘法？

举个例子，现在有两个多项式：

$$A(x) = x^3 + 2x^2 + 3x + 4\\ B(x) = 4x^3 + 3x^2 + 2x + 1$$

我们要如何去计算这两个多项式的乘积$C(x)=A(x)\cdot B(x)$？

显然，直接计算不是不行 —— 我们不妨用更一般的形式来写：

$$A(x)=\sum_{i=0}^{n}a_ix^i\\ B(x)=\sum_{i=0}^{n}b_ix^i$$

那么$C(x)$ 可以表示为

$$C(x)=\sum_{i=0}^{2n}c_ix^i,\\ c_i=\sum_{j=0}^{i}a_jb_{i-j}$$

很不幸，这也是个$\mathcal{O}(n^2)$ 的算法。我们也需要寻找更好的算法。

但是，回头看一下上面的式子 —— 有什么特殊的地方？乘积多项式$C(x)$ 的系数$c_i$

$$c_i=\sum_{j=0}^{i}a_jb_{i-j}$$

正好是多项式$A(x)$ 和$B(x)$ 系数的卷积！

卷积这件事会带来很多有趣的性质，卷积定理告诉我们 —— 时域相乘等于频域卷积，时域卷积等于频域相乘，这也让我们回到了傅里叶变换的问题上来。

不过暂且先抛开傅里叶变换，我们先把目光放到多项式上来。

# 系数表示和点值表示

我们前面采用了多项式的系数来表示它，不过这不算是唯一方法 —— 我们可以通过多项式在一些点上的取值来表示它。

为了后面介绍的方便，我们引入两个术语 ——Evaluation（估计）和 Interpolation（插值）。Evaluation 是指通过多项式的系数来计算多项式在某些点上的取值，而 Interpolation 则是指通过多项式在某些点上的取值来计算多项式的系数。

中学的知识就已经说明了，一个$n$ 次多项式可以通过$n+1$ 个点来唯一确定。因此，我们可以通过多项式在$n+1$ 个点上的取值来表示它：

$$(x_i,y_i),\ i=0,1,\cdots,n\rightarrow A(x)$$

那么利用这种表示，我们也可以计算多项式的乘法：对于$A(x)$ 上的一点$(x_i,y_{a,i})$ 和$B(x)$ 上的一点$(x_i,y_{b,i})$，我们可以计算出$C(x)$ 上的一点$(x_i,y_{a,i}y_{b,i})$。那么计算$2N$ 个点，我们便可以通过拉格朗日插值法确定$C(x)$。

不过，这样的计算量依旧是$\mathcal{O}(n^2)$ 的。

可以优化吗？

# 利用奇偶函数的性质

很中学地，我们注意到，奇函数和偶函数有着这样的性质：

$$f_{odd}(x)=-f_{odd}(-x)\\ f_{even}(x)=f_{even}(-x)$$

对于奇函数和偶函数，我们只需要计算$f(x)$ 就能得到$f(-x)$ 的值 —— 这让我们的计算量减少了一半。

而每个多项式都可以分解为两个部分的和：

$$\begin{aligned} f(x)&=\sum_{i=0}^{n}a_ix^i=\underbrace{\sum_{i=0}^{n/2}a_{2i}x^{2i}}_{P_{even}(x^2)}+x\underbrace{\sum_{i=0}^{n/2}a_{2i+1}x^{2i}}_{P_{odd}(x^2)}\\ &=P_{even}(x^2)+xP_{odd}(x^2) \end{aligned}$$

（这里的表述有点小问题，但是不影响意思的表达……）

那么我们计算得到$P_{even}(x_i^2)$ 和$P_{odd}(x_i^2)$ 之后，我们就可以得到$f(x_i)$ 和$f(-x_i)$ 的值了：

$$\begin{aligned} f(x_i)=P_{even}(x_i^2)+x_iP_{odd}(x_i^2)\\ f(-x_i)=P_{even}(x_i^2)-x_iP_{odd}(x_i^2) \end{aligned}$$

这样我们把一个$n$ 次多项式的计算量减少到了两个$\frac{n}{2}$ 次多项式的计算量，而且看起来是能够递归计算的。

能吗？

# 复数？

似乎不能。

还是举个例子。函数$f(x)=x^3+x^2-x+1$，我们需要选择四个点才能够唯一确定它。为了尽可能利用上面的性质，我们选择$x_1,-x_1,x_2,-x_2$ 这四个$x$ 对应的点。

那我们便需要计算$P_{even}(x^2)=x^2+1$ 和$P_{odd}(x^2)=x^2-1$ 两个函数在$x_1,x_2$ 处的取值。拿$P_{even}(x^2)$ 为例，为了尽可能利用递归的性质，$x_2$ 应该等于$-x_1$。

但是，$P_{even}(x^2)$ 传入的参数是$x^2$，而$x_2=-x_1$，这样就会导致$P_{even}(x_2^2)=P_{even}(x_1^2)$，这样就无法利用递归的性质了。

解决办法很简单 —— 复数。

如果我们需要计算的多项式次数为$d$，我们选择$n=2^{\lceil \log_2d\rceil}$（便于二分递归），我们要选的点便是$x^n=1$ 的$n$ 个根：

$$\omega_i=e^{-2j\pi i/n}$$

在这些根中，$\omega_i$ 和$\omega_{n/2+i}$ 是相反的：

$$\omega_{n/2+i}=e^{-2j\pi (n/2+i)/n}=e^{-j\pi}e^{-2j\pi i i/n}=-\omega_i$$

这样，可以保证$P_{even}(x^2)$ 和$P_{odd}(x^2)$ 的参数可以是一对相反数，便可以利用递归的性质了。

我们可以通过代码来实现这个过程：

def evaluation(coeffs):

    n = len(coeffs)

    if n == 1:

        # 对于 0 次也就是常数，其值是恒定的

        return coeffs

    even, odd = coeffs[::2], coeffs[1::2]

    # 分别递归处理偶次项和奇次项

    even_eval, odd_eval = evaluation(even), evaluation(odd)

    omega = np.exp(-2j * np.pi / n * np.arange(n // 2))

    # 在进行 Evaluation 的点中，\omega_i 和 \omega_{n/2+i} 是相反的

    return np.concatenate(

        [even_eval + omega * odd_eval, even_eval - omega * odd_eval]

    )

不难发现这是一个$\mathcal{O}(n\log n)$ 的过程。

# Evaluation 和 DFT

至少，我们现在获得了一种快速计算多项式值的方式 —— 那么这和这篇文章的标题有什么关系？

记$\omega=e^{-2j\pi/n}$，我们现在得到了一系列的值：

$$F=[f(\omega^0),f(\omega^1),\cdots,f(\omega^{n-1})]^T$$

如果我们使用传统的方法计算，这些值从何而来？记

$$f(x)=\sum_{i=0}^{n-1}a_ix^i$$

我们可以把$F$ 的计算过程展开成一个矩阵乘法：

$$\begin{bmatrix} f(\omega^0)\\ f(\omega^1)\\ \vdots\\ f(\omega^{n-1}) \end{bmatrix} = \underbrace{\begin{bmatrix} 1&1&1&\cdots&1\\ 1&\omega&\omega^2&\cdots&\omega^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega^{n-1}&\omega^{2(n-1)}&\cdots&\omega^{(n-1)^2} \end{bmatrix}}_{M} \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix}$$

回过头看 DFT 的表述：

$$\begin{aligned} F(n)=&\sum_{k=0}^{N-1}f_N(k)e^{-j\frac{2\pi}{N}nk} =&\sum_{k=0}^{N-1}f_N(k)\omega^{nk} \end{aligned}$$

记$f_N=[a_0,a_1,...,a_{n-1}]^T$，将 DFT 的结果序列记作$F=[F(0),F(1),...,F(n-1)]^T$，如果将$F$ 的计算过程展开成一个矩阵乘法：

$$F=M_{DFT}f_N$$

这里的$M_{DFT}$ 和上面的$M$ 是一致的。

也就是说，这个对多项式的 Evaluation 过程，其实就是 DFT 的过程 —— 我们对多项式的系数进行了 DFT，得到了多项式在$\omega^0,\omega^1,...,\omega^{n-1}$ 这些点上的取值（相当于系数的频谱）。

更加值得庆祝的是，我们现在已经有了一个$\mathcal{O}(n\log n)$ 的 DFT 算法了，这就是我们所说的快速傅里叶变换（FFT）。

# Interpolation 与 IDFT

既然 Evaluation 和 DFT 是一致的，而 Interpolation 和 Evaluation 是相反的过程，那 Interpolation 和 IDFT 是否是一致的？

显然是一致的。依据 IDFT 的定义，我们能够写出：

$$\begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix} = \frac{1}{n} \underbrace{\begin{bmatrix} 1&1&1&\cdots&1\\ 1&\omega^{-1}&\omega^{-2}&\cdots&\omega^{-(n-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega^{-(n-1)}&\omega^{-2(n-1)}&\cdots&\omega^{-(n-1)^2} \end{bmatrix}}_{M^{-1}} \begin{bmatrix} f(\omega^0)\\ f(\omega^1)\\ \vdots\\ f(\omega^{n-1}) \end{bmatrix}$$

不难验证$M$ 和$M^{-1}$ 互为逆矩阵，上式中左侧为多项式系数，右侧为逆 DFT 矩阵和多项式在$\omega^0,\omega^1,...,\omega^{n-1}$ 这些点上的取值，这个过程也就是插值。

更妙的是，和 Evaluation 相比，两者只是相差一个$n$ 的倍数和$\omega$ 的符号，从程序上来说，几乎没有区别。

def ifft(x):

    n = len(x)

    if n == 1:

        return x

    even, odd = coeffs[::2], coeffs[1::2]

    # 分别递归处理偶次项和奇次项

    even_eval, odd_eval = ifft(even), ifft(odd)

    omega = np.exp(2j * np.pi / n * np.arange(n // 2))

    return np.concatenate(

        [even_eval + omega * odd_eval, even_eval - omega * odd_eval]

    ) / n

# 回到多项式乘法问题

回到原先的多项式乘法……

我们现在已经可以通过 FFT 来计算多项式在$\omega^0,\omega^1,...,\omega^{n-1}$ 这些点上的取值了，通过两次 FFT，我们可以计算得到多项式$A(x)$ 和$B(x)$ 在这些采样点上的取值，然后相乘，再通过一次 IFFT，我们就可以得到多项式$C(x)$ 的系数。

整个过程的时间复杂度还是$\mathcal{O}(n\log n)$。

# End

嗯，终于学会 FFT 了，可喜可贺，可喜可贺。