快速傅里叶变换（FFT）算法学习

前言

人生如逆旅，我亦是行人。

---

一、介绍

算法的世界多么广大，我们可以将算法大致分为两类：

• 第一类是较为有用的算法：比如一些经典的图算法，像 DFS 和 BFS（深度 / 广度优先算法），这些算法应用在很多方面，他们非常高效，
  

• 第二类算法是那些极具美感的算法：例如当我们第一次看到汉诺塔的递归实现算法时的状态，你肯定会觉得这个算法贼牛逼，甚至会被它的美感所震撼到，但这些算法或许并没有那么有用或者高效，不过我们还是会研究它，就像没事干了我们还是会看小说；
  

这些伟大而精妙的算法会激发我们的灵感，促使我们发散性思考，下面介绍一种同时属于以上两类的算法：快速傅里叶变换（FFT）。

FFT 算法毫无疑问是上个世纪发明的最伟大的算法之一。很多我们现代生活所依赖的科技，比如无线通讯和 GPS ，以至于任何和信号处理有关的算法，都依赖于 FFT 的深刻洞见。

人们常常不能感知 FFT 的美感，因为它常常有很长很长的前置知识要求， 比如离散傅里叶变换、时域/频域转换、甚至更多。

yysy（有一说一），想要真正搞清 FFT 的应用，你也不得不把这些前置要求都学好。

快速傅里叶变换 ( fast Fourier transform ), 即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称 FFT 。

---

二、学习

• 简单的问题：给定两个多项式，计算二者的乘积。
• 我们的任务就是设计高效的算法解决这个问题。
  

1. 数学上，这个问题可以直接通过重复使用乘法分配律解决。
   
   但是在计算机中，首先需要解决的问题是，如何存储一个多项式？

• 显然，最自然的做法就是储存多项式的系数。我们只需要把系数映射到一个数组中。

系数应该按如下图所示的顺序储存。

因为这样一来数组的第 k 个数字正好对应多项式的 k 阶项系数。

这种表示方法为多项式的系数表示方法。

一般情况下，给定两个 d 阶多项式，二者乘积应为 2d 阶多项式。

并且这一算法，如果用 naive 的乘法分配律来算，时间复杂度为 O(d ^ 2) 。因为多项式 A 的每一项都会跟多项式 B 的每一项相乘。

如何改进这个算法？

---

例子

我们以最简单的多项式，即一阶多项式/线性函数为例。

我们可以用两个系数来表示线性函数，即截距（0阶系数）和斜率（1阶系数）（注：一对系数一对一决定了唯一的一条直线），还有没有别的直线的表示也可以唯一确定一个直线呢？

表示其实很多，我们这里使用直线的两点表示。

• 美丽的几何学告诉我们：两点确定一条直线。

• 事实上：高阶多项式也有类似的性质。 （即：任意 d 阶多项式都由 d+1 个点唯一确定的）

例如下图中的三个点确定了左边这个二次函数：

如果图中有四个点，那么我们就能确定左边这个唯一的三次函数：

接下来证明一下这个结论：

• 假设我们有 p（x） 这个 d 阶多项式，通过了给定的 d+1 个点，我们需要证明：系数是唯一。
• 因此，我们将这 d+1 个点带入多项式的方程，得到 d+1 个方程。线性方程组当然可以写作 矩阵-向量 表示。
  
• 矩阵 M 可逆，只要这些 点 x_0、… 、x_d 不重复（即两两都不同（下面为 范德蒙德矩阵）），只需证明 M 的行列式不为零（范德蒙德行列式）
  
  从而我们知道系数p_0、… 、p_d 是被 d+1 个点唯一确定的，从而多项式也是唯一的。

我们现在有两种表示多项式的方法：

• 第一种是系数表示法；
• 第二种是 d+1 个点表示法，称为值表示法；（利用值表示法，多项式乘法变得不能再简单了）

---

回到之前的问题

乘出来的多项式 C（x）必然为 4 阶，所以需要 5 个点表示。我们现在只需要挑出五个点并计算已有的两个多项式的值，然后把对应每个点的两个值相乘，得到 C 在每个点的函数值。

根据我们前面的学习，五个点可以确定唯一的多项式，

这个多项式的系数表示我们也算出来了，利用值表示法，我们计算多项式乘法的时间一下子从 O（d^2）缩短到了 O（d）。

接下来开始快速的多项式乘法算法了。

---

问：

• 给定两个系数表示的 d 阶多项式，我们已经知道了值表示法计算多项式乘法更快，所以我们可以先计算两个多项式在 2d + 1 个点上的值，然后将函数值一对对地乘起来，从而得到乘出来的多项式的值表示，最后，将值表示转换回系数表示。

• 主要的方法我们大致已经清楚，现在的问题在于还有些步骤我们没搞清楚，
• 事实上我们缺的就是：如何将系数转换成值表示，以及反过来，如何将值表示转换成系数？
• 这一过程，实际上就是 FFT 重要之处。

---

四、求值（evaluation）

让我们先关注一个方向：从系数表示到值表示。这个问题称为求值（evaluation）。

给定 d 阶多项式和 n 个点，我们想计算多项式在这 n 个点上的值（n >= d+1）

• 我们可以采取粗暴的方式：直接在 x 轴上挑取 n 个点，一个一个计算函数值，但是如果这样计算的话，将会变得十分麻烦。 即每个点的计算都是 O(d)，加起来还是 O(nd) >= O(d^2)。这样问题又变得像最初那样复杂，现在的问题是能不能整点更快的？

• 将问题简化一下：考虑一个简单的多项式，比如：P(x) = x ^ 2，在 n=8 个点进行求值。
• 那么问题又来了，我们选那八个点更具有代表性呢？
• 有没有这样的一对点：当你知道一个点的函数值时，另一个点的也可以立刻知道？
• 答案是：互为相反数的数对。（前提是：函数 P 是一个偶函数）

• 如果是 P(x) = x^3 咋办，？
• 事实上这个技巧依旧是对的，不过我们要在 -x 的函数值前面加个负号，

• 所以在这两个偶/奇函数的情形中，我们只需要计算一半数量的点就好了，因为剩下一半从奇偶性中就可以的得到了。

---

再来看看奇偶这个点子对一般的多项式还适不适用，将多项式分成奇/偶函数两个部分，把奇函数部分的 x 提出来一个，那么两个括号里都是两个偶多项式函数。

这样的分解使得我们对于一对相反数计算函数值时，只需要计算其中一个，同时我们可以把两部分都看作 x^2 的函数，你可以看到两个子函数的阶数都降到了 2，比原来的 5 阶不知道高到哪去了。

• 推广这一想法，如果我们有一个 n-1 阶多项式，我们想计算它在 n 个点的函数值，我们可以将多项式分为 奇 / 偶两部分，每部分都只有 n/2 - 1 阶
• 对于这两个只有一半阶的多项式，我们怎么计算他们在对应点的值呢？ 这又是一个求值问题，不过这次我们需要计算我们所给的点的平方的函数值。因为我们一开始取了一对对相反数，所以平方之后正好只剩 n/2 个点了，有没有一点像递归算法。
• 总结：计算多项式 P 在 n 个点上的值，这 n 个点是一对对的相反数。我们把多项式分成（如上所述）的两个部分，从而我们有两个阶为 n/2 - 1 的多项式，每个多项式也只需要求 n/2 个点的值。我们只需要递归地求这两个小多项式的值，利用下方这两个公式带回去，就可以得到原多项式的 n 个的值。最后我们就得到了原多项式的值的表达。

如果上述的一切都行得通，那么我们的时间复杂度仅为 O(nlogn) ，因为两个子问题都是原问题规模的一半，而且我们只需要线性时间来计算原多项式的值。

• （一个问题）在递归步骤：我们假设了每个多项式我们都使用相反数来对其进行求值。注意对子问题而言，每个求值点都是平方数，所以都是正的，递归不成立。
• 如何将新的求值点也搞成相反数对？

• 唯一的方法：就是将这些点都取成复数。这也正是 FFT 算法最奇思妙想的地方。

我们需要专门挑一些复数，使得它们平方之后，依旧是正负成对出现的。

现在又有一个问题来了：该怎么取这些复数呢？

---

• 我们可以通过一个逆向工程来得到答案。

---