FFT 快速傅立叶变换

July 29, 2017

C++, ACM

快速傅立叶变换是离散傅立叶变换的加速版。傅立叶的一个用途是用来计算多项式乘法。先讲一下多项式乘法。

多项式乘法 #

多项式有两种表示方法：

系数表示法 #

$$A = a0 + a1x^1 + a2x^2 + … + anx^n$$

点值表示法 #

A ={ (x1, y1), (x2, y2),...,(x3, y3) }

点值表示法相当与用点值来求得各项的系数，所以要准确表示一个多项式就必须至少与多项式未知数的数量一致才可以。

传统的系数表示法求多项式乘法的方式复杂度是O(n^2)。如果使用点值来计算多项式乘法，就会有一些可优化的地方。首先两个n次多项式相乘得到的结果多项式有2n项，这个是确定的，也就是说要用点值表示法就要有2n个点。如何用点值表示法求多项式相乘呢？我们取点的时候就取x相同的点值，这样对应的y值相乘的结果就是结果多项式的y。所以我们在取两个多项是的点的时候就要取2n个点，然后得到2n个结果点。然后剩余的任务就是将这个点值表示法转换成系数表示法得到了我们要的多项式。这个由点值转换为系数的方式称之为插值。

但是这样的算法，在第一步的时候就崩溃了，因为求一个点的值就要至少O(n)，我们要求2n×2个点复杂度就O(n^2)了。然后就需要一些数学手段。

秦九韶算法 #

这个算法的中心思想是通过选取特殊的点来简化计算，什么点是特殊的呢，用复数点。

介绍这个之前需要先复习一下复数。高中数学的基础还是比大学的更有用。复数的基本表示方法：z = x + yi。在复平面上的从原点到点(x,y)的向量是其对应的几何解释。利用直角坐标与极坐标之间关系，复数还有一种三角表示：z = r(cosθ + isinθ)，θ代表极角，r是复数向量的模|z|=sqrt(x^2+y^2)。另一种表示法是指数表示法，见下图 fft-1 . 然后有一个棣莫弗公式能够将两个复数相乘简化：

设两个复数（用三角形式表示）:

        z1 = r1(cosθ1 + isinθ1)
        z2 = r2(cosθ2 + isinθ2)

那么相乘的结果就是：

        z1*z2 = r1*r2[cos(θ1+θ2) + isin(θ1+θ2)]

这个过程很好证明，三角函数就可以了。
这个公式可以进行推广，数学归纳法，得到：
        
        z1*z2*...*zn = r1*r2*...*rn[cos(θ1+θ2+...+θn) + isin(θ1+θ2+...+θn)]

这个公式再变化一下，另z1=z2=...=zn， 就得到复数的幂次公式：
        
        z^n = r1^n[cos(nθ) + isin(nθ)]

然后我们用指数形式换一下:
        
        z^n = r1^n e^(nθi)

然而现在还不够，下面引入复数开根计算。

对于一个复数做开n次根计算会得到n个复数。这n个复数的几何意义： z^(1/n) 的 n 个值是以原点为中心，r^(1/n) 为半径的圆的内接正 n 边形的 n 个顶点。

开根公式：

z^(1/n) = r^(1/n)[cos(θ+2k*pi)/n + isin(θ+2k*pi)/n] , k=0,1,2,...,n-1 z^(1/n) = r^(1/n)* e^(i(θ+2k*pi)/n), k=0,1,2,...,n-1

现在回到我们的主题，在求n次多项式乘法的时候，我们要选取2n个点。这2n个点我们选取2n个复数点。我们令w^(2n) = 1，这2n个复数点我们就选在w的位置上。从上面我们知道了对于复数1 开2n次方根会得到2n个复数的，这些复数就是我们选的x。这些复数也叫做单位复根，因为是复数1+0i这个单位复数的根。

下面再进一段简单的公式规范。

由于我们要计算的结果多项式是2n项，所以我们先将A和B这两个多项式补成2n次多项式，没有的项系数为0就好了。

下面我们令: n = 2n;

对于我们取的单位复根wi，其对应的yi = a0wi^0 + a1wi^1 + ... + a(n-1)wi^(2n)

我们定义：y向量是系数向量a的离散傅立叶变换(DFT)。

这个定义其实就是一个矩阵相乘嘛，多项式值计算的向量表示。我们要做的事情就是快速计算出y向量。因为计算y向量在原来的算法中是个O(n^2)的算法。下面就是快速傅立叶变换(FFT)

快速傅立叶变换的主要思路十分简单，就是分治。由于其使用分治算法，所以其有一个特殊的要求，就是项数必须是2的幂次。之所以有这个要求，是因为这个分治策略是基于下面的公式：

对于多项式
    A(x) = a0 + a1x^1 + a2x^2 + ... + a(n-1)x^(n-1),
    这里保证n是2的幂，也就是总的项数是2的幂，也就是系数向量的长度是2的幂。

现在利用A的奇数位系数和偶数位系数分别构造两个多项式：
    A[0](x) = a0 + a2x^1 + a4x^2 + ... + a(n-2)x^(n/2-1)
    A[1](x) = a1 + a3x^1 + a5x^2 + ... + a(n-1)x^(n/2-1)

那么 A(x) = A[0](x^2) + xA[1](x^2);
这个的证明很简单，不证明了。
这样求A(x) 就转化为求 A[0](x^2) + xA[1](x^2) 。

在上面我们将了这么多就是为了选取特殊的x，那么现在就是要发挥其特殊性的时候了。

再来梳理几点，

x^n = 1 （复数）
x = 1^(1/n) = e^(i2kpi/n)

从上面我们构造的两个式子可以很明显的知道，需要知道x和x^2 之间的神秘关系嘛。

这个时候我们需要在引入一段证明：

我们令：w_n 表示 1的n次复根，即 w_n = 1^(1/n)

我们可以得到一个关系： (w_2n)^2 = 1^(1/n)
即我们得到了一个关系： (w_n) = (w_2n)^2    
也即：
        w_(n/2) = (w_n)^2

上面我们所用的x是1的n次单位复根，即w_n。那么x^2 就是1的(n/2)次单位复根。通过这个关系我们可以得到：x^(n/2) = 1^(1/2) = e^(ipi) = cos(pi)+isin(pi) = -1注意这里指数中有复数，拆回到三角形式计算才是正确的。这里我的证明貌似有些乱，这个引理叫做：相消引理。

上面我们说过，这n个单位复根是正n边形的顶点。而通过上面这个式子，我们又发现我们平方以后将这个正n边形变成了正n/2边形。也就意味着，x的真正的取值数量被我们缩减到一半了。根据分治递归的思维，你知道会一直这么缩减下去，于是复杂度就降低到了O(lgn). 那么是怎么缩减的呢？表现在公式上是什么样子？这里n是2的幂，下面我们令n=4. 由于需要表示单位复根的次数，现在的编辑格式不允许上下标。所以引入一张图，其中：

W的下标表示单位复根的次数，上标表示这是第i个取值。这里通过公式可以发现其重复的顺序是前n/2 与后n/2一致。新的完整公式如下：

其实不通过这个公式，从上面的几何意义中我们也可以找到规律。因为点的值就跟k有关嘛，正n边形缩小为正n/2边形，其角度的分子变成2倍，那么对应的角度就会呈现一轮接一轮的规律。例如原本是正8边形，这八个点分布在8个方向上，角度分别是0，45，90，135，180，… 然后我们将每个的角度都乘以2.，就变成了0，90，180，270，0….

这部分就是折半引理

总结 #

到这里我们介绍了其中的一些数学原理，虽然这个对于写代码来说并没有多少用处。但是知道原理搞起来更爽，不断的探寻世界的真理，如果其他的东西不好探究，至少某些数学的原理还是好探究的。

这只是整个过程的一半呢，从我们计算多项式乘法的步骤来看，

 计算点值（求值）
 点乘
 计算计算系数（插值）

现在我们只是讲完了计算点值的理论还没有写计算点值的分治代码。但是接下来不打算直接写代码，因为理论还没有结束，过程非常巧妙。

点乘的部分就不用说了，就是扫一遍就好了。

插值 #

求值的过程其实就是一个矩阵乘法这个我们说过了, 如下，其中下标两个数分别表示n次单位复根 和第i个单位复根

我们知道这里的n个变量之间是有关系的，所以我们可以进行消掉一些变量。这时候我们引入一个特殊的变量叫做：主单位复根，并使用该复根替代其他的变量。主单位复根是指k取1时的那个复根，值为：W1 = cos(2pi/n) + isin(2pi/n)。这个可以直接计算出来的，所以这个矩阵就变成了一个新的范徳蒙德矩阵。新式子如下：

在引入一个求和引理：

对于任何n>=1 和不能被n整除的非零整数k，有：
SIGMA{i=0,n-1}{(Wn,k)^i} = 0;

这里引理暂时看不懂没有关系。这个引理主要用来求这个矩阵的逆矩阵的。因为要插值可以直接用逆矩阵来搞。

我觉得这部分的证明就没太有意思了，因为我看的证明是特么直接给出结论然后证明的，跟不证明没有什么区别。我们只要知道这个矩阵的逆矩阵很有规律，而且跟原矩阵有很大关系，待会儿写代码的时候会说。

从程序设计的角度理解 #

还是问这个问题，就是快速傅立叶变换到底做什么以及到底做了什么？答：快速傅立叶变换是用来快速求向量卷积的。在本例中就是求多项式乘法。两个向量的卷积仍然是一个向量，长度为原向量长度之和-1。

对于程序设计来说，知道输入和输出，就可以拿封装好的代码拿来用了。