快速傅里叶变换及其实现

第1章引言

离散傅里叶变换虽然在数学层面很有用，但其算法的时间复杂度较高，在算法层面并不实用。继而，后续研究者又提出了快速傅里叶变换（Fast Fourier Transform，FFT）算法，这才彻底解决了问题。

多项式乘法早在中学数学中就已经学过，例如：

\[\left( 1+2x+3x^2 \right \left( 4+5x+6x^2 \right =4+13x+28x^2+27x^3+18x^4 \]

3*3的多项式乘法，则需要进行9次乘法计算，和少量的合并同类项计算。也就是说，如果想要计算一次N*N的多项式乘法，则需要进行\(N^2\次乘法计算，和少量的合并同类项计算，即：这种计算多项式乘法的算法的时间复杂度为\(\varTheta(N^2\。

第2章多项式的系数表达与点值表达

想要研究多项式，就需要先把多项式写出来。在本章中，我们研究多项式的两种表达方式：系数表达与点值表达。

2.1 系数表达

这种写法还可以再省略一些：由于每个系数后面的\(x^n\写不写出来都一样，所以可以只写出每一项的系数，并构成一个向量：

\[1+2x+3x^2\,\,\Rightarrow \,\,\left( 1, 2, 3 \right \\ 4+5x+6x^2\,\,\Rightarrow \,\,\left( 4, 5, 6 \right \]

\[\left( 1, 2, 3 \right \otimes \left( 4, 5, 6 \right =\left( 4, 13, 28, 27, 18 \right \]

这是一种全新的向量间的乘法运算，称为：卷积（Convolution），用符号\(\otimes\表示。

2.2 点值表达

我们都知道：两点确定一条直线，而直线是一个包含常数项和一次项的多项式；所以，我们也可以说：两点确定一个一次多项式。那么，三点能不能确定一个二次多项式呢？四点又能不能确定一个三次多项式呢？更多的点呢？

引理1（点值表达的唯一性）：对于任意的n个点构成的点集：\(\left\{ \left( x_0, y_0 \right , \left( x_1, y_1 \right , \left( x_2, y_2 \right , ..., \left( x_{n-1}, y_{n-1} \right \right\}\，如果\(x_0\ne x_1\ne x_2\ne ...\ne x_{n-1}\，则该点集能够唯一确定一个\(n-1\次多项式。

\[\text{设}n-1\text{次多项式：}A\left( x \right =a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1} \\ \text{将}\left( x_0, y_0 \right \text{，}\left( x_1, y_1 \right \text{，}\left( x_2, y_2 \right \text{，}...\text{，}\left( x_{n-1}, y_{n-1} \right \text{依次带入，得方程组：} \\ \begin{cases} a_0+a_1x_0+a_2{x_0}^2+...+a_{n-1}{x_0}^{n-1}=y_0\\ a_0+a_1x_1+a_2{x_1}^2+...+a_{n-1}{x_1}^{n-1}=y_1\\ a_0+a_1x_2+a_2{x_2}^2+...+a_{n-1}{x_2}^{n-1}=y_2\\ ...\\ a_0+a_1x_{n-1}+a_2{x_{n-1}}^2+...+a_{n-1}{x_{n-1}}^{n-1}=y_{n-1}\\ \end{cases} \\ \text{将上式写为矩阵方程形式：} \\ \left[ \begin{matrix} 1& x_0& x_{0}^{2}& \cdots& x_{0}^{n-1}\\ 1& x_1& x_{1}^{2}& \cdots& x_{1}^{n-1}\\ 1& x_2& x_{2}^{2}& \cdots& x_{2}^{n-1}\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& x_{n-1}& x_{n-1}^{2}& \cdots& x_{n-1}^{n-1}\\ \end{matrix} \right] \left[ \begin{array}{c} a_0\\ a_1\\ a_2\\ \cdots\\ a_{n-1}\\ \end{array} \right] =\left[ \begin{array}{c} y_0\\ y_1\\ y_2\\ \cdots\\ y_{n-1}\\ \end{array} \right] \\ \text{最左边的矩阵为范德蒙德矩阵，且由于任意的}x_i\ne x_j\text{，所以其行列式：} \\ \prod_{0\leqslant i<j\leqslant n-1}{\left( x_j-x_i \right}\ne 0 \\ \text{所以，矩阵方程存在唯一解，得证。} \]

也就是说，我们不仅可以用n个系数唯一的表达一个n-1次多项式，还可以用这个多项式上的n个点唯一表达。

已知一个多项式的点值表达，怎么得到其系数表达？
点值表达存在的意义是什么？

第一个问题很简单，我们已经知道，想要得到一个多项式的点值表达，需要满足两个条件：

每个点的x互不相同

第二个问题也不难解决，不是有很多点吗，那就用这些点解上面证明过程中的那个线性方程组，解出来的\(a_0, a_1, a_2, ..., a_{n-1}\就是系数表达了。

事实上，我们在这里犯了一个错误：两个二次多项式做乘法，结果应当是一个四次多项式，而四次多项式需要用五个点才能唯一表示，而上面只有三个点，显然是不够的。这个问题提醒了我们：在计算多项式乘法的时候，不能只看现在这个多项式是几次的，还应当看乘法的结果多项式是几次的，后者的次数才能决定一开始要取几个点。

先看从系数表达转点值表达的过程，想做这件事，就要选n个点，每个点依次带入多项式算一遍，而多项式里面全是各种高次幂，时间复杂度早已不可接受了。如果使用霍纳法则，将计算一个多项式的时间复杂度优化到\(\varTheta(N\，那最终也是一个时间复杂度为\(\varTheta(N^2\的算法。

\[A\left( x \right =\sum_{i=0}^{n-1}{\left( y_i\frac{\prod_{j\ne i}{\left( x-x_j \right}}{\prod_{j\ne i}{\left( x_i-x_j \right}} \right} \]

不难看出，这个公式的时间复杂度为\(\varTheta(N^2\。

所以，现在的目标是：找到一对算法，能够更快的在多项式的系数表达和点值表达之间进行转换。

第3章欧拉公式

欧拉公式是一个非常有名的公式，其将复数域上的指数函数和三角函数联系在了一起。让我们从\(e^{i\theta}\的麦克劳林级数展开开始：

\[e^{i\theta}=\sum_{n=0}^{\infty}{\frac{\left( i\theta \right ^n}{n!}}=1+i\theta -\frac{\theta ^2}{2!}-\frac{i\theta ^3}{3!}+\frac{\theta ^4}{4!}+\frac{i\theta ^5}{5!}-... \\ \text{将上式按奇偶顺序重新排列并整理得到：} \\ =\left( 1-\frac{\theta ^2}{2!}+\frac{\theta ^4}{4!}-... \right +i\left( \theta -\frac{\theta ^3}{3!}+\frac{\theta ^5}{5!}-... \right \\ =\sum_{n=0}^{\infty}{\frac{\left( -1 \right ^n}{\left( 2n \right !}\theta ^{2n}+i\sum_{n=0}^{\infty}{\frac{\left( -1 \right ^n}{\left( 2n+1 \right !}\theta ^{2n+1}}} \\ \text{上式的左右两边分别是}cos\theta \text{和}sin\theta \text{的麦克劳林级数：} \\ =\cos \theta +i\sin \theta \]

引理2（欧拉公式）：

\[e^{i\theta}=\cos \theta +i\sin \theta \]

引理3（欧拉恒等式）：

\[e^{\pi i}=-1 \]

第4章 n次单位复数根

这一章中，我们研究n次单位复数根。什么是n次单位复数根呢？其指的是以下n次方程的根：

\[\omega ^n=1 \]

根据欧拉公式，我们可以凑出n次单位复数根的一般形式：

引理4（n次单位复数根）：

\[\text{方程：}\omega ^n=1\text{的所有根依次为：}\left( e^{\frac{2\pi i}{n}} \right ^k, k=0, 1, 2, ..., n-1\text{；记作：}\omega _{n}^{k} \]

\[\text{将}\left( e^{\frac{2\pi i}{n}} \right ^k\text{带入原方程得：} \\ \left( e^{\frac{2\pi i}{n}} \right ^{kn}=e^{2\pi ik} \\ \text{由欧拉公式：} \\ e^{2\pi ik}=\cos 2\pi k+i\sin 2\pi k=1 \\ 得证。 \]

由复数在复平面上的极坐标表示可知：各个n次单位复数根是复平面上单位圆的各个n等分点，故其具有一些特殊的性质。请看：

引理5（折半引理）：\(\omega _{n}^{k}=\omega _{n/2}^{k/2}\

\[\omega _{n}^{k}=e^{\frac{2\pi ik}{n}}=e^{\frac{2\pi i\frac{k}{2}}{\frac{n}{2}}}=\omega _{n/2}^{k/2} \]

引理6：\(\omega _{n}^{n}=1\

\[\omega _{n}^{n}=e^{\frac{2\pi in}{n}}=e^{2\pi i}=\left( -1 \right ^2=1 \]

引理7：\(\omega _{n}^{n/2}=-1\

\[\omega _{n}^{n/2}=e^{\frac{2\pi i\frac{n}{2}}{n}}=e^{\pi i}=-1 \]

引理8：\(\omega _{n}^{-k}=\overline{\omega _{n}^{k}}\

\[\because \omega _{n}^{k}=e^{\frac{2\pi ik}{n}}=\cos \left( \frac{2\pi k}{n} \right +i\sin \left( \frac{2\pi k}{n} \right \\ \text{又}\because \omega _{n}^{-k}=e^{\frac{-2\pi ik}{n}}=\cos \left( -\frac{2\pi k}{n} \right +i\sin \left( -\frac{2\pi k}{n} \right =\cos \left( \frac{2\pi k}{n} \right -i\sin \left( \frac{2\pi k}{n} \right \\ \therefore \omega _{n}^{-k}=\overline{\omega _{n}^{k}} \]

n次单位复数根及其若干引理的用途将在后续章节展开。

第5章离散傅里叶变换与离散傅里叶逆变换

使用所有的n次单位复数根将一个n-1次多项式从系数表达转为点值表达的过程，就称为离散傅里叶变换（Discrete Fourier Transform，DFT）。

那么，离散傅里叶逆变换又是什么呢？顾名思义，其是一个变换回来的过程：

将离散傅里叶变换得到的点值表达转为系数表达的过程，就称为离散傅里叶逆变换（Inverse Discrete Fourier Transform，IDFT）。

事实上，这两个猜想都是正确的。其算法分别被称为：快速傅里叶变换（Fast Fourier Transform，FFT）和快速傅里叶逆变换（Inverse Fast Fourier Transform，IFFT）。在下面的两章中，我们分别研究这两种算法。

第6章快速傅里叶变换

6.1 快速傅里叶变换的数学原理

既然是一种优化算法，那就不能一上来就带入点进行计算，我们需要先做一些准备：

将多项式按奇偶顺序重排并整理：

\[\text{设}n-1\text{次多项式：}A\left( x \right =a_0+a_1x+a_2x^2+a_3x^3...a_{n-1}x^{n-1} \\ \text{按奇偶顺序重排此多项式，得：} \\ A\left( x \right =\left( a_0+a_2x^2+...+a_{n-2}x^{n-2} \right +x\left( a_1+a_3x^2+...+a_{n-1}x^{n-2} \right \]

将多项式分解为两个子多项式：

\[\text{对于：}A\left( x \right =\left( a_0+a_2x^2+...+a_{n-2}x^{n-2} \right +x\left( a_1+a_3x^2+...+a_{n-1}x^{n-2} \right \\ \text{设：}A^{\left[ 0 \right]}\left( x \right =a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1} \\ \text{且：}A^{\left[ 1 \right]}\left( x \right =a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1} \\ \text{则：}A\left( x \right =A^{\left[ 0 \right]}\left( x^2 \right +xA^{\left[ 1 \right]}\left( x^2 \right \]

现在，我们将n次单位复数根，即\(\omega _{n}^{k}\带入上式：

\[\text{对于：}A\left( x \right =A^{\left[ 0 \right]}\left( x^2 \right +xA^{\left[ 1 \right]}\left( x^2 \right \\ \text{将}x=\omega _{n}^{k}\text{带入得：} \\ A\left( \omega _{n}^{k} \right =A^{\left[ 0 \right]}\left( \omega _{n}^{2k} \right +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n}^{2k} \right \]

\[\text{对于：}A\left( \omega _{n}^{k} \right =A^{\left[ 0 \right]}\left( \omega _{n}^{2k} \right +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n}^{2k} \right \\ 1. \text{当}k<\frac{n}{2}\text{时：} \\ A\left( \omega _{n}^{k} \right =A^{\left[ 0 \right]}\left( \omega _{n}^{2k} \right +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n}^{2k} \right \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right \left( \text{由折半引理} \right \\ 2. \text{当}k\geqslant \frac{n}{2}\text{时：} \\ \text{设：}k=k^{\prime}+n/2 \\ \text{则：}A\left( \omega _{n}^{k} \right =A\left( \omega _{n}^{k^{\prime}+n/2} \right \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{2k^{\prime}+n} \right +\omega _{n}^{k^{\prime}+n/2}A^{\left[ 1 \right]}\left( \omega _{n}^{2k^{\prime}+n} \right \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{2k^{\prime}} \right +\omega _{n}^{k^{\prime}+n/2}A^{\left[ 1 \right]}\left( \omega _{n}^{2k^{\prime}} \right \left( \text{由引理}6 \right \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{2k^{\prime}} \right -\omega _{n}^{k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n}^{2k^{\prime}} \right \left( \text{由引理}7 \right \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k^{\prime}} \right -\omega _{n}^{k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k^{\prime}} \right \left( \text{由折半引理} \right \]

对比这两种情况得到的两个结果，可以发现：这两个结果唯一的差别就是一个正负号。也就是说，每当在\(k<n/2\这半边算出一个点，就可以立即在\(k+n/2\处再算出一个点。这样一来，计算量直接减少一半。这还没完，在计算\(A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right\和\(A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right\的时候，这两个多项式不仅长度减半，而且可以继续用这个性质，计算量继续减少一半。这个过程什么时候停下来呢？不难想象：当一个只有两项的多项式被拆成两个只有一项的多项式时，这个过程就停了，因为只有一项的多项式就是一个常数，其值不仅不需要计算，而且与自变量的取值无关。

6.2 一个手工计算快速傅里叶变换的实例

上一节中，我们研究了快速傅里叶变换的数学原理，但如果只看公式的话，读者可能仍然不理解这个算法是怎么进行的。所以，这一节中，我们就通过一个例子，来真正计算一次快速傅里叶变换。

按照快速傅里叶变换的算法流程，我们首先要做的是：将这个多项式的项数补齐至2的整数次幂。由于这个多项式的项数是4，其已经是2的整数次幂了，所以，不需要补齐。

\[A\left( x \right =1+2x+3x^2+4x^3 \\ \text{设：}A^{\left[ 0 \right]}\left( x \right =1+3x \\ \text{且：}A^{\left[ 1 \right]}\left( x \right =2+4x \\ \text{则：}A\left( x \right =A^{\left[ 0 \right]}\left( x^2 \right +xA^{\left[ 1 \right]}\left( x^2 \right \]

此时，我们已经将一个四项的多项式分解为两个两项的多项式了。但这还不够，我们还需要将这两个多项式继续分解为四个只有一项的多项式：

\[\text{对于：}A^{\left[ 0 \right]}\left( x \right =1+3x \\ \text{设：}A^{\left[ 0 \right] \left[ 0 \right]}\left( x \right =1 \\ \text{且：}A^{\left[ 0 \right] \left[ 1 \right]}\left( x \right =3 \\ \text{则：}A^{\left[ 0 \right]}\left( x \right =A^{\left[ 0 \right] \left[ 0 \right]}\left( x^2 \right +xA^{\left[ 0 \right] \left[ 1 \right]}\left( x^2 \right \\ \text{对于：}A^{\left[ 1 \right]}\left( x \right =2+4x \\ \text{设：}A^{\left[ 1 \right] \left[ 0 \right]}\left( x \right =2 \\ \text{且：}A^{\left[ 1 \right] \left[ 1 \right]}\left( x \right =4 \\ \text{则：}A^{\left[ 1 \right]}\left( x \right =A^{\left[ 1 \right] \left[ 0 \right]}\left( x^2 \right +xA^{\left[ 1 \right] \left[ 1 \right]}\left( x^2 \right \]

与输入无关，而分别恒等于其常数项的值。

\[\text{对于：}A^{\left[ 0 \right]}\left( x \right =A^{\left[ 0 \right] \left[ 0 \right]}\left( x^2 \right +xA^{\left[ 0 \right] \left[ 1 \right]}\left( x^2 \right \\ \text{要得到}A^{\left[ 0 \right]}\left( x \right \text{的}FFT\text{，我们就需要将}\omega _{2}^{0}\text{和}\omega _{2}^{1}\text{带入}A^{\left[ 0 \right]}\left( x \right \text{；} \\ \text{由快速傅里叶变换的性质，我们实际上要做的是：} \\ 1. \text{计算}\alpha =A^{\left[ 0 \right] \left[ 0 \right]}\left( \omega _{1}^{0} \right \\ 2. \text{计算}\beta =\omega _{2}^{0}A^{\left[ 0 \right] \left[ 1 \right]}\left( \omega _{1}^{0} \right \\ 3. A^{\left[ 0 \right]}\left( \omega _{2}^{0} \right =\alpha +\beta \text{；}A^{\left[ 0 \right]}\left( \omega _{2}^{1} \right =\alpha -\beta \\ \text{由于}A^{\left[ 0 \right] \left[ 0 \right]}\left( \omega _{1}^{0} \right \text{，}A^{\left[ 0 \right] \left[ 1 \right]}\left( \omega _{1}^{0} \right \text{的值与输入无关，就是}1\text{和}3\text{，且}\omega _{2}^{0}=1\text{，所以：} \\ \alpha =1\text{；}\beta =3\text{；}A^{\left[ 0 \right]}\left( \omega _{2}^{0} \right =4\text{；}A^{\left[ 0 \right]}\left( \omega _{2}^{1} \right =-2 \]

同理：

\[\text{对于：}A^{\left[ 1 \right]}\left( x \right =A^{\left[ 1 \right] \left[ 0 \right]}\left( x^2 \right +xA^{\left[ 1 \right] \left[ 1 \right]}\left( x^2 \right \\ \text{要得到}A^{\left[ 1 \right]}\left( x \right \text{的}FFT\text{，我们就需要将}\omega _{2}^{0}\text{和}\omega _{2}^{1}\text{带入}A^{\left[ 1 \right]}\left( x \right \text{；} \\ \text{由快速傅里叶变换的性质，我们实际上要做的是：} \\ 1. \text{计算}\alpha =A^{\left[ 1 \right] \left[ 0 \right]}\left( \omega _{1}^{0} \right \\ 2. \text{计算}\beta =\omega _{2}^{0}A^{\left[ 1 \right] \left[ 1 \right]}\left( \omega _{1}^{0} \right \\ 3. A^{\left[ 1 \right]}\left( \omega _{2}^{0} \right =\alpha +\beta \text{；}A^{\left[ 1 \right]}\left( \omega _{2}^{1} \right =\alpha -\beta \\ \text{由于}A^{\left[ 1 \right] \left[ 0 \right]}\left( \omega _{1}^{0} \right \text{，}A^{\left[ 1 \right] \left[ 1 \right]}\left( \omega _{1}^{0} \right \text{的值与输入无关，就是}2\text{和}4\text{，且}\omega _{2}^{0}=1\text{，所以：} \\ \alpha =2\text{；}\beta =4\text{；}A^{\left[ 1 \right]}\left( \omega _{2}^{0} \right =6\text{；}A^{\left[ 1 \right]}\left( \omega _{2}^{1} \right =-2 \]

\[\text{对于：}A\left( x \right =A^{\left[ 0 \right]}\left( x^2 \right +xA^{\left[ 1 \right]}\left( x^2 \right \\ \text{要得到}A\left( x \right \text{的}FFT\text{，我们就需要将}\omega _{4}^{0}\text{，}\omega _{4}^{1}\text{，}\omega _{4}^{2}\text{，}\omega _{4}^{3}\text{带入}A\left( x \right \text{；} \\ \text{由快速傅里叶变换的性质，我们实际上要做的是：} \\ 1.\text{计算}\alpha _0=A^{\left[ 0 \right]}\left( \omega _{2}^{0} \right \\ 2.\text{计算}\beta _0=\omega _{4}^{0}A^{\left[ 1 \right]}\left( \omega _{2}^{0} \right \\ 3.\text{计算}\alpha _1=A^{\left[ 0 \right]}\left( \omega _{2}^{1} \right \\ 4.\text{计算}\beta _1=\omega _{4}^{1}A^{\left[ 1 \right]}\left( \omega _{2}^{1} \right \\ 5.A\left( \omega _{4}^{0} \right =\alpha _0+\beta _0\text{；}A\left( \omega _{4}^{2} \right =\alpha _0-\beta _0 \\ 6.A\left( \omega _{4}^{1} \right =\alpha _1+\beta _1\text{；}A\left( \omega _{4}^{3} \right =\alpha _1-\beta _1 \\ \text{而：}A^{\left[ 0 \right]}\left( \omega _{2}^{0} \right =4\text{，}A^{\left[ 1 \right]}\left( \omega _{2}^{0} \right =6\text{，}A^{\left[ 0 \right]}\left( \omega _{2}^{1} \right =-2\text{，}A^{\left[ 1 \right]}\left( \omega _{2}^{1} \right =-2\text{，都是已经计算好的} \\ \text{且}\omega _{4}^{0}=1\text{，}\omega _{4}^{1}=e^{\frac{2\pi i}{4}}=\cos \frac{\pi}{2}+i\sin \frac{\pi}{2}=i \\ \text{所以：}\alpha _0=4\text{；}\beta _0=6\text{；}\alpha _1=-2\text{；}\beta _1=-2i \\ \text{所以：}A\left( \omega _{4}^{0} \right =10\text{；}A\left( \omega _{4}^{2} \right =-2\text{；}A\left( \omega _{4}^{1} \right =-2-2i\text{；}A\left( \omega _{4}^{3} \right =-2+2i \\ \text{又：}\omega _{4}^{0}=1\text{；}\omega _{4}^{1}=i\text{；}\omega _{4}^{2}=\left( \omega _{4}^{1} \right ^2=-1\text{；}\omega _{4}^{3}=\left( \omega _{4}^{1} \right ^3=-i \\ \text{所以，}A\left( x \right \text{的}FFT\text{为：} \\ \left\{ \left( 1,10 \right ,\left( i,-2-2i \right ,\left( -1,-2 \right ,\left( -i,-2+2i \right \right\} \]

至此，我们就完成了一次快速傅里叶变换的计算。读者可以将4个自变量依次带入多项式，来验证结果的正确性。

6.3 过程更简略的手工计算实例

这一节中，我们再重新算一次上面这个多项式的快速傅里叶变换；但这一次，我们省去一切不必要的过程，只保留真正需要计算的部分，看看是什么体验。

1 2 3 4

这里的1 2 3 4表示的就是上面的\(A\left( x \right =1+2x+3x^2+4x^3\

接下来，将系数按奇偶顺序分组：

1   2   3   4
1   3 | 2   4
1 | 3 | 2 | 4

这里的1 | 3 | 2 | 4依次对应着上面的\(A^{\left[ 0 \right] \left[ 0 \right]}\left( x^2 \right \text{，}A^{\left[ 0 \right] \left[ 1 \right]}\left( x^2 \right \text{，}A^{\left[ 1 \right] \left[ 0 \right]}\left( x^2 \right \text{，}A^{\left[ 1 \right] \left[ 1 \right]}\left( x^2 \right\，由于这四个多项式的值与输入无关，所以其值分别就是1、3、2、4。

以连续的两个系数为一组，每一组中，相邻的两个系数之间做一对计算，需要用到系数：\(\omega _{2}^{0}=1\：

1 3 | 2 4

1 + 1 * 3 = 4
1 - 1 * 3 = -2
2 + 1 * 4 = 6
2 - 1 * 4 = -2

4 -2 6 -2

第二次倒推以连续的四个系数为一组，每一组中，第n个系数和第n+2个系数之间做一对计算，需要用到系数：\(\omega _{4}^{0}=1\text{；}\omega _{4}^{1}=i\：

4 -2 6 -2

4 + 1 * 6 = 10
4 - 1 * 6 = -2
-2 + i * -2 = -2-2i
-2 - i * -2 = -2+2i

10 -2-2i -2 -2+2i

至此，快速傅里叶变换就计算完成了（n次单位复数根的值不需要算出）。由此可见，一旦读者熟练掌握了快速傅里叶变换的计算原理，就可以使用这种非常简洁的计算过程进行快速傅里叶变换的计算了。

第7章快速傅里叶逆变换

7.1 快速傅里叶逆变换的数学原理

对于多项式：

\[A\left( x \right =\sum_{j=0}^{n-1}{a_jx^j} \]

\[y_i=A\left( \omega _{n}^{i} \right =\sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{i} \right ^j} \]

此时，我们需要构造一个新的多项式，这个多项式的系数由\(y_i\给出，并将单位复数根全部取倒数后带入：

\[\text{设：}B\left( x \right =\sum_{i=0}^{n-1}{y_ix^i} \\ z_k=B\left( \omega _{n}^{-k} \right =\sum_{i=0}^{n-1}{y_i\left( \omega _{n}^{-k} \right ^i} \]

\[z_k=B\left( \omega _{n}^{-k} \right =\sum_{i=0}^{n-1}{y_i\left( \omega _{n}^{-k} \right ^i} \\ =\sum_{i=0}^{n-1}{\left( \left( \sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{i} \right ^j} \right \left( \omega _{n}^{-k} \right ^i \right} \\ =\sum_{i=0}^{n-1}{\left( \left( \sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{j} \right ^i} \right \left( \omega _{n}^{-k} \right ^i \right} \\ =\sum_{i=0}^{n-1}{\sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{j} \right ^i\left( \omega _{n}^{-k} \right ^i}} \\ =\sum_{i=0}^{n-1}{\sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{j-k} \right ^i}} \\ =\sum_{j=0}^{n-1}{a_j\sum_{i=0}^{n-1}{\left( \omega _{n}^{j-k} \right ^i}} \]

此时，需要分两种情况讨论：

\[1.\text{当}j=k\text{时：} \\ \text{此时，}\left( \omega _{n}^{j-k} \right ^i=1^i=1 \\ \text{则：}a_k\sum_{i=0}^{n-1}{1}=na_k \\ 2.\text{当}j\ne k\text{时：} \\ \text{此时，}\sum_{i=0}^{n-1}{\left( \omega _{n}^{j-k} \right ^i}\text{是一个首项为}\left( \omega _{n}^{j-k} \right ^0=1\text{，公比为}\omega _{n}^{j-k}\text{的等比数列，由等比数列求和公式得：} \\ \sum_{j=0}^{n-1}{a_j\sum_{i=0}^{n-1}{\left( \omega _{n}^{j-k} \right ^i}} \\ =\sum_{j=0}^{n-1}{a_j\cdot \frac{1\left( 1-\left( \omega _{n}^{j-k} \right ^n \right}{1-\omega _{n}^{j-k}}} \\ =\sum_{j=0}^{n-1}{a_j\cdot \frac{1-1^{j-k}}{1-\omega _{n}^{j-k}}}\left( \text{由引理}6 \right \\ =\sum_{j=0}^{n-1}{a_j\cdot 0} \\ =0 \\ \text{综上：}z_k=na_k\text{，即：}a_k=\frac{z_k}{n} \]

为了研究这个问题，让我们回到快速傅里叶变换推导过程的起点：

\[A\left( x \right =A^{\left[ 0 \right]}\left( x^2 \right +xA^{\left[ 1 \right]}\left( x^2 \right \]

\[1.\text{当}k<\frac{n}{2}\text{时：} \\ A\left( \omega _{n}^{k} \right =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right \\ 2.\text{当}k\geqslant \frac{n}{2}\text{时：} \\ A\left( \omega _{n}^{k} \right =A\left( \omega _{n}^{k^{\prime}+n/2} \right =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k^{\prime}} \right -\omega _{n}^{k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k^{\prime}} \right \]

即：只需要计算\(A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right\和\(\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right\就行了。

\[\text{对于：}A\left( \omega _{n}^{-k} \right =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k} \right +\omega _{n}^{-k}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k} \right \\ 1.\text{当}k<\frac{n}{2}\text{时：} \\ A\left( \omega _{n}^{-k} \right =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k} \right +\omega _{n}^{-k}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k} \right \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{-k} \right +\omega _{n}^{-k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{-k} \right \left( \text{由折半引理} \right \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{-k} \right +\overline{\omega _{n}^{k}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{-k} \right \left( \text{由引理}8 \right \\ 2.\text{当}k\geqslant \frac{n}{2}\text{时：} \\ \text{设}k=k^{\prime}+n/2 \\ \text{则：}A\left( \omega _{n}^{-k} \right =A\left( \omega _{n}^{-k^{\prime}-n/2} \right \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k^{\prime}-n} \right +\omega _{n}^{-k^{\prime}-n/2}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k^{\prime}-n} \right \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k^{\prime}} \right +\omega _{n}^{-k^{\prime}+n/2}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k^{\prime}} \right \left( \text{由引理}6 \right \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k^{\prime}} \right -\omega _{n}^{-k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k^{\prime}} \right \left( \text{由引理}7 \right \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{-k^{\prime}} \right -\omega _{n}^{-k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{-k^{\prime}} \right \left( \text{由折半引理} \right \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{-k^{\prime}} \right -\overline{\omega _{n}^{k^{\prime}}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{-k^{\prime}} \right \left( \text{由引理}8 \right \]

可以发现：如果将单位复数根的倒数带入，那么在进行快速傅里叶变换的过程中，只有一个系数的差别（\(A^{\left[ 0 \right]}\和\(A^{\left[ 1 \right]}\括号里面的差别可以忽略，因为分解后的多项式最终将与输入无关）。也就是说，只需要对快速傅里叶变换的计算过程做两处微小的改动，我们就可以再一次利用快速傅里叶变换去计算所有\(z_k\的值，从而计算出所有\(a_k\的值了。这两处改动分别为：

将每个快速傅里叶变换的结果除以n

7.2 一个手工计算快速傅里叶逆变换的实例

上一节中，我们研究了快速傅里叶逆变换的数学原理。这一节中，我们就使用上一章得到的快速傅里叶变换的计算结果10 -2-2i -2 -2+2i，进行一次快速傅里叶逆变换的手工计算，如果计算结果能够回到1 2 3 4，就说明我们的计算是正确的。

接下来，对系数进行分组：

10   -2-2i   -2      -2+2i
10   -2    | -2-2i   -2+2i

接下来，进行第一次倒推，需要用到系数：\(\overline{\omega _{2}^{0}}=1\：

10 -2 | -2-2i -2+2i

10 + 1 * -2 = 8
10 - 1 * -2 = 12
-2-2i + 1 * -2+2i = -4
-2-2i - 1 * -2+2i = -4i

8 12 -4 -4i

接下来，进行第二次倒推，需要用到系数：\(\overline{\omega _{4}^{0}}=1\text{；}\overline{\omega _{4}^{1}}=-i\：

8 12 -4 -4i

8 + 1 * -4 = 4
8 - 1 * -4 = 12
12 + -i * -4i = 8
12 - -i * -4i = 16

4 8 12 16

上面得到的4 8 12 16，分别就是\(z_0, z_1, z_2, z_3\。最后，我们将其都除以\(n=4\，就能得到\(a_0, a_1, a_2, a_3\了：

可见，计算结果完全符合预期。

第8章快速傅里叶变换的实现

8.1 对齐至2的整数次幂的算法

快速傅里叶变换的第一步是将多项式的项数对齐至2的整数次幂，所以，我们需要根据输入的项数，来找到需要对齐到的项数。这一需求的朴素算法是使用一个循环，并使用一个从1开始，不断自乘2的数字和输入项数作比较，直至这个数字已经大于等于输入项数时，算法终止。这个算法很简单，读者可以自行尝试。

unsigned __nextPow2(unsigned N
{
    N--;

    N |= N >> 1;
    N |= N >> 2;
    N |= N >> 4;
    N |= N >> 8;
    N |= N >> 16;

    return N + 1;
}

这个算法不难理解，其要点在于：

0b101

0b111

0b1000

0b101

0b111

最高n位

如果N已经是2的整数次幂，那么，算法直接返回N就行了。但是，判断一个数字是不是2的整数次幂需要额外的代价，且很明显，这个判定的失败率是很高的，因为绝大多数的整数都不是2的整数次幂。所以，干脆就不要判定这件事了，而是将N减去1。如果N是2的整数次幂，减去1后就会丢失其最高位的1，并变成一个全是1的数字，在经过多次（无用的）位或运算后，又被加上1，回到了原值。例如0b1000，减去1后会变成0b111，其丢失了原数字最高位的1，最终，0b111又会因为加1而回到0b1000并返回。另一方面，如果N并不是2的整数次幂，那就说明N除了最高位的1以外，在低位还有1，这样一来，减去1就不会使N丢失最高位的1，所以，其结果不受影响

8.2 位逆序算法

仔细观察分组前后，各个系数的索引值的二进制表示，这里以8个系数为例：

000 001 010 011 100 101 110 111  // 分组前
000 100 010 110 001 101 011 111  // 分组后

不难发现：分组前后的每一对索引值都是位逆序的。

那么，怎么实现位逆序呢？朴素的算法是：通过一个循环，将待转换的数字不断右移1位，同时将转换后的数字不断左移1位，并将两个数字的最低位对接即可。读者可以自行尝试。

unsigned __bitReverse(unsigned N, unsigned bitWidth
{
    N = ((0xaaaaaaaa & N >> 1 | ((0x55555555 & N << 1;
    N = ((0xcccccccc & N >> 2 | ((0x33333333 & N << 2;
    N = ((0xf0f0f0f0 & N >> 4 | ((0x0f0f0f0f & N << 4;
    N = ((0xff00ff00 & N >> 8 | ((0x00ff00ff & N << 8;

    N = ((N >> 16 | (N << 16 >> (32 - bitWidth;

    return N;
}

这个算法不难理解，其要点在于：

0xaaaaaaaa

0b1010...

0x55555555

0b0101

类似的，0xcccccccc是形如0b1100...的位掩码；而0x33333333是形如0b0011的位掩码；在这两个位掩码，以及后续的左右移位和位或的作用下，N会以每两位为一组发生交换。以此类推，N又会以每4、8位为一组进行交换
最后，N需要以每16位为一组，完成最后一次交换，此时就不需要位掩码了，直接交换即可。至此，N的所有位完成了逆序
上述算法完成的是32位无符号整数的位逆序，而实际输入的数字很可能并没有这么多位。例如：0b001的位逆序应该是0b100，而按照上述算法，最终的结果是：0b100...(后面还有29个0，这不是我们需要的。所以，最后还需要做一次右移，将多余的0去掉

8.3 快速傅里叶变换与快速傅里叶逆变换的实现

快速傅里叶逆变换需要在最后对所有的结果除以n

我们可以使用一个要么是1，要么是-1的数字来同时区别这两处不同。请看：

vector<complex<double>> __FFT(const vector<complex<double>> &coefList, double conjNum
{
    vector<complex<double>> FFTList(coefList.size(;

    unsigned bitWidth = __builtin_ctz(FFTList.size(;

    for (unsigned idx = 0; idx < FFTList.size(; idx++
    {
        FFTList[idx] = coefList[__bitReverse(idx, bitWidth];
    }

    for (unsigned N = 2; N <= FFTList.size(; N *= 2
    {
        for (unsigned startIdx = 0; startIdx < FFTList.size(; startIdx += N
        {
            complex<double> curOmega(1.;
            complex<double> mulOmega(cos(2 * M_PI / N, conjNum * sin(2 * M_PI / N;

            for (unsigned leftIdx = startIdx, rightIdx = startIdx + N / 2; leftIdx < startIdx + N / 2; leftIdx++, rightIdx++
            {
                auto leftNum  = FFTList[leftIdx] + curOmega * FFTList[rightIdx];
                auto rightNum = FFTList[leftIdx] - curOmega * FFTList[rightIdx];

                FFTList[leftIdx]  = leftNum;
                FFTList[rightIdx] = rightNum;

                curOmega *= mulOmega;
            }
        }
    }

    if (conjNum == -1.
    {
        for (auto &FFTNum: FFTList
        {
            FFTNum /= FFTList.size(;
        }
    }

    return FFTList;
}

__FFT函数用于计算快速傅里叶变换以及快速傅里叶逆变换。当conjNum = 1.时，其处于快速傅里叶变换模式；而当conjNum = -1.时，其处于快速傅里叶逆变换模式（由于这个函数不作为对外接口，所以没有对conjNum使用布尔值或枚举变量等编程手段限制其他错误的值，读者如果对此感到介意，可以自行实现一个更严谨的接口）。
coefList为系数列表，所有的系数已经由主调函数从double类型转为了complex<double>类型，并已经进行了对齐处理；conjNum已在上文中说明，其只会传入1.或-1.。
__builtin_ctz完成的，其返回输入数字从最低位到第一个1之间的0的数量。分组操作通过循环进行，其将coefList中的系数重排至FFTList列表中。
第一重循环用于遍历N的取值，N从2开始，以不断自乘2的方式递增，直至与多项式的项数一致时终止。
startIdx存放的是当前分组的起始索引值；而分组的长度（决定了startIdx的循环增量）是恒等于N的。比如，第一次倒推时，以两个数字为一组；第二次倒推时，以四个数字为一组；以此类推。
curOmega变量用于存放\(\omega _{n}^{k}\或\(\overline{\omega _{n}^{k}}\的当前值，其被初始化为1；而mulOmega变量用于存放\(\omega _{n}^{1}\或\(\overline{\omega _{n}^{1}}\，其用于在每一对系数计算完成后，将curOmega自乘一次mulOmega。
mulOmega的值基于欧拉公式：
\[\omega _{n}^{1}=e^{\frac{2\pi i}{n}}=\cos \frac{2\pi}{n}+i\sin \frac{2\pi}{n}
\\
\overline{\omega _{n}^{1}}=e^{\frac{2\pi i}{n}}=\cos \frac{2\pi}{n}-i\sin \frac{2\pi}{n}
\]

conjNum用于控制上式中的正负号。
leftIdx维护，初始化为startIdx，即当前分组的起始索引值；右半边的索引值由rightIdx维护，初始化为startIdx + N / 2，即当前分组右半边的第一个索引值；这两个索引值同步向前递增，从而访问到分组内的每一对系数。在循环体中，我们同时计算并更新一对系数，然后更新curOmega。
8.4 卷积的实现
卷积的实现是前面所有准备工作的汇总。让我们先梳理一下实现思路：
形参方面，卷积的输入是两个不保证等长的vector<double>

在进行快速傅里叶变换之前，我们需要先准备好足够多的点来表示结果多项式，即：需要将输入的两个系数列表都用0扩充到足够的长度。那么，需要多少个系数呢？这里需要做一个简单的计算：一个长度为\(n\的系数列表，表示的是一个\(n-1\次多项式；而另一个长度为\(m\的系数列表，表示的是一个\(m-1\次多项式；这两个多项式相乘的结果是一个\(n+m-2\次多项式；而这样的多项式一共有\(n+m-1\项。此外，根据快速傅里叶变换的要求，系数列表的长度必须是2的整数次幂，所以，我们还需要将\(n+m-1\这个数字对齐到2的整数次幂，作为两个系数列表扩充后的长度
快速傅里叶变换需要的系数列表是vector<complex<double>>类型的，而输入的系数列表是vector<double>类型的，需要进行转换
当两个系数列表都准备好后，进行两次快速傅里叶变换，将两个多项式从系数表达转为点值表达；然后，通过一个循环进行点值表达下的多项式乘法；最后，再进行一次快速傅里叶逆变换，将点值表达转为系数表达
快速傅里叶逆变换的输出是vector<complex<double>>类型的系数列表，而我们最终需要的是vector<double>类型的系数列表，需要进行转换。此外，还需要舍去由于对齐到2产生的系数扩充

vector<double> calcVectorConvolution(const vector<double> &leftCoefList, const vector<double> &rightCoefList
{
    unsigned coefSize  = leftCoefList.size( + rightCoefList.size( - 1;
    unsigned alignSize = __nextPow2(coefSize;

    vector<complex<double>> leftAlignCoefList(alignSize;
    vector<complex<double>> rightAlignCoefList(alignSize;

    for (unsigned idx = 0; idx < leftCoefList.size(; idx++
    {
        leftAlignCoefList[idx].real(leftCoefList[idx];
    }

    for (unsigned idx = 0; idx < rightCoefList.size(; idx++
    {
        rightAlignCoefList[idx].real(rightCoefList[idx];
    }

    auto leftFFTList  = __FFT(leftAlignCoefList, 1.;
    auto rightFFTList = __FFT(rightAlignCoefList, 1.;

    for (unsigned idx = 0; idx < leftFFTList.size(; idx++
    {
        leftFFTList[idx] *= rightFFTList[idx];
    }

    auto resFFTList  = __FFT(leftFFTList, -1.;

    vector<double> resCoefList(coefSize;

    for (unsigned idx = 0; idx < resCoefList.size(; idx++
    {
        resCoefList[idx] = resFFTList[idx].real(;
    }

    return resCoefList;
}

形参方面，leftCoefList和rightCoefList是两个多项式的系数表达。
coefSize用于存放结果多项式的项数，其计算公式已经由上文讨论过；alignSize用于存放将coefSize对齐到2的整数次幂后的多项式的项数，其决定了快速傅里叶变换需要的列表长度。
alignSize作为长度生成leftAlignCoefList和rightAlignCoefList，这两个系数列表用于快速傅里叶变换；并将leftCoefList和rightCoefList中的系数分别放入这两个列表的前面部分。
leftAlignCoefList和rightAlignCoefList从系数表达转为点值表达，存放在leftFFTList和rightFFTList中；然后，使用一个循环进行点值表达下的多项式乘法，将rightFFTList乘入leftFFTList中；最后，再进行一次快速傅里叶逆变换，将leftFFTList从点值表达转为系数表达，存放在resFFTList中。
resFFTList中存放的是alignSize个complex<double>类型的系数，而我们需要的是这个列表中前coefSize个double类型的系数，所以，函数的最后一段用于提取这部分系数并返回。读者可以自行验证：在resFFTList中，除了我们提取出的部分外，其余部分无论是实部还是虚部，都是0。
第9章 讨论
第三章中研究的欧拉公式的推导过程不能作为其证明过程。这是因为，麦克劳林级数展开需要求出函数的高阶导数，而复数域下的指数函数以及三角函数的导数均依赖于欧拉公式，从而造成循环论证。欧拉公式的严格证明超出了本文的范围。
已经存在不要求n为2的整数次幂的快速傅里叶变换算法，但其超出了本文的范围。
快速傅里叶变换在自然科学，计算机科学等诸多领域都有着广泛的应用，希望本文能够为读者提供帮助。
樱雨楼
2023.2
            声明：本文由用户投稿上传，本站不保证内容的真实与正确性，并且不承担相关法律责任，如有侵权请提供版权资料并联系删除！
            
                编程笔记
                »
                快速傅里叶变换及其实现
            
            
                
                    赞同 (42)
                
                or
                
                    
                    分享 (0)

快速傅里叶变换及其实现

第1章引言

第2章多项式的系数表达与点值表达

2.1 系数表达

2.2 点值表达

第3章欧拉公式

第4章 n次单位复数根

第5章离散傅里叶变换与离散傅里叶逆变换

第6章快速傅里叶变换

6.1 快速傅里叶变换的数学原理

6.2 一个手工计算快速傅里叶变换的实例

6.3 过程更简略的手工计算实例

第7章快速傅里叶逆变换

7.1 快速傅里叶逆变换的数学原理

7.2 一个手工计算快速傅里叶逆变换的实例

第8章快速傅里叶变换的实现

8.1 对齐至2的整数次幂的算法

8.2 位逆序算法

8.3 快速傅里叶变换与快速傅里叶逆变换的实现

8.4 卷积的实现

第9章讨论

相关文章

Hi，您需要填写昵称和邮箱！

第1章 引言

第2章 多项式的系数表达与点值表达

2.1 系数表达

2.2 点值表达

第3章 欧拉公式

第4章 n次单位复数根

第5章 离散傅里叶变换与离散傅里叶逆变换

第6章 快速傅里叶变换

6.1 快速傅里叶变换的数学原理

6.2 一个手工计算快速傅里叶变换的实例

6.3 过程更简略的手工计算实例

第7章 快速傅里叶逆变换

7.1 快速傅里叶逆变换的数学原理

7.2 一个手工计算快速傅里叶逆变换的实例

第8章 快速傅里叶变换的实现

8.1 对齐至2的整数次幂的算法

8.2 位逆序算法

8.3 快速傅里叶变换与快速傅里叶逆变换的实现

8.4 卷积的实现

第9章 讨论

相关文章

Hi，您需要填写昵称和邮箱！

第1章引言

第2章多项式的系数表达与点值表达

第3章欧拉公式

第5章离散傅里叶变换与离散傅里叶逆变换

第6章快速傅里叶变换

第7章快速傅里叶逆变换

第8章快速傅里叶变换的实现

第9章讨论