十分简明易懂的FFT（快速傅里叶变换）
FFT前言
快速傅里叶变换 (fast Fourier transform) ，即一种基于highly efficient 计算discrete Fourier transform (DFT) 的统一方法 ，简称 FFT 。该算法由 J.W. 库利 和 T.W. 图基 于 1965 年提出 。通过使用 FFT 算法 ，计算机实现 discrete Fourier transform 所需 的 multiplication operations 显著减少 ，尤其 当 被转换 的 sampling points N 增加时 ，FFT 算法 的 计算量节省 将更加明显 。
FFT（Fast Fourier Transformation） 是一种基于 discrete Fourier transform (DFT) 的 高效算法 。即为一种高效的 Fourier transform method 。
它是根据 discrete Fourier transform (DFT) 的 特性 对 discrete Fourier transform algorithm 进行改进而获得的一种方法 。
——百度百科
Fast Fourier Transform (缩略体形式)，中文译名：快速傅里叶变换 (用于加速多项式乘法运算)。
传统的高精度乘法算法所需时间为 O(n^2)；然而通过使用 Fast Fourier Transform 算法可以在 O(n \log_2 n) 的时间复杂度内完成这一过程。
虽然简洁高效但晦涩难懂。
通常较为晦涩难懂。
建议为学习复数与三角函数的基础知识做初步了解（即使缺乏基础也没关系）。
多项式的系数表示法和点值表示法
FFT本质上是一个基于系数形式，并且在所需的时间为O(n\log_2n)内将一个多项式转换为其点值形式的算法
多项式的系数表示和点值表示可以互相转换
系数表示法
一个最高次数为n−1且包含n个项的多项式f(x)可以用求和公式表示为f(x)=\sum^{n−1}_{i=0}a_ix^i。同样地，在这种情况下也可以通过列出各项对应的系数来表示f(x)，即f(x)=\{a_0,a_1,a_2,\dots,a_{n−1}\}。这就是所谓的系数表达方式，在平时的学习中广泛采用这种方法。
点值表示法
把多项式放到平面直角坐标系里面，看成一个函数
把n个不同的x代入，会得出n个不同的y，在坐标系内就是n个不同的点
那么这n个点精确地确定了这个多项式，并且仅此一个满足条件的多项式
其原因较为简单明了：通过将n个方程式联立求解，我们能够形成一个包含n个方程的n元方程组；在此过程中，每个方程中的各个系数均能够被确定下来。
那么f(x)还可以用f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))\}来表示
这就是点值表示法
高精度乘法下两种多项式表示法的区别
对于使用系数表示 的两个多项式而言，在处理它们时需要将它们相乘。假设这两个多项式分别定义为A(x)和B(x)。我们需要逐一取A中的每一个系数与B中的每一个系数进行计算。随后通过这些计算来确定其时间复杂度为O(n²)
但两个用点值表示 的多项式相乘，只需要O(n)的时间
什么意思呢？
分别定义为f(x)=\{(x_i, f_i)|i=0, 1, 2, ..., n-1\}和g(x)=\{(x_i, g_i)|i=0, 1, 2, ..., n-1\}其中i表示第i个离散数据点的位置。其乘积函数h=x↦h_x即为h_x = f_x \cdot g_x对所有x∈D成立其中$D={x_i|i=0, 1, 2, ..., n-1})
所以这里的时间复杂度只有一个枚举的O(n)
突然感觉高精度乘法能O(n)暴艹 一堆题？
但是朴素 的系数表示法转点值表示法的算法还是O(n^2)的，逆操作类似
基本
难道高精度乘法只能O(n^2)了吗？
DFT前置知识&技能
复数
毕竟高中有所以不多说
将形如a+bi（其中a,b皆为实数值）的量定义为复数量 ，其中a被称作其实部部分而b则被视为虚部部分。若一复数量具有零值的虚部，则可视为实数值量；若其实部部分为零但虚部不为零，则该复数量被称为纯虚数值量。
复体域可被视为实数值域上的代数闭包体系，在此体系中任何具复数值系数的一元多项式方程必存在至少一个根。
有关这一概念的历史发展，则可追溯至意大利米兰学者吉罗拉莫·卡尔达诺于16世纪初首次提出相关理论。
此后随着让·达朗贝尔、亚伯拉罕·棣莫弗、莱昂哈德·欧拉以及卡尔·弗里德里希·高斯等人深入研究与推广，在随后的发展过程中这一概念逐渐被数学界所认可。
初中数学老师通常会告诉我们√(-1)不存在，在实数范围内。
扩展到复数集C中，并定义i为虚数单位且满足i²=-1。
一个复数z可以用形式z=a+bi（其中a,b属于R）来表示。
在这里：
-
a被称为该复数的实部
-
b被称为该复数的虚部
-
i被称为虚数单位
例如，在这个表达式中， -
a是一个实数值，
-
b也是一个实数值，
-
i则是用来表示虚的部分并满足特定性质的一个特殊符号。
- 在复数集中就可以用i表示负数的平方根，如\sqrt{-7}=\sqrt{7}i
还可以把复数看成复 平面直角坐标系上的一个点 ，比如下面
x轴就是实数集中的坐标轴，y轴就是虚数单位i轴
这个点(2,3)对应的复数就是2+3i，也可以想象它作为向量表示为(2,3)。此外，在极坐标系中这一特定位置还可以用(\sqrt{13},\theta)来描述。对于一个任意的复数z=a+bi来说，在其极坐标形式中其实可以被重新定义为其模与辐角的形式：即r=|z|=\sqrt{a^2+b^2}；而它的辐角\theta = \arctan{\frac{b}{a}}(其中a≠0). 这样一来该复数的共轭复数就可以被重新定义为其实部不变而虚部符号相反的情况：即\overline{z}=a-bi
复数的运算
复数不同于点或向量，并与实数一样能够执行基本的算术运算。
设两个复数分别为 z_1=a+bi, z_2=c+di，
那么 z_1 + z_2 = (a + c) + (b + d)i
以及 z_1 z_2 = (ac - bd) + (ad + bc)i
复数不同于点或向量，并与实数一样能够执行基本的算术运算。
设两个复数分别为 z_{1}=a+b\text{i}, z_{2}=c+d\text{i}，
那么 z_{1} + z_{2} = (a + c) + (b + d)\text{i}
以及 z_{1} z_{2} = (ac - bd) + (ad + bc)\text{i}。
复数相加 也满足平行四边形法则
这张是从网上盗的
即AB+AD=AC
复数形式的乘法运算还具有一个显著的关键特性。
(a_{1},\theta_{1})*(a_{2},\theta_{2})=(a_{1}a_{2},\theta_{1}+\theta_{2})
即其模长的乘积与极角之和
DFT（离散傅里叶变换）
- 一定注意从这里开始所有的 n都默认为 2的整数次幂
对于任意系数多项式转点值来说，则是方便地选取任意n个x值来进行计算。然而直接计算每一个x_k的不同幂次会导致算法复杂度上升至O(n²)，这显然是一个较高的效率问题。选择一组巧妙安排好的x坐标点能够显著减少运算量。这些特定位置上的点并非随意选取而是经过精心设计以确保高效运算
如果我们选择一些变量x使得每个x满足某些幂次后等于1的话，则无需进行全部的幂运算
- 傅里叶说：这个圆圈上面的点都可以做到
With the origin as the center, draw a unit circle with radius 1. Thus, any point on the unit circle can be raised to a power to obtain 1. Fourier mentions it should also be divided into n equal parts; for example, when n=8.
橙色点即为n=8时要取的点，从(1,0)点开始，逆时针从0号开始标号，标到7号
记编号为k的点代表的复数值为\omega_n^k，那么由模长相乘，极角相加 可知(\omega_n^1)^k=\omega_n^k
其中\omega_n^1称为n次单位根 ，而且每一个\omega都可以求出 （我三角函数不好）
\omega_n^k=\cos{k\over n}2π+i\sin{k\over n} 2π
或者说也可以这样解释这条式子
注意sin^2\theta+cos^2\theta=1什么的，就容易理解了
那么\omega^0_n,\omega^1_n,...,\omega^{n-1}_n即为我们要代入的x_0,x_1,...,x_{n-1}
单位根的一些性质
推FFT的过程中需要用到\omega的一些性质
\omega^k_n=\omega^{2k}_{2n}
它们表示的点（或向量）表示的复数是相同的
证明
\omega^k_n=cos{k\over n}2π+isin{k\over n} 2π=cos{2k\over 2n}2π+isin{2k\over 2n} 2π=\omega^{2k}_{2n}
\omega^{k+{n \over 2}}_n=-\omega_n^k
它们表示的点以原点为中心对称；所表示的复数其实部互为相反数；所表示的向量大小相等且方向相反
证明
\omega^{n\over 2}_n=cos{{n\over 2}\over n}2\pi+isin{{n\over 2}\over n}2\pi=cos\pi+isin\pi=-1
（这个东西和e^{ix}=cosx+isinx与e^{i\pi}+1=0有点关系，我不会就不讲了 ）
\omega^0_n=\omega^n_n
- 都等于1，或1+0i
FFT（快速傅里叶变换）
尽管DFT创造出了大量非常厉害的ω, 用作多项式x值的代入；然而仅仅将此类特殊x值用于变换的过程被称为DFT。然而实际上还是要采用直接而繁琐的方法——即代入单位根进行计算。
尽管DFT创造出了大量非常厉害的ω, 用作多项式x值的代入；然而仅仅将此类特殊x值用于变换的过程被称为DFT。然而实际上还是要采用直接而繁琐的方法——即代入单位根进行计算。
- DFT还是暴力 O(n^2)…
然而DFT可以通过递归的方式来实现，并由此引出了快速傅里叶变换FFT（Fast Fourier Transform）。我们通常会设定一个多项式A(x)。即：首先设一个多项式A(x)
A(x) = \sum^{n-1}_{i=0} a_i x^i = a_0 + a_1 x + a_2 x^2 + \dots + a_{n-1} x^{n-1}
按A(x)下标的奇偶性 把A(x)分成两半，右边再提一个x
A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})
=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})
两边好像非常相似，于是再设两个多项式A_1(x),A_2(x)，令
A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{{n\over 2}-1}A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{{n \over 2}-1}
比较 明显得出
A(x)=A_1(x^2)+xA_2(x^2)
令k<{n\over 2}并以\omega^k_n替代x进入多项式A(x)中。注意：以下步骤中，请重点考虑单位根的相关性质。
A(\omega^k_n)=A_1((\omega^k_n)^2)+\omega^k_nA_2((\omega^k_n)^2) 等于 A_1(\omega^{2k}_n)+\omega^k_n A_2(\omega^{2k}_n) ，这又等于 A_1(\omega^k_{n over 2})+\omega^k_n A_2(\omega^k_{n over 2}) 。
当我们考虑将\omega^{k+\frac{n}{2}}_{\, n}代入到表达式中时，
我们可以得到以下结果：
A\left( \omega _{\, n}^{k+\frac{n}{2}} \right) = A_{1}\left( \left( \omega _{\, n}^{2} \right)^{kn + n} \right) + \left( \cos kn + i \sin kn \right) A_{1}\left( ...
- 发现了什么？
对于这两个多项式 A(\omega^k_n) 和 A(\omega^{k+\frac{n}{2}}_n}, 它们的其余部分具有相同的系数, 仅相差一个符号。
也就是说, 如果已知 $A_1} (\omega^\frac{k} {{\frac{n} {}}} ) 和 A_1} (\omega^\frac{k} {{\frac{n} {}}} ) 的值, 我们便能够 共同得到 A (ωkn) 和 A (ω^{k+\frac{n}{₂}}} _n) 的值。
基于此, 我们便可以采用 递归分治的方法 来实现 FFT 的计算过程.
在回溯过程中仅需扫描当前序列的一半即可推导出另一半的解
当n=1时只有一个常数项，则直接返回该常数值
算法的时间复杂度具有O(n\log_2n)的阶
IFFT（快速傅里叶逆变换）
思考一下是否正确？除了掌握Fast Fourier Transform (FFT)之外还需学会Inverse Fast Fourier Transform (IFFT)（快速傅里叶逆变换）。将两个多项式进行卷积运算并完成两次Fast Fourier Transform (FFT)计算以得到其点值形式的信息。然而，在系数形式下这些信息并不直观有用。
怎么把点值表示的多项式快速 转回系数表示法？
IDFT暴力O(n^2)做？其实也可以用FFT用O(n\log_2n)的时间搞
您是否想过傅里叶为何选择了使用ω^k_n作为x的值而非其他数值呢？这一选择的原因确实存在但目前仍难以理解其背后的逻辑。然而个人认为只需要掌握一个核心结论即可。
一个多项式在分治的过程中乘以单位根的共轭复数，在分而治之的过程中被除以 n 即相当于原始多项式的每一项系数
意思就是说FFT和IFFT可以一起搞
朴素版FFT板子
c++内置了一个名为模板complex的库
a.real()即返回复数a$的实部值
#include<complex>
#define cp complex<double>
void fft(cp *a,int n,int inv)//inv是取共轭复数的符号
{
if (n==1)return;
int mid=n/2;
static cp b[MAXN];
fo(i,0,mid-1)b[i]=a[i*2],b[i+mid]=a[i*2+1];
fo(i,0,n-1)a[i]=b[i];
fft(a,mid,inv),fft(a+mid,mid,inv);//分治
fo(i,0,mid-1)
{
cp x(cos(2*pi*i/n),inv*sin(2*pi*i/n));//inv取决是否取共轭复数
b[i]=a[i]+x*a[i+mid],b[i+mid]=a[i]-x*a[i+mid];
}
fo(i,0,n-1)a[i]=b[i];
}两个多项式a,b相乘再转系数多项式c，通常只用打这么一小段
cp a[MAXN],b[MAXN];
int c[MAXN];
fft(a,n,1),fft(b,n,1);//1系数转点值
fo(i,0,n-1)a[i]*=b[i];
fft(a,n,-1);//-1点值转系数
fo(i,0,n-1)c[i]=(int)(a[i].real()/n+0.5);//注意精度显然, FFT 专用于计算 n 为 2 的整数次幂时的多项式形式. 因此, 在 FFT 运算前, 必须将 n 设置为符合这一条件.
这个板子看着好像 很优美，但是
递归常数太大 ，要考虑优化…
FFTの优化——迭代版FFT
这个图也是盗的
这个很容易发现点什么吧？
- 每个位置分治后的最终位置为其二进制翻转后得到的位置
这样处理的话, 我们可以先对原序列进行变换, 将每个元素放置到其最终位置上; 然后逐步向上合并. 只需一句代码即可完成对第i位置上的元素进行预处理, 并计算其逆序数rev[i].
fo(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));至于蝴蝶变换它死了 其实是我不会
真·FFT板子
void fft(cp *a,int n,int inv)
{
int bit=0;
while ((1<<bit)<n)bit++;
fo(i,0,n-1)
{
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次（就是没交换）
}
for (int mid=1;mid<n;mid*=2)//mid是准备合并序列的长度的二分之一
{
cp temp(cos(pi/mid),inv*sin(pi/mid));//单位根，pi的系数2已经约掉了
for (int i=0;i<n;i+=mid*2)//mid*2是准备合并序列的长度，i是合并到了哪一位
{
cp omega(1,0);
for (int j=0;j<mid;j++,omega*=temp)//只扫左半部分，得到右半部分的答案
{
cp x=a[i+j],y=omega*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=x-y;//这个就是蝴蝶变换什么的
}
}
}
}这个板子好像不是那么好背
至少这个板子已经很优美 了
FFT后记
本人版权意识薄弱……
部分知识源自上述提到的资源中的一部分内容。这些资源包括：
- <https://www.cnblogs.com/RabbitHu/p/FFT.html]：该文详细介绍了快速傅里叶变换（FFT）的基本原理及其应用；
- <https://www.cnblogs.com/zwfymqz/p/8244902.html?mType=Group#_label3]：该文探讨了深度学习模型在图像识别任务中的性能优化方法；
- **<https://blog.cdn.net/ggn_2015/article/details/68922404/**：该文章深入分析了数据预处理技术在机器学习项目中的重要性。
此外，请参考上述资源以获得更详细的信息。
NTT我来了