引子

​ 数字信号处理FFT讲的和什么一样，完全听不懂还好这个有翻译的3Blue1Brown讲解FFT，不然就gg了，由于视频讲的太妙了，让我不总结一下不行。

Step 1

​ 如果让你求两个多项式相乘后的多项式即求$C(x):C(x)=A(x)*B(x)$，你会想到什么方法？On的暴力？没错，但是我们能否改进呢：

引理1：

  • 最少已知一个多项式的m+1个点，则可以唯一确定这个多项式（m为多项式的阶数，下面假设$n=2^{\lceil{log(m+1)}\rceil}$​​​​）

由于多项式可写成矩阵相乘的形式，则可通过范德蒙行列式和克莱默法得到上述结论。

Step 2

现在问题则变成了：

​ 已知两个n阶多项式的n+1个点，求这两个多项式的乘积后的多项式

现在暴力则变成了O(n)，对吗？但是还有些问题，比如如何找这n+1个点及已知n+1个点如何转换为多项式。后者我们先按下不表，先解决找点的问题。

如果随便找n+1个点，每个点x计算y(x)也是O(n)总复杂度还是O(n^2)没有变，我们需要更高效的找点算法：

假设函数$y=x^2+1$，如果让你找8个点是不是会倾向于找对称点呢？比如寻找$[(1,2),(2,5),…(4,17)]$这样你就知道了$[(-1,2),..(-4,17)]$寻找8个点就变成了找4个点的问题，规模小了一半！类比于奇函数也有这种性质。

我们知道任何函数都可以表示偶函数和奇函数的和，多项式也不例外。

假设有多项式:

可以变形为：

为方便理解，我变量替换了一下，t=x^2，可以看出这里$P_e$和$P_t$​​都是偶函数，由于：

故满足上述条件，此变换可以将寻找的点数减半。我们只需求出$P_e(t)$和$P_o(t)$两个多项式即可。

看到这里，显然可以发现这是一个递归过程，因为$P_e(t)$和$P_o(t)$明显可以在像$P(x)$一样拆分。

但是遗憾的是$t^\prime$​​ 在实数域取负值，也就是说无法利用O(1)对称，似乎陷入了瓶颈？接着往下看，这就是FFT的由来。

FFT

​ 我们现在需要的是寻找一个$x_1、x_2$有$x_1^2=-x_2^2$，显然实数域不存在，故拓展到负数域，由于可以随便取初始点，方便起见令$t^n=1$，由欧拉公式可知，这样寻找的n个点则是在负数域单位圆周上等间隔分布即第k个点$w_k=e^{\frac{2\pi}{n}*k}$​​，这样无论递归到何时，t都存在对称点的情况，同时推导P(x)的公式也需要跟着改变：

由于

根据第i个点$w^i$求函数值如下：

FFT代码


void FFT(vector<complexx> &A, int tempn,int flag)
{
    // flag作用后面会解释
    if (tempn == 1)
        return;
    int M = tempn/2;
    vector<complexx> A0, A1;
    for (int i = 0; i <= tempn; i+=2)
    {
        A0.push_back(A[i]);
        A1.push_back(A[i+1]);
    }
    FFT(A0, M,flag);
    FFT(A1, M,flag);
    auto W = complexx(cos(1.0 * Pi / M),flag*sin(1.0 * Pi / M));
    auto w = complexx(1.0, 0.0);
    for (int i = 0; i < M; ++i)
    {
        A[i] = A0[i] + w * A1[i];
        A[i + M] = A0[i] - w * A1[i];
        w= w*W;
    }
    A0.clear();
    A1.clear();
}

IFFT

现在还剩最后一个问题，如何将已知的函数值转换回对应的多项式系数呢？

已知多项式可以被以下矩阵乘法表示

则取逆变换有：

由范德蒙行列式的特性求逆可转换为以下形式：

可见IFFT与FFT变动并不大，代码甚至只需要修改W的正负（flag的作用），然后最后除n即可

参考