Fourier Transform，如此美妙！

本文概述

1. 函数可以做坐标轴，像是 x, y, z 轴那样

2. 如果坐标轴正交，还可以做投影之类的事

3. 信号在函数空间中投影

4. 傅里叶变换 & 离散傅里叶变换

5. fft 和 ntt 算法

6. 用傅里叶圆画画！

7. 傅里叶变换其他应用

前言

信号，无处不在，比如音乐频谱：

简单来说，傅里叶变换就是将时域信号函数 f(t) 转化成频域信号函数 F(ω)，而在频域中，不仅信号的频率成分一目了然，而且还有一些美妙的数学性质我们可以利用，比如时域中进行卷积等价于频域中逐点相乘，后面会详细介绍

傅里叶级数告诉我们，信号可以展开成正弦曲线：

或者换成等价的三角形式，也是大部分高数教材上的样子：

我们可以从坐标分解的角度去理解傅里叶级数，其相当于是把无数个两两正交的正余弦函数作为坐标轴，然后让信号映射到该空间中，信号向量在某一维度上投影的大小就是其在该维度的坐标，类似我们熟悉的 x-y-z 正交坐标系。这些函数坐标轴张成了一个高维空间，存在于这个高维空间中的向量就是信号，维度越高，即不同频的正余弦函数越多，信号在此空间中的投影就越精细，越接近原信号

函数也可以作坐标轴？

在 ℝ³ 空间中，若选择一组基底 (e₁, e₂, e₃)，其中

则任意向量 v(x, y, z) ∈ ℝ³ 都可以唯一表示为：

类似的，若选择一组函数作为基底，比如 (1, x, x²)

则这 3 条抽象的函数坐标轴能张成一个 V = {a + bx + cx² ∣ a, b, c ∈ ℝ} 的多项式空间，此空间包含了所有最高项次数不超过 2 的多项式

任意 f ∈ V 都可以表示为：

f 在此坐标系中的坐标就是 (a, b, c)

正交性

平日我们生活的空间中，最多只能找出三个两两正交的坐标轴，但是函数不同，可以存在无数个两两正交的函数坐标轴。函数中定义的正交是指，在确定区间上内积为零，即满足：⟨u₁, u₂⟩ = ∫_a^bu₁(x) ⋅ u₂(x) dx = 0，也称为：函数 u₁(x)、u₂(x) 在该区间正交

上述的 e₁, e₂, e₃ 两两正交，是一组正交基底，即：

一组正交函数基底比如 ，其在 [ − 1, 1] 上两两正交：

三角正交基，其在傅里叶级数中被应用：

不同频率的正余弦函数在 x ∈ [0, 2π] 正交：

cos  可以看成 sin  偏移 π/2 相位，同理，cos (mx) 和 cos (nx)、sin (mx) 和 sin (nx)、sin (mx) 和 cos (nx) 也在此区间正交

复指数正交基，三角正交基的一种等价正交基，在复傅里叶级数中应用：

由欧拉公式 e^ix = cos x + isin x 可知，三角函数和 e 有天然的联系：

复指数函数之间同样正交：

投影与正交基展开

正交形式便于我们计算和使用投影

在普通向量空间中，选定一组基底后，空间中任意向量都可表示成基向量的线性组合，函数也类似，基本上可以理解成一个无限维向量 (f(x₀), f(x₁), …, f(x_n), …)，向量或者是函数 f 在某个基函数 u_k 上的投影向量为：

其中 c_k 就是 f 在基函数 u_k 上投影的系数

也就是说，对于一组正交基函数，我想知道函数向量在某一维的大小(即该方向的系数)，我只需要计算其在该维度的投影就可以了，相应的，若我知道了每个方向上的系数大小，那我也可以重新拼回原向量

对于一组正交基函数 u₁, u₂, …, u_m 张成的空间 V，函数 f 在此空间上的投影为：

我们给它一个名字叫作：正交基展开，即将一个函数或向量唯一的用一组正交基函数或基向量表示

这种展开理论上可以无损（但是大部分时候不需要无损

在 x, y, z 轴张成的空间 V 即 ℝ³ 中，若向量 v 本身属于 ℝ³，则

而当向量 v 是一个比当前空间维度更高维的向量时，可以预想，若 v 在此空间强行投影，则会损失在空间 V 之外的部分，因为高维信息根本无法被表达，只能得到一个有损向量：

函数也同理

如你所料，傅里叶级数就是正交基展开，其思路就是将复杂函数转换成简单正交基函数，也就是正余弦函数

另外，傅里叶变换也是正交基展开，其使用傅里叶基，而此基性质相当微妙，让我们得以快速变换

Fourier Transform

来到正题，有了上述前置知识，傅里叶变换就相当好理解了

其非常之精妙，式子极其简洁且优雅：

两个式子分别是正变换(FT)：将函数由时域 f(t) 转为频域 F(ω)，和逆变换(IFT)：将函数由频域转回时域，其中 u_ω(t) 为基函数

如你所料也如同上文所说，傅里叶变换其实就是将 f(t) 投影到正交基函数上，逆变换也就是使用基函数和其对应的系数重新拼回原函数

有个问题你可能已经注意到，为什么前面的系数是 ？据上方投影章节所述，，这里的系数不该是 吗？非也，此处变量 t 定义在整个实数轴，积分上下限都是无穷，致平方长度非有限，因此普通正交定义已不再适用，内积 ⟨u_ω, u_ω⟩ 不可直接计算。为解决此问题，需使用更高级的数学工具：广义正交定义，由此可以推出 ⟨f, u_ω⟩ = 2πF(ω)，即 ，在此不展开

Discrete Fourier Transform

在实际应用中，对着无数点疯狂积分显然不现实，计算机只能处理有限采样点

选点很 easy，就假设从 f(t) 中均匀选出 n 个点，那之后的问题是，要选多少个正交基函数 u_k(t)？注意到选出的这些点可以看作一个 n 维向量，而 n 维向量需用 n 个线性无关的基底才能唯一表示，因此需要 n 个，而每个基函数也得选 n 个点作为 n 维基向量

于是就有，离散傅里叶变换(DFT)：

其中，正交基向量：

单位根 ω_n = e^2πi/n，傅里叶矩阵为：

整体系数

即

傅里叶矩阵求逆非常之便捷，观察傅里叶矩阵(或其共轭)发现，其和正交矩阵只差一个系数(严格来说复数域中应叫酉矩阵)

而复数域中的正交矩阵有一个很重要的性质 M^− 1 = M^H，即正交矩阵的逆矩阵等于其共轭转置

并且傅里叶矩阵是对称矩阵，满足 W^T = W

因此，W 的逆矩阵为：

也就有，离散傅里叶逆变换(IDFT)：

矩阵 W 的第 k 行是第 k 个傅里叶基，信号 f(t) 在 u_k 上的系数就是

进而

进行 DFT 时，计算每个系数 F_k 需要 O(n)，IDFT 还原时，计算每个 f(t_j) 需要 O(n)，所以时间复杂度均为 O(n²)

Fast Fourier Transform

前置知识：单位根

将复平面的单位圆等分成 N 份，圆上的点所代表的复数就是 N 个 N 次单位根

其恰好是 ω_N = e^2πi/N 的 N 个不同的幂，也组成了 z^N = 1 的所有解：

单位复根有一些重要性质：

折半性

消去性

对称性

可以从视频中获得更直观的理解 FFT(快速傅里叶变换)：优美的分治算法

前文为了便于理解，系数在前，而系数在哪并不会影响傅里叶变换的根本性质，为了便于计算卷积和多项式乘法，下文统一采用以下定义

正变换

逆变换

观察式子

可以发现每个分量都能看作一个多项式 A(x)

于是就有：

也就是说，等价于求多项式 A(x) 在这 N 个点处的值

我们将 A(x) 按奇偶拆开成两个多项式，偶数项组成的多项式记为 A^[0](x)，另一个多项式记为 A^[1](x)

则原多项式可以写成：

将 ω_N^− k 代入，得到：

利用单位根折半性，上式可以写作：

而利用单位根对称性，我们可以得到

令两个子多项式的 DFT 结果分别为：

则

这一步被称之为蝶形运算

对比 y_k 和 y_k + N/2，发现其结果只是子多项式相差负号，因此只需计算 y₀^[0], y₁^[0], …, y_{N/2 − 1}^[0] 和 y₀^[1], y₁^[1], …, y_{N/2 − 1}^[1] 的值，就能算出所有的 y₀, y₁, …, y_N − 1，也就是，求 N 项式在 N 个点处的值转化成了求两个 N/2 项式在 N/2 个点处的值。原问题被拆成了两个规模折半的子问题，为了便于敲代码，还需将 N 变成 2 的幂次，可直接补充系数为零的高次项，易知这种操作并不会影响最终结果，且添加的项数不会超过 N，对时间复杂度也影响不大

递归的过程形成了一个完全二叉树，发挥人类智慧，发现 A[i] 在叶子中所在的位置恰好是 i 的二进制反转，于是可以迭代求解。每层的时间复杂度都是 O(N)，总时间复杂度 O(Nlog N)，具体看代码

CODE

const double PI = acos(-1);
using cp = complex<double>;

void fft(vector<cp>& a, int mode){
    int n = a.size();
    for(int i=1, j=0; i<n; i++){
        int bit = n >> 1;
        while(j & bit){
            j ^= bit;
            bit >>= 1;
        }
        j ^= bit;
        if(i < j) swap(a[i], a[j]);
    }
    for(int len=2; len<=n; len<<=1){
        cp wlen(cos(2*PI/len), -mode*sin(2*PI/len));
        int half = len >> 1;
        for(int i=0; i<n; i+=len){
            cp w(1, 0);
            for(int j=0; j<half; j++){
                auto x = a[i+j];
                auto y = a[i+j+half] * w;
                a[i+j] = x+y; a[i+j+half] = x-y;
                w *= wlen;
            }
        }
    }
    if(mode == -1) for(auto& x : a) x/=n;
}

Number Theoretic Transform

事实上，u_k(x) 不止一种选法

设 p 为质数，g 为 p 的一个原根，令变换长度 N 满足 N ∣ (p − 1) 且 N = 2^m，还可以定义模 p 意义下的单位根

其同样具有进行 FFT 所需的所有性质

对称性，

折半性，

消去性，

因此也就有了整数版本的傅里叶变换，即 NTT（数论变换）。其解决了普通 FFT 所带来的浮点误差，一般计算也更快

其正变换定义为

其逆变换定义为

其中 N^− 1 表示 N 模 p 的乘法逆元

傅里叶变换与卷积

频域中有个很好的性质，时域中的卷积等价于频域中的逐点相乘，即

而我们知道卷积的计算是 O(n²) 的，而逐点相乘只需 O(n)，因此可以加速卷积，这也为快速计算多项式乘法提供基础

证明：

定义

则

令 u = t − τ，则 t = u + τ，dt = du，则