快速傅里叶变换（FFT）学习笔记

KobeBeanBryantCox · 2025-03-27 11:38:30 · 算法·理论

快速傅里叶变换（FFT）学习笔记

前言

这是我第二篇学习笔记，上一篇 SAM 的在这里。

我一般不写学习笔记，除非太难理解或者算法本身太过巧妙。

不同于 SAM，这次写学习笔记的原因是因为其太过巧妙了。

前置知识：矩阵，弧度制，三角函数，代数基础，函数基础，复数概念。

参考资料：这里。本学习笔记详细了视频的内容，然后那个视频还是有错误的，我将在这里指正其错误。

如果任何一句话您看不懂，请先尝试找 DeepSeek，它很聪明。如果再看不懂，请评论或者私信我。

01. 引入

如果我们有两个多项式 A(x)=x^2+3x+2 和 B(x)=2x^2+1。

我们要把他相乘得到 C(x)=A(x)B(x)，我们平时是怎么做的？

在数学上，我们使用乘法分配律，如下：

这是数学上的做法。

在计算机中，首先一个重要的事情是，我们如何存储这个多项式？

最自然的做法是用数组存储多项式的系数，如下：

• \sin\theta=\dfrac{\text{对边}}{\text{斜边}}=\dfrac b r\Rightarrow b=r\sin\theta
• 所以 z=a+bi=r\cos\theta+i\cdot r\sin\theta=r(\cos\theta+i\sin\theta)。

因为有欧拉公式：e^{i\theta}=\cos\theta+i\sin\theta，我们可以得到 z=re^{i\theta}。

所以 z^n=r^ne^{in\theta}。

现在来解 z^n=1。

因为 1=1\cdot e^{i\cdot0}，所以 z^n=1\cdot e^{i\cdot 0}=r^ne^{in\theta}。

为了使其相等，我们要满足两个条件：

1. 模相等：r^n=1；
2. 辐角相差 2\pi 的倍数（弧度制下 2\pi 是 360 度，转了一圈）：n\theta=0+2\pi k,k\in\mathbb{Z}。

解得：

• 模 r=1；
• 辐角 \theta=\dfrac{2\pi k}n。

由于 r=1，所以它们都在单位圆上。

根据解得的 \theta，相邻两个辐角差为 \dfrac{2\pi(k+1)}n-\dfrac{2\pi k}n=\dfrac{2\pi}n。

直观理解就是把 2\pi 的圆周均匀分成了 n 份，所以其等距分布。

我们这就证明了刚开始得出的结论。

此时 z=re^{i\theta}=1\cdot e^{i\cdot\frac{2\pi k}n}，这就是 z^n=1 的解了，其中 k\in[0,n-1]。

我们定义 \omega_n=e^{\frac{2\pi i}n}，我们称这个 \omega_n 叫做单位根。

不难发现 [\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}] 就是 1 的 n 个 n 次方根，因为代入它们后，会发现其与上文的 z=1\cdot e^{i\cdot\frac{2\pi k}n},k\in[0,n-1] 是相等的。

所以我们的 FFT 的第一层递归此时就是要在 1,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1} 处取点。

回到原来的问题，为什么这些 \omega_n^k 就能用来递归呢，具体的，为什么其一定有成对的相反数呢？

上图是 n=16 时的单位圆和单位根（上文强调过，n 只能取 2 的整数次幂）。

我们可以很清楚地看见 \omega_{16}^{j+\frac n 2}=-\omega_{16}^j（注意向量的复平面与平面直角坐标系的区别，复平面的负号是要求与原点对称），这个叫做单位根的对称性。

（其实根据这个还能推出单位根的周期性，不过不重要。）

上文说过，第 i 层递归用的数应该是第 i-1 层用的数的平方。

所以第二层用的数就是 [(\omega_{16}^0)^2,(\omega_{16}^1)^2,\cdots,(\omega_{16}^{7})^2]（只取到 7 是因为 8\sim15 的点与 0\sim7 的点分别关于原点对称，平方后是同一个值。）

你会发现这些东西居然是 [\omega_{16}^0,\omega_{16}^2,\omega_{16}^4,\cdots,\omega_{16}^{12},\omega_{16}^{14}]。

然后它居然等价于 [\omega_8^0,\omega_8^1,\cdots\omega_8^7]（可以画画图理解，会发现两者等价）。

（\omega_{2n}^{2k}=\omega_n^k 其实叫做折半引理，不过这也不重要。）

然后这些居然也是具有对称性的。再平方后也是同理的，直到只有一个值。

这样就选完数了。FFT 也就是这样了（是不是很简单），下面来对着代码总结 FFT 的过程。

FFT 递归版代码

别看我们说了这么多，其实代码超级短（12 行）！

我们一点一点来。

typedef double lf;
const lf PI=acos(-1);
void FFT(complex<lf>*a,int n)
{
}

这是定义 FFT 函数，表示将 a 这个系数表达式转化成值表达式，返回在 a 中。

complex 是 c++ STL 提供的复数模板。

传入的时候所有 a_k 的虚部都是 0i。

这里 n 是长度，传入的是大于应该有的次数的最小 2 的整数幂。

if(n==1)return;

长度为 1 直接退出。

int m=n>>1;
complex<lf>a0[m],a1[m];
for(int i=0;i<n;i+=2)a0[i>>1]=a[i],a1[i>>1]=a[i+1];

把其分成更小的长度 m，然后把 a 分为奇偶两部分。

FFT(a0,m),FFT(a1,m);

递归处理奇偶两部分。

complex<lf>wn={cos(2*PI/n),sin(2*PI/n)},w={1,0};

定义单位根 w_n。注意这里的 w 是单位根的幂次，会在循环中不断乘 w_n。

事实上，你使用：

complex<lf>wn=exp(complex<lf>{0,2*PI/n});

来定义单位根也是可以的（根据欧拉公式）。

但是很多题解用的都是第一种。我们就跟随潮流吧。（或许是因为精度的误差？但两种都能过模板题。有知道的麻烦告知我。）

for(int i=0;i<m;i++,w*=wn)a[i]=a0[i]+w*a1[i],a[i+m]=a0[i]-w*a1[i];

然后你把它们合起来就好了（根据先前计算的公式）。

为什么合起来的时候一定是前一半后一半呢？

先前的公式（去前文看）：

\fbox{$\begin{array}{l}f(x_j)&=f_e(x_j^2)+x_j\times f_o(x_j^2)\\f(-x_j)&=f_e(x_j^2)-x_j\times f_o(x_j^2)\end{array}$}

我们知道此时 x_j=\omega^j，所以如下：

\fbox{$\begin{array}{l}f(\omega^j)&=f_e(\omega^{2j})+\omega^j\times f_o(\omega^{2j})\\f(-\omega^j)&=f_e(\omega^{2j})-\omega^j\times f_o(\omega^{2j})\end{array}$}

我们前文也说过 -\omega^j=\omega^{j+\frac n 2}=\omega^{j+m}（对称性），所以：

\fbox{$\begin{array}{l}f(\omega^j)&=f_e(\omega^{2j})+\omega^j\times f_o(\omega^{2j})\\f(\omega^{j+m})&=f_e(\omega^{2j})-\omega^j\times f_o(\omega^{2j})\end{array}$}

我们 f 的下标用的就是 \omega 的上标，所以一定是前一半后一半地存。

总代码如下：

typedef double lf;
const lf PI=acos(-1);
void FFT(complex<lf>*a,int n)
{
    if(n==1)return;
    int m=n>>1;
    complex<lf>a0[m],a1[m]; // 这样开数组其实有点不太规范的
    for(int i=0;i<n;i+=2)a0[i>>1]=a[i],a1[i>>1]=a[i+1];
    FFT(a0,m),FFT(a1,m);
    complex<lf>wn={cos(2*PI/n),sin(2*PI/n)},w={1,0};
    for(int i=0;i<m;i++,w*=wn)a[i]=a0[i]+w*a1[i],a[i+m]=a0[i]-w*a1[i];
}

FFT 非递归代码

为了阅读的连续性，我不打算先讲 IFFT，心急的读者可以跳到后文中“插值”那一部分。但不建议这么做。

由于递归时动态开 complex 数组巨大的常数，有些题是会 TLE 的（虽然模板题过了），我们需要学习非递归版的。

我们来看一下我们的递归过程：

step 1: | 0   1   2   3   4   5   6   7 |
step 2: | 0   2   4   6 | 1   3   5   7 |
step 3: | 0   4 | 2   6 | 1   5 | 3   7 |
step 4: | 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 |

其中 step 是我们递归的层数。

我们观察一下第一层和第四层：

0 1 2 3 4 5 6 7
0 4 2 6 1 5 3 7

我们写出其二进制：

000 001 010 011 100 101 110 111
000 100 010 110 001 101 011 111

发现什么没有？我们发现第一层的二进制反过来就变成了第四层的二进制。

所以我们可以根据第一层的二进制，来预处理出来第四层的二进制。

void init(int len)
{
    int b=log2(len);
    for(int i=1;i<len;i++)rev[i]=(rev[i>>1]>>1)+((i&1)<<(b-1));
}
/*
这个 rev 就是 i 的二进制倒过来后的二进制的数
为什么可以这样递推呢？

我们考虑 2*i=(i<<1)
那么反过来的二进制 rev[2*i]=(rev[i]>>1)
所以可以写成 rev[i]=(rev[i>>1]>>1)，这里 i 为偶数

考虑到 2*i+1=(i<<1)+1
反过来的二进制是 rev[2*i+1]=(rev[i]>>1)+1
所以可以写成 rev[i]=(rev[i>>1]>>1)+(1<<(b-1))，这里 i 为奇数，其中 b 是 log2(len)

把它们合起来，就变成了最终形式：rev[i]=(rev[i>>1]>>1)+((i&1)<<(b-1))
*/

然后根据递归的合并过程，我们用迭代模拟递归，即枚举长度 1,2,4，然后两两合并。具体的看代码。

void FFT(complex<lf>*a,int n)
{
    for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]); // 按照刚刚的把 a 的顺序调整到最后一层
    for(int len=1;len<=(n>>1);len<<=1) // 枚举合并后的长度的一半
    {
        complex<lf>wn={cos(PI/len),sin(PI/len)}; // 单位根。注意此时由于枚举的长度是一半，所以 PI 不需要 *2
        for(int i=0;i<n;i+=(len<<1)) // 枚举合并的每一个区间的起始位置
        {
            complex<lf>w={1,0};
            for(int j=0;j<len;j++,w*=wn) // 枚举区间的每一个元素
            {
                complex<lf>a0=a[i+j],a1=w*a[i+j+len]; // 好像这里叫做蝴蝶操作，不过这不重要，没什么用
                a[i+j]=a0+a1,a[i+j+len]=a0-a1; // 正常计算即可
            }
        }
    }
}

这样非递归版的居然比递归版的快了近 1 倍。

插值

其实完整的 FFT 并没有做完。

我们刚刚只是把系数表达式转化成了值表达式。

现在把其对应位置相乘得到了答案。

但是答案依然是值表达式，我们还需要转化成系数表达式，我们称这一步为插值。

一般把这一步叫做 IFFT（Inverse Fast Fourier Transform，快速傅里叶逆变换）。

做法的解释有很多，我用矩阵来解释一下。

用矩阵的角度看待 FFT

想一想求值的时候我们是在干什么，是不是找一堆 x 然后代入求 f(x)？

就是这些式子：

实际上可以写成矩阵的形式：

我们知道此时 x_k=\omega^k，可以改写成如下形式：

事实上那个 n\times n 的矩阵被称为离散傅里叶变换（DFT）矩阵，不过这不重要。

用矩阵的角度看待 IFFT

事实上本质上我们要根据很多 (x,f(x)) 的二元组解出所有 p_k。

这不就直接把 DFT 求个逆就好了吗：

这个逆矩阵我让 DeepSeek 帮我们求了，所以改写成如下：

诶这是不是跟之前的 DFT 超级像？

IFFT 代码

所以可以直接把之前的 FFT 拿过来，然后改一行代码，再加一行代码即可：

void IFFT(complex<lf>*a,int n)
{
    for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int len=1;len<=(n>>1);len<<=1)
    {
        complex<lf>wn={cos(PI/len),-sin(PI/len)}; // 这里改成虚部 - 的
        for(int i=0;i<n;i+=(len<<1))
        {
            complex<lf>w={1,0};
            for(int j=0;j<len;j++,w*=wn)
            {
                complex<lf>a0=a[i+j],a1=w*a[i+j+len];
                a[i+j]=a0+a1,a[i+j+len]=a0-a1;
            }
        }
    }
    for(int i=0;i<n;i++)a[i]={a[i].real()/n+0.5,a[i].imag()}; // 再加上这一行，矩阵最前方的 1/n
}

定义单位根时加负号是因为 \omega^{-1}=e^{\frac{-2\pi i}{n}}=\cos(\frac{2\pi}n)-i\sin(\frac{2\pi}n)。

有没有注意到 +0.5 这个操作？那是因为复数的精度误差太大了！

这代码这么像，两个东西显然可以合并啊：

int rev[N];
void init(int len)
{
    int b=log2(len);
    for(int i=1;i<len;i++)rev[i]=(rev[i>>1]>>1)+((i&1)<<(b-1));
}
void FFT(complex<lf>*a,int n,int type=1)
{
    for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int len=1;len<=(n>>1);len<<=1)
    {
        complex<lf>wn={cos(PI/len),type*sin(PI/len)};
        for(int i=0;i<n;i+=(len<<1))
        {
            complex<lf>w={1,0};
            for(int j=0;j<len;j++,w*=wn)
            {
                complex<lf>a0=a[i+j],a1=w*a[i+j+len];
                a[i+j]=a0+a1,a[i+j+len]=a0-a1;
            }
        }
    }
}
void IFFT(complex<lf>*a,int n){FFT(a,n,-1);for(int i=0;i<n;i++)a[i]={a[i].real()/n+0.5,a[i].imag()};}

这就是全部的 FFT 了。

03. 总结

• 我们用多项式乘法提出了 FFT。
• 我们思考把多项式用对应点的值表示。
• 多项式的值需要选取合适的点。
• 我们首先想到了利用相反数。
• 但是这样不方便递归，所以我们想到了复数。
• 我们使用 1 的 n 个 n 次方根，这样平方完，依然是正负成对出现的。
• 这也是 FFT 的精髓。
• 当我们反过来做插值的时候，我们发现 IFFT 居然可以用 FFT 来解决。
• 这样就做完了所有的问题。
• 不得不感叹 FFT 真的太精巧了。

04. 模板题代码

模板题在这里。

#include<bits/stdc++.h>
#define Code using
#define by namespace
#define wjb std
Code by wjb;
int in()
{
    int k=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9')k=k*10+c-'0',c=getchar();
    return k*f;
}
void out(int x)
{
    if(x<0)putchar('-'),x=-x;
    if(x<10)putchar(x+'0');
    else out(x/10),putchar(x%10+'0');
}
typedef double lf;
const lf PI=acos(-1);
const int N=(1<<21)+20;
int rev[N];
void init(int len)
{
    int b=log2(len);
    for(int i=1;i<len;i++)rev[i]=(rev[i>>1]>>1)+((i&1)<<(b-1));
}
void FFT(complex<lf>*a,int n,int type=1)
{
    for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int len=1;len<=(n>>1);len<<=1)
    {
        complex<lf>wn={cos(PI/len),type*sin(PI/len)};
        for(int i=0;i<n;i+=(len<<1))
        {
            complex<lf>w={1,0};
            for(int j=0;j<len;j++,w*=wn)
            {
                complex<lf>a0=a[i+j],a1=w*a[i+j+len];
                a[i+j]=a0+a1,a[i+j+len]=a0-a1;
            }
        }
    }
}
void IFFT(complex<lf>*a,int n){FFT(a,n,-1);for(int i=0;i<n;i++)a[i]={a[i].real()/n+0.5,a[i].imag()};}
complex<lf>a[N],b[N];
int main()
{
    int n=in(),m=in();
    for(int i=0;i<=n;i++)a[i]={in(),0};
    for(int i=0;i<=m;i++)b[i]={in(),0};
    int len=1;while(len<=n+m)len<<=1;
    init(len);
    FFT(a,len),FFT(b,len);
    for(int i=0;i<len;i++)a[i]*=b[i];
    IFFT(a,len);
    for(int i=0;i<=n+m;i++)out(a[i].real()),putchar(' ');
    return 0;
}

尾言

\sout{\LaTeX} 太难打啦啊啊！

在代码的海洋里，我们都是追逐算法星光的旅人。快速傅里叶变换在时光里生长，每一轮迭代都凝结着前人智慧的结晶。或许此刻你正为迭代实现的位逆序而困惑，为单位根的选择而辗转，但若将视角拉长，这些多项式终将编织成照亮前路的银河。

对所有 OIer 而言，那些调试到凌晨的夜晚，那些反复优化的常数，那些 AC 后颤抖的指尖，都将沉淀为生命里珍贵的底色。无论未来是站在领奖台的聚光灯下，还是转身走向更广阔的天地，请记得：你曾在 OI 的星空中留下过属于自己的轨迹。

愿我们永远保持对未知的好奇，像 FFT 处理多项式般优雅地接纳所有可能的挑战。无论未来走向何方，这段为 OI 燃烧的岁月，都将是我们回望时最璀璨的星辰。前路漫漫亦灿灿，愿所有努力终不被辜负，AC 常伴，未来可期。