这是一篇从博客原站迁移过来的文章

原发布时间：2021/2/28
1. 为了适应新站的编辑系统进行了部分改写
2. 对内容进行了部分修改

傅里叶变换可以完成时域与频域的相互转换，在竞赛中常用于大数乘法与多项式乘法。

傅里叶变换长什么样

（连续）傅里叶变换/逆变换公式如下：

$$\displaystyle F(\omega)=\mathcal{F}[f(t)]=\int_{-\infty}^{+\infty}f(t)e^{-i\omega t}dt$$

$$\displaystyle f(t)=\mathcal{F}^{-1}[F(\omega)]=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(\omega)e^{i\omega t}d\omega$$

其中$$F(\omega)$$为$$f(t)$$的象函数，$$f(t)$$为$$F(\omega)$$的象原函数

这是一个无穷区间反常积分，众所周知计算机是没法处理连续区间数据的，所以计算机使用的是另一种傅里叶变换——离散傅里叶变换（DFT）：

$$\displaystyle X(k)=\mathrm{DFT}[x(n)]=\sum_{n=0}^{N-1}x(n)e^{-i\frac{2\pi n}{N}k}$$

$$\displaystyle x(n)=\mathrm{IDFT}[X(k)]=\frac{1}{N}\sum_{n=0}^{N-1}X(k)e^{i\frac{2\pi n}{N}k}$$

其中$$X(k)$$为$$x(n)$$的象函数，$$x(n)$$为$$X(k)$$的象原函数

离散傅里叶变换本质上就是从象原函数中等距抽取有限数量的采样点，将积分转换成求和处理。
而我们的快速傅里叶变换（FFT）就是基于离散傅里叶变换算法的优化加速。

从傅里叶级数到傅里叶变换

傅里叶变换源于法国数学家及物理学家——让·巴普蒂斯·约瑟夫·傅里叶的发现：

能将满足一定条件的某个函数表示成正弦/余弦函数或者它们的积分的线性组合

这就是傅里叶级数，也是高数书上所告诉我们的：

$$\displaystyle f(x)=\frac{a_0}{2}+\sum_{n=1}^\infty(a_n\cos nx+b_n\sin nx)$$

$$\displaystyle a_0=\frac{1}{\pi}\int_{-\pi}^\pi f(x)dx$$

$$\displaystyle a_n=\frac{1}{\pi}\int_{-\pi}^\pi f(x)\cos nx dx$$

$$\displaystyle b_n=\frac{1}{\pi}\int_{-\pi}^\pi f(x)\sin nx dx$$

由欧拉公式：

$$\displaystyle\cos x=\frac{e^{ix}+e^{-ix}}{2}$$

$$\displaystyle\sin x=\frac{e^{ix}-e^{-ix}}{2i}$$

我们可以把上面这个式子写成更加简单的复数形式：

$$\displaystyle f(x)=\sum_{n=-\infty}^\infty c_ne^{inx}$$

$$\displaystyle c_n=\frac{a_n-ib_n}{2}=\frac{1}{2\pi}\int_{-\pi}^\pi f(x)e^{-inx}dx$$

以上是对于f(x)是一个周期为$$2\pi$$的函数而言的，对于周期为$$2l$$的函数，只需要进行一次放缩变换即可：

$$\displaystyle f(x)=\sum_{n=-\infty}^\infty c_ne^{i\frac{n\pi x}{l}}$$

$$\displaystyle c_n=\frac{a_n-ib_n}{2}=\frac{1}{2l}\int_{-l}^l f(x)e^{-i\frac{n\pi x}{l}}dx$$

傅里叶变换是连续的，而傅里叶级数仍然是离散的，所以我们要先设法把傅里叶级数给变成连续的（当然变成连续的之后就不能叫傅里叶级数了）。

我们现在重新思考一下变化的$$n$$代表了什么。因为$$n$$会在正弦/余弦函数中以$$x$$的系数出现，也就是作为正弦/余弦函数的角速度，由于角速度与频率之间有关系：

$$\displaystyle f=\frac{\omega}{2\pi}=\frac{n}{2\pi}$$

所以对应到傅里叶级数中，我们就可以知道每个$$a_n$$、$$b_n$$实质上就是将$$f(x)$$用正弦/余弦函数分解后，频率为$$\frac{n}{2\pi}$$的正弦/余弦函数的系数。至于为什么我们要把这些系数用$$c_n$$来表示？其实就是纯属为了省纸。两者在算数意义上是一致的，也仅仅是算数意义上一致，写成复数形式好看而已，不必深究。

为了得到连续的频率，我们只需要把求和换成积分，弃用$$c_n$$（因为已经不是离散了），把它换成一个连续的函数$$F(t)$$即可：

$$\displaystyle f(x)=\int_{-\infty}^\infty F(t)e^{itx}dt$$

$$\displaystyle F(t)=\frac{1}{2\pi}\int_{-\pi}^\pi f(x)e^{-itx}dx$$

走到这一步，我们得到的式子看起来已经和傅里叶变换差不多了，但还是有一点点不一样，比如系数$$\frac{1}{2\pi}$$的位置，但是这个很好理解，因为$$f(x)$$、$$F(t)$$都出现在了对方的积分式里，由积分的性质可知这个系数是可以转移给对方的。

那么只剩下把$$F(t)$$的积分上下限换成无穷了。其实想想就知道这个积分上下限当我们把求和换成积分后，就一定得是无穷，同样是因为积分的性质，两者要互相转换的话，那么积分的上下限是必须要一致的，否则会得到错误的答案。至于为什么无穷求和就行，我认为这应该就是离散和连续的不同吧。（本段纯属推测，我没有学过相关的知识，但是理解还是可以的）
所以我们最终得到了如下式子：

$$\displaystyle F(t)=\int_{-\infty}^\infty f(x)e^{-itx}dx$$

$$\displaystyle f(x)=\frac{1}{2\pi}\int_{-\infty}^\infty F(t)e^{itx}dt$$

这下成功我们从傅里叶级数开始，得到了傅里叶变换。我们把傅里叶变换抽取有限点计算，就得到了离散傅里叶变换。

离散傅里叶变换

我们复习一下离散傅里叶变换的公式：

$$\displaystyle X(k)=\mathrm{DFT}[x(n)]=\sum_{n=0}^{N-1}x(n)e^{-i\frac{2\pi n}{N}k}$$

$$\displaystyle x(n)=\mathrm{IDFT}[X(k)]=\frac{1}{N}\sum_{n=0}^{N-1}X(k)e^{i\frac{2\pi n}{N}k}$$

其中$$X(k)$$为$$x(n)$$的象函数，$$x(n)$$为$$X(k)$$的象原函数

只要你认识公式，我相信你可以直接实现出一份$$O(n^2)$$的DFT算法，只需要对于每个$$X(k)$$跑一次循环即可。
但是在竞赛中，这还是太慢了，所以我们需要用到一些小技巧来升级DFT。

n次单位根

在数学上有这么一类特殊的方程：$$x^n=1$$，由于n次方程有n个复数根，易得：

$$\displaystyle x=e^{i\frac{2\pi k}{n}},k=(0,1,2,...,n-1)$$

形如$$e^{i\frac{2\pi k}{n}}$$的东西就叫n次单位根，我们用$$\varepsilon_k$$来表示，这些根都位于复平面的单位圆上，把单位圆均分成n份。
n次单位根有一些特殊的性质，而这些性质可以用来加速DFT：

性质1

$$\displaystyle \varepsilon_k=\varepsilon_{k\mod n}$$

这个性质可以把式子转换成三角函数，由周期性证得。

性质2

$$\displaystyle \varepsilon_k=(\varepsilon_1)^k$$

这个性质可以把式子转换成自然底数幂证得。

性质3

$$\displaystyle \varepsilon_{k+\frac{n}{2}}=-\varepsilon_k$$

推导：$\begin{equation*}\begin{split}\varepsilon_{k+\frac{n}{2}}&=e^{i\frac{2\pi(k+\frac{n}{2})}{n}}\\&=e^{i\frac{2\pi k+\pi n}{n}}\\&=e^{i\frac{2\pi k}{n}}\cdot e^{i\pi}\\&=-e^{i\frac{2\pi k}{n}}\\&=-\varepsilon_k\end{split}\end{equation*}$

n次本原单位根

这个东西和n次单位根比起来好像就是多了两个字，但是两者的意思却完全不一样。
如果n个n次单位根可以由一个n次单位根的0到(n-1)次幂求得，那么我们称这个n次单位根为n次本原单位根，我们记为$$\omega_n$$。

如所有的4次单位根为：$$1$$、$$i$$、$$-1$$、$$-i$$。其中$$i$$和$$-i$$通过求0到3次幂可以导出所有的4次单位根，而$$1$$和$$-1$$却不行，所以4次本原单位根为$$i$$和$$-i$$。

对于某个n来说，它的n次本原单位根可能不止一个，但是对于所有的n来说，有两个本原单位根都是一样的，那就是$$\varepsilon_1$$和$$\varepsilon_{-1}$$，即$$e^{i\frac{2\pi}{n}}$$和$$e^{-i\frac{2\pi}{n}}$$。
原因也很简单，由n次单位根性质2，当然可以导出所有的n次单位根，只不过$$\varepsilon_1$$和$$\varepsilon_{-1}$$只是顺着数和倒着数的区别。

在这里我们需要用到一个n次本原根的性质：

$$\omega_N^{dk}=\omega_{\frac{N}{d}}^k$$

推导：$\begin{equation*}\begin{split}\omega_N^{dk}&=e^{i\frac{2\pi}{N}dk}\\&=e^{i\frac{2\pi}{\frac{N}{d}}k}\\&=\omega_{\frac{N}{d}}^k\end{split}\end{equation*}$

由DFT到FFT

有了上面的前置知识，我们开始对DFT进行加速。

我们令$$\omega_N=\varepsilon_{-1}=e^{-i\frac{2\pi}{N}}$$，所以DFT公式可以写作：

$$\displaystyle X(k)=\sum_{n=0}^{N-1}x(n)\omega_N^{nk}$$

将$$X(k)$$展开，得：

$$\displaystyle X(k)=x(0)\omega_N^{0k}+x(1)\omega_N^{1k}+x(2)\omega_N^{2k}+...+x(N-1)\omega_N^{(N-1)k}$$

我们得到了一个(N-1)次多项式，当然这里的变量是$$\omega_N^k$$，所以我们令：

$$\displaystyle A(\omega_N^k)=x(0)\omega_N^{0k}+x(1)\omega_N^{1k}+x(2)\omega_N^{2k}+...+x(N-1)\omega_N^{(N-1)k}$$

然后我们可以添加常数为0的项来把这个多项式的项数变成偶数。我们再设两个式子：

$$\displaystyle P(\omega_N^k)=x(0)\omega_N^{0k}+x(2)\omega_N^{1k}+x(4)\omega_N^{2k}+...+x(N-2)\omega_N^{\frac{N-2}{2}k}$$

$$\displaystyle Q(\omega_N^k)=x(1)\omega_N^{0k}+x(3)\omega_N^{1k}+x(5)\omega_N^{2k}+...+x(N-1)\omega_N^{\frac{N-2}{2}k}$$

这时如果我们拿$$(\omega_N^k)^2$$，即$$\omega_N^{2k}$$来替换$$P(\omega_N^k)$$和$$Q(\omega_N^k)$$的变量，再给$$Q(\omega_N^k)$$乘上一个$$\omega_N^k$$，我们又能得到：

$$\displaystyle P(\omega_N^{2k})=x(0)\omega_N^{0k}+x(2)\omega_N^{2k}+x(4)\omega_N^{4k}+...+x(N-2)\omega_N^{(N-2)k}$$

$$\displaystyle \omega_N^kQ(\omega_N^{2k})=x(1)\omega_N^{1k}+x(3)\omega_N^{3k}+x(5)\omega_N^{5k}+...+x(N-1)\omega_N^{(N-1)k}$$

不难发现：

$$\displaystyle A(\omega_N^k)=P(\omega_N^{2k})+\omega_N^kQ(\omega_N^{2k})$$

又由n次本原单位根性质，得：

$$\displaystyle A(\omega_N^k)=P(\omega_{\frac{N}{2}}^k)+\omega_N^kQ(\omega_{\frac{N}{2}}^k)$$

看见这个$$\frac{N}{2}$$，相信大家应该知道我们的加速思路是往哪里走了吧。
我们考虑二分算法，所以用$$\omega_N^{k+\frac{N}{2}}$$代替$$\omega_N^k$$，得：

$$\displaystyle A(\omega_N^{k+\frac{N}{2}})=P(\omega_{\frac{N}{2}}^k)+\omega_N^{k+\frac{N}{2}}Q(\omega_{\frac{N}{2}}^k)$$

再利用n次单位根性质3，得：

$$\displaystyle A(\omega_N^{k+\frac{N}{2}})=P(\omega_{\frac{N}{2}}^k)-\omega_N^kQ(\omega_{\frac{N}{2}}^k)$$

看！$$A(\omega_N^k)$$和$$A(\omega_N^{k+\frac{N}{2}})$$之间只差了一个负号！
所以我们只需要算出$$\omega_N^k$$、$$P(\omega_{\frac{N}{2}}^k)$$和$$Q(\omega_{\frac{N}{2}}^k)$$就能得到$$A(\omega_N^k)$$和$$A(\omega_N^{k+\frac{N}{2}})$$。

又因为P多项式的系数是A多项式的偶次项，Q多项式的系数是A多项式的奇次项，所以我们可以拆分A多项式的系数，就能得到P、Q多项式。而且由于P、Q多项式的构型又与A多项式完全一致，所以我们求解P、Q多项式的时候又可以把它们当成新的A多项式来求解。
这就是一个标准的二分算法。

至此，我们成功把时间复杂度降到了$$O(n\log_2n)$$，把DFT变成了FFT。
我们可以基于这种递归思路来实现FFT了，代码如下：

void FFT(std::complex<double>* src, int n) {
    if (n == 1) // 只有一个元素，无需计算
        return;
    int mid = n >> 1;
    static std::complex<double> dst[MAXN];
    for (int i = 0; i < mid; ++i) { // 前半为偶次项，后半为奇次项
        dst[i] = src[i << 1];
        dst[i + mid] = src[i << 1 | 1];
    }
    memcpy(src, dst, sizeof(dst));
    FFT(src, n >> 1); // 计算偶次项
    FFT(src + mid, n >> 1); // 计算奇次项
    for (int i = 0; i < mid; ++i) { // 计算原式
        // n次本原单位根的i次幂
        std::complex<double> w(cos(2 * PI * i / n), -sin(2 * PI * i / n));
        dst[i] = src[i] + w * src[i + mid];
        dst[i + mid] = src[i] - w * src[i + mid];
    }
    memcpy(src, dst, sizeof(dst));
}

快速傅里叶逆变换（IFFT）

公式为：

$$\displaystyle x(n)=\mathrm{IDFT}[X(k)]=\frac{1}{N}\sum_{n=0}^{N-1}X(k)e^{i\frac{2\pi n}{N}k}$$

我们发现FFT和IFFT的不同仅为：

1. 把FFT的输出当成输入给IFFT计算
2. IFFT的输出需要乘一个系数$$\frac{1}{N}$$
3. IFFT的n次本原单位根取的是$$\varepsilon_1$$

仅此而已，所以我们只需要稍微改动一下FFT的代码，就可以让它同时满足计算FFT和IFFT的需要：

// op = 1: FFT, op = -1: IFFT
void FFT(std::complex<double>* src, int n, int op) {
    if (n == 1) // 只有一个元素，无需计算
        return;
    int mid = n >> 1;
    static std::complex<double> dst[MAXN];
    for (int i = 0; i < mid; ++i) { // 前半为偶次项，后半为奇次项
        dst[i] = src[i << 1];
        dst[i + mid] = src[i << 1 | 1];
    }
    memcpy(src, dst, sizeof(dst));
    FFT(src, n >> 1); // 计算偶次项
    FFT(src + mid, n >> 1); // 计算奇次项
    for (int i = 0; i < mid; ++i) { // 计算原式
        // n次本原单位根的i次幂，op控制共轭
        std::complex<double> w(cos(2 * PI * i / n), -op * sin(2 * PI * i / n));
        dst[i] = src[i] + w * src[i + mid];
        dst[i + mid] = src[i] - w * src[i + mid];
    }
    memcpy(src, dst, sizeof(dst));
}

当然在递归函数里我们没法把最后一步乘以$$\frac{1}{N}$$解决，所以需要你在函数外部计算。

还可以更快？

当然！因为这份代码的常数有点大，我们需要用一些技巧来让它更快。

不要用std::complex

有些情况下STL也是会被卡的，特别是在没开O2优化时。我们最好自己实现一个复数类，只需要实现乘法、加法和减法操作即可。

struct cpx {
    double r, i;
    cpx(double a = 0.0, double b = 0.0): r(a), i(b) {}
    cpx operator+(const cpx& tmp) {
        return cpx(r + tmp.r, i + tmp.i);
    }
    cpx operator-(const cpx& tmp) {
        return cpx(r - tmp.r, i - tmp.i);
    }
    cpx operator*(const cpx& tmp) {
        return cpx(r * tmp.r - i * tmp.i, r * tmp.i + i * tmp.r);
    }
};

把递归变成非递归

不像普通的二分算法，只需要移动计算区间即可，FFT的计算是专门针对奇偶项的，每次二分会对数组进行重排，很明显不能直接展开。
所以我们还是动手实际分析一下。

当n=2时，是我们递归操作的最小单元，此时可知：

$$X[0]=x[0]+\omega_2^0\cdot x[1]$$

$$X[1]=x[0]-\omega_2^0\cdot x[1]$$

把上式转换成一张图，如下：
我们能很直观地看出X[0]和X[1]都由x[0]和x[1]转化而来，且哪个乘上了$$\omega_2^0$$也知道了。我们称这种结构为蝶形结，这也是为什么这种变换叫蝴蝶变换——它看起来就像一只蝴蝶。

如果你可以理解上面这张图，那么我们现在更上一层，看看n=4的情况：
左半部分还是我们熟悉的n=2的情况，但是似乎下标顺序好像有些不同？
我们再看右半部分：中间部分就是下标0~3分成了偶数和奇数部分！
没错，这个中间层就是第一次奇偶分项得到的下标结果。同样，右半部分也是两项得一项的原式计算部分。

我们再来看看n=8的情况：
有点眼花？我们还是一层一层来，看右边的中间层，很明显下标也是奇偶分项的；再看左边的中间层：它被分成了两个小单元，每个单元依旧是上一层按项序号顺序（不是下标）奇偶分项。

把整张图从右往左看，就是整个递归分项的过程；从左往右看，就是回溯计算的过程。
现在我们得到了最左边（最底层）的下标顺序，那么我们就可以用循环代替递归来逐层向上计算，得到最终的X[k]了！
但是好像这个下标顺序好像也看不出什么端倪来……我们把它们转化成二进制：

$$0:(0)_{10}\to(000)_2$$

$$1:(4)_{10}\to(100)_2$$

$$2:(2)_{10}\to(010)_2$$

$$3:(6)_{10}\to(110)_2$$

$$4:(1)_{10}\to(001)_2$$

$$5:(5)_{10}\to(101)_2$$

$$6:(3)_{10}\to(011)_2$$

$$7:(7)_{10}\to(111)_2$$

看起来还是没啥端倪，那我们把序号也一起二进制化：

$$(000)_2:(0)_{10}\to(000)_2$$

$$(001)_2:(4)_{10}\to(100)_2$$

$$(010)_2:(2)_{10}\to(010)_2$$

$$(011)_2:(6)_{10}\to(110)_2$$

$$(100)_2:(1)_{10}\to(001)_2$$

$$(101)_2:(5)_{10}\to(101)_2$$

$$(110)_2:(3)_{10}\to(011)_2$$

$$(111)_2:(7)_{10}\to(111)_2$$

好像下标就是序号的二进制倒位！
于是我们能得到新的算法代码：

// 要求n为2的倍数，len为n的位数-1
int n, len;
int rvs[n];
// 递推计算二进制倒位
for (int i = 1; i < n; ++i)
    rvs[i] = (rvs[i >> 1] >> 1) | ((i & 1) << (len - 1));
// FFT(op == true) or IFFT(op == false)
void FFT(cpx* A, int n, int len, bool op) {
    // Reverse
    for (int i = 0; i < n; ++i)
        if (i < rvs[i])
            std::swap(A[i], A[rvs[i]]);
    // Start calculating
    for (int i = 1; i < n; i <<= 1) {
        // FFT need sin(PI / i), but IFFT need -sin(PI / i)
        cpx wn(cos(PI / i), sin(PI / i) * (op ? 1.0 : -1.0));
        for (int j = 0; j < n; j += (i << 1)) {
            cpx w(1.0, 0.0);
            for (int k = 0; k < i; ++k, w = w * wn) {
                cpx x = A[j + k], y = w * A[i + j + k];
                A[j + k] = x + y;
                A[i + j + k] = x - y;
            }
        }
    }
    // If is IFFT
    if (!op)
        for (int i = 0; i < n; ++i) {
            A[i].r /= n;
            A[i].i /= n;
        }
}

至此，您已经学会了快速傅里叶算法（存疑）！