快速傅里叶变换FFT(Fast Fourier Transform)

快速傅里叶变换(Fast Fourier Transform,FFT)用来计算离散傅里叶变换（Discrete Fourier Transform,DFT）及其逆变换（IDFT），本质上多项式乘法实际上是多项式系数向量的卷积（convolution）
（以上百度）

看了很多网上关于FFT的讲解，有一些是直接忽略了公式的推导，另外一些有推导，但是推导中的细节却没有讲清楚。本着不懂就学的心态，我把FFT的思维和推导细节用公式讲清楚，方便后人能更细致地学习FFT。

在了解FFT之前，需要有一些前置的知识，以下为目录。

欧拉函数
复数及单位根
多项式的系数、点值表达和多项式相乘
FFT的思维过程（Cooley-Tukey 算法）
4.1 DFT
4.2 IDFT
FFT的实现过程以及FNT

一、欧拉函数

$\tag{1} e^{i\theta}=cos\theta+isin\theta$ 其中i为虚数单位即 $\tag{2} i^2 = -1$ 即虚数单位 $\sqrt{-1}$

二、复数及单位根

复数形式： $a+bi$ ，其中i为虚数单位

复数乘法：对于两个复数 $A=a_0+b_0i$ 和 $B=a_1+b_1i$ ，则 $A*B = a_0a_1+b_0b_1i^2+(a_0b_1+a_1b_0)i=a_0a_1-b_0b_1+(a_0b_1+a_1b_0)i$
由于欧拉公式（见公式1）令 $a=rcos\theta , b=rsin\theta$ 则复数 $\tag{3} a+bi = r(cos\theta + isin\theta)=re^{i\theta}$
其中 $r$ 为该复数所在复平面圆的半径， $\theta$ 为该复数在复平面中的幅角。则两个复数为 $A=r_Ae^{i\theta_{A}}, B=r_Be^{i\theta_{B}}$ ，即 $\tag{4}A*B=r_Ar_Be^{i(\theta_{A}+\theta_{B})}$
根据4式可得，两个复数的相乘可以看作是幅角相加，模长相乘。

单位根：对于满足 $w^n=1$ 方程的复数，我们称其为n次单位根。由于根据复数乘法，我们可知：复数相乘为幅角相加，模长相乘。则对于每个单位根，模长为1，幅角的n倍为0。即 $r=1,sin{(n \theta )} = 0$ （易得）。

为了保证幅角的n倍始终为0，由于 $sin(2k\pi)=0,k=0,1,2,...,n$ 这个性质，我们可以将单位根表示为 $e^{\frac{2k\pi i}{n}},k=0,1,2,...,n$ ，其中 $r=1,\theta=\frac{2k\pi }{n},k=0,1,2...,n$ 。

我们发现无论k取值， $\theta$ 的n倍始终为0。

记 $w_n=e^{\frac{2k\pi}{n}}$ ，则n次单位根可以表示为 $w_n^{0},w_n^{1},w_n^{2},...,w_n^{n-1}$

三、多项式的系数、点值表达和多项式相乘

多项式的系数表达：给定一个多项式 $\tag{5} f(x)=\sum_{i=0}^{n}a_ix^i =a_0x^0+a_1x^1+a_2x^2+...+a_{n-1}x^{n-1}+a_nx^n$
其中 $\mathbf {a} = (a_0,a_1,a_2,...,a_n)$ 为系数向量

多项式的点值表达：给定一个多项式如公式（5），我们将其离散化，设取 $n+1$ （这里为什么是n+1项，将在第四节中讲到）互不相同的值 $P = \lbrace p_0,p_1,p_2,p_3,...,p_{n-1},p_n\rbrace$ ，将其代入可得 $\tag{6} A(x)=\lbrace(p_0,f(p_0)),(p_1,f(p_1)),(p_2,f(p_2)),...,(p_n,f(p_n))\rbrace$ ，则 $A(x)$ 为 $f(x)$ 在 $P$ 上的点值表达。

多项式系数表达的乘法：给定两个多项式 $\tag{6} A(x)=\sum_{i=0}^na_ix^i=a_0x^0+a_1x^1+a_2x^2+...+a_{n-1}x^{n-1}+a_nx^n$
$\tag{7} B(x)=\sum_{i=0}^nb_ix^i=b_0x^0+b_1x^1+b_2x^2+...+b_{n-1}x^{n-1}+b_nx^n$
则多项式系数表达的乘法为 $\tag{8} C(x)=A(x)*B(x)=\sum_{i=0}^{2n}c_ix^i$
其中:
$\tag{9} c_i = \sum_{j+k=i,0\leq j,k\leq n} a_jb_kx^i$
有式（9）可得，计算复杂度为 $\omicron(n^2)$

多项式点值表达的乘法：给定两个多项式如式（6）与（7），则其在 $P=\lbrace p_0,p_1,p_2,...,p_{n-1},p_n \rbrace$ 上的点值表达分别为：
$\tag{10} A'(x)=\lbrace(p_0,A(p_0)),(p_1,A(p_1)),(p_2,A(p_2)),...,(p_n,A(p_n))\rbrace$
$\tag{11} B'(x)=\lbrace(p_0,B(p_0)),(p_1,B(p_1)),(p_2,B(p_2)),...,(p_n,B(p_n))\rbrace$
则多项式点值表达的乘法为
$\tag{12} C'(x) =A'(x)*B'(x)=\lbrace (p_0,A(p_0)*B(p_0)),(p_1,A(p_1)*B(p_1)),(p_2,A(p_2)*B(p_2)),...,(p_n,A(p_n)*B(p_n))\rbrace$
可见，当我们已知 $A'(x)和B'(x)$ 即可在 $\omicron(n)$ 的复杂度下求得结果多项式的点值表达。

四、FFT的思维过程（Cooley-Tukey 算法）

对于一个多项式的乘法，根据上述前置知识的补充，我们可以得知：降低多项式乘法复杂度的方法是将常见的多项式系数表达转变为多项式的点值表达，做完点值表达的乘法后，最后再将点值表达转化为系数表达，即可完成多项式乘法。
所以问题转变为：
1.如何将多项式系数表达转变为多项式点值表达
2.如何将多项式点值表达转变为多项式系数表达
由此引出了离散傅里叶变换 DFT（Discrete Fourier Transformation）和逆离散傅里叶变换 IDFT(Inverse Discrete Fourier Transformation)

4.1 离散傅里叶变换DFT(Discrete Fourier Transformation)

离散化多项式的一种方法是将值代入到多项式中，依次求出点值。显而易见，这种方法的复杂度是 $\omicron(n^2)$ 的，这与我们降低复杂度的想法是冲突的。
于是我们引入了FFT的经典算法——Cooley-Tukey 算法，来降低离散化的复杂度，这是一个基于分治策略的算法。

给定一个多项式 $\tag{13} A(x)=\sum_{i=0}^na_ix^i=a_0x^0+a_1x^1+a_2x^2+...+a_{n-1}x^{n-1}+a_nx^n$
我们将其根据奇偶项分成两个项数相同的多项式(将多项式补充到 $2^k项,n\leq 2^k$ ，补充项数系数为0。PS：为什么是 $2^k$ 项呢，后续将会提及)：
$\tag{14} A_0(x)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}x^{2i}=a_0x^0+a_2x^2+...+a_nx^n$
$\tag{15} A_1(x)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}x^{2i+1}=x\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}x^{2i}=a_1x^1+a_3x^3+...+a_{n-1}x^{n-1}$
显而易见： $A(x)=A_0(x)+A_1(x)$

在进一步之前：我们需要返回单位根的知识点。根据n次单位根的表达 $w_n^k,k=0,1,2,...,n-1$ ，我们可以获得一个等式
$\tag{16} w_n^{\frac{n}{2}+k}=w_n^{\frac{n}{2}}·w_n^k=-w_n^k$

我们将其代入式子（14），（15）可得：
$\tag{17} A_0(w_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}w_n^{2ki}=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}w_{\frac{n}{2}}^{ki}$
$\tag{18} A_1(w_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}w_n^{2ki+1}=w_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}w_n^{2ki}=w_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}w_{\frac{n}{2}}^{ki}$
即 $\tag{19} A(w_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}w_{\frac{n}{2}}^{ki}+w_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}w_{\frac{n}{2}}^{ki}$
至此我们发现原本需要 $n$ 次代入值的等式降低到了 $\frac{n}{2}$ 次，依次递归下去，则我们只需要递归 $log_2{n}$ 次即可在 $\omicron(1)$ 的复杂度下求得式子，即我们求得 $n+1$ 个点值对的复杂度为 $\omicron(nlog_2n)$ ，是可以接受的复杂度。

为了更加严谨的证明，以下过程供还有疑问的读者参考
由于式子（16）可得 $w_n^k=-w_n^{\frac{n}{2}+k}$
则 $\tag{20} A(w_n^{\frac{n}{2}+k})=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}w_{\frac{n}{2}}^{(\frac{n}{2}+k)i}+w_n^{\frac{n}{2}+k}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}w_{\frac{n}{2}}^{(\frac{n}{2}+k)i}=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}w_{\frac{n}{2}}^{ki}-w_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}w_{\frac{n}{2}}^{ki}$
其中求和中的 $w_{\frac{n}{2}}^{(\frac{n}{2}+k)i}$ 直接被替换为 $w_{\frac{n}{2}}^{ki}$ 的原因是，经过平方以后，负号被抵消。
复杂度公式则为 $T(n)=T(\frac{n}{2})+\omicron({n})=\omicron({nlog_2n})$

以上为Cooley-Tukey离散傅里叶变换DFT的思路。

我们回到为什么我们要将多项式补充到 $2^k$ 项的问题上，我们发现，根据Cooley-Tukey 分治的策略，我们需要每次分治的两部分都有同样的项数，则总项数必须为 $2^k$ 项。

4.2 逆离散傅里叶变换IDFT(Inverse Discrete Fourier Transformation)

经过DFT，我们将多项式的系数表达转换为多项式的点值表达。在完成乘法运算以后，我们为了获取系数的变换，需要将多项式的点值表达转换为多项式的系数表达。这时我们使用的方法是逆离散傅里叶变换IDFT，他是DFT的逆。

求解IDFT的过程实际上是一个求解线性方程的问题，给出 $n$ 个线性方程为：
$\tag{21}\begin{cases} a^0(w_n^0)^0+a^1(w_n^0)^1 + ... +a^{n-1}(w_n^0)^{n-1} = A(w_n^0)\\ a^0(w_n^1)^0+a^1(w_n^1)^1 + ... +a^{n-1}(w_n^1)^{n-1} = A(w_n^1)\\...\\a^0(w_n^{n-1})^0+a^1(w_n^{n-1})^1 + ... +a^{n-1}(w_n^{n-1})^{n-1} = A(w_n^{n-1}) \end{cases}$ 矩阵形式如下：
$\tag{22} \begin{bmatrix} (w_n^0)^0 & (w_n^0)^1&...&(w_n^0)^{n-1} \\ (w_n^1)^0 & (w_n^1)^1&...&(w_n^1)^{n-1}\\...\\(w_n^{n-1})^0 & (w_n^{n-1})^1&...&(w_n^{n-1})^{n-1} \end{bmatrix} \begin{bmatrix}a^0\\a^1\\...\\a^{n-1}\end{bmatrix}= \begin{bmatrix}A(w_n^0)\\A(w_n^1)\\...\\A(w_n^{n-1})\end{bmatrix}$
假设上述矩阵为 $\mathbf V$ ，则对于矩阵 $\mathbf D$ 中值 $d_{ij}=-v_n^{ij}$
$\tag{23} D=\begin{bmatrix} (w_n^{-0})^0 & (w_n^{-0})^1&...&(w_n^{-0})^{n-1} \\ (w_n^{-1})^0 & (w_n^{-1})^1&...&(w_n^{-1})^{n-1}\\...\\(w_n^{-(n-1)})^0 & (w_n^{-(n-1)})^1&...&(w_n^{-(n-1)})^{n-1} \end{bmatrix}$
设两个矩阵相乘以后的结果为 $\mathbf E=\mathbf D ·\mathbf V$
$\tag{24} e_{ij} =\sum_{k=0}^{n-1}d_{ik}v_{kj}=\sum_{k=0}^{n-1}w_n^{-ik}w_n^{kj}=\sum_{k=0}^{n-1}w_n^{k(j-i)}$
当 $j=i$ 时， $e_{ij}=n$
当 $j\neq i$ 时， $e_{ij}=\frac{1-(w_n^{j-i})^n}{1-(w_n^{j-i})}=0$
（其中由于 $w_n^{j-i}$ 为 $n$ 次单位根，又因为 $n$ 次单位根的 $n$ 次为1，所以上式成立）
所以 $\mathbf E=n\mathbf I,\mathbf D=\frac{1}{n}\mathbf V^{-1}(其中\mathbf I 为单位矩阵)$
则 $\tag{25} \begin{bmatrix}a^0\\a^1\\...\\a^{n-1}\end{bmatrix}= \frac{1}{n}\begin{bmatrix} (w_n^{-0})^0 & (w_n^{-0})^1&...&(w_n^{-0})^{n-1} \\ (w_n^{-1})^0 & (w_n^{-1})^1&...&(w_n^{-1})^{n-1}\\...\\(w_n^{-(n-1)})^0 & (w_n^{-(n-1)})^1&...&(w_n^{-(n-1)})^{n-1} \end{bmatrix}\begin{bmatrix}A(w_n^0)\\A(w_n^1)\\...\\A(w_n^{n-1})\end{bmatrix}$
则IDFT便是将 $w_n^{i}=w_n^{-i}$ 对结果再做一次DFT，即可获得最后的系数。

回到为什么要有 $n+1$ 取值的原因在于：为了多项式的点值表达转换为多项式的系数表达，我们需要系数矩阵可逆，为了使得系数矩阵可逆，需要满足矩阵满秩，为了满足矩阵满秩，我们需要使得矩阵的行数大于等于列数。则对于有 $n+1$ 个系数的多项式，必须要有 $n+1$ 个以上的取值

四、FFT的实现过程以及FNT

在具体实现FFT的过程中，还需要考虑到对于每一次递归我们如何进行合理的划分。于是这里引入bitreverse算法，也叫做蝴蝶变换。

16个数进行DFT变换划分过程

如图所示，这是16个数进行DFT变换划分的过程。网上的步骤说的很复杂，简单来说，对于DFT的第

i

次划分，二进制第

i

位上为

1

的分为一类，为

0

的分为一类，即可完成。
举例说明：

对于步骤一来说，第1次划分，二进制第1位（即右边第一位）为0的分到左边，为1的分到右边。我们可以发现左边都为偶数，右边都为奇数。（因为奇数二进制第1位为1）
对于步骤二的左半部分，第2次划分，二进制第2位（即右边第二位）为0的分到左边，为1的分到右边。同理，右半部分，第2次划分，也是二进制第2位（即右边第二位）为0的分到左边，为1的分到右边。
同理，对每次的分组进行划分，即可将所有的数字划分到对应的组中。

通过这种划分方法，我们同时还能总结出另外一个规律，对于对于 $2^k$ 个数字中的任意一个位置的数字，假设原本位置为 $x$ ，二进制反转的函数为 $inverse(x)$ ，则最后其所在的位置为 $inverse(x)$ （第一个位置为0），其中 $x$ 为 $k$ 位二进制。
举例说明：
对于 $16（2^4）$ 个数字中的 $14（1110）$ ，则其 $4$ 位二进制的反转为 $0111（7）$ ，则其最后的位置为第 $7$ 位（ps：图中没有继续将所有的数字划分到每组一个，读者可自行划分检验）

这里可以补充一个写法：假如我们将原数组定义为 $A$ ，经过反转后的数组定义为 $B$ ，则 $A[reverse(i)]=B[i]$ 。
又因为如果 $i$ 是偶数，则 $i=(i>>1)<<1$ ，则对于 $reverse(i)=reverse(i>>1)>>1$ ，但考虑到如果是 $i$ 是奇数，则 $((i>>1) << 1) | 1$ ，则对于 $reverse(i)=reverse(i>>1)>>1|(1<<(k-1))$ 其中 $k$ 为满足 $2^k>n$ 的最小值。
综合写可以写成 $reverse(i)=reverse(i>>1)>>1 | (i\&1) << (k-1)$
通过这个写法，我们可以直接写出所有数字经过DFT划分后的结果。

参考：从多项式乘法到快速傅里叶

最后编辑于：2021.03.22 17:16:00

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 203,271评论 5赞 476
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,275评论 2赞 380
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,151评论 0赞 336
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,550评论 1赞 273
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,553评论 5赞 365
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,559评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,924评论 3赞 395
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,580评论 0赞 257
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,826评论 1赞 297
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,578评论 2赞 320
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,661评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,363评论 4赞 318
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,940评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,926评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,156评论 1赞 259
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,872评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,391评论 2赞 342