9.3 快速傅里叶变换

一、概述

线性代数的很多应用都很难在短时间内介绍清楚，一般都会在完善理论和新增的应用之间做权衡。通常是理论介绍比较多，但是这一节不同，本节介绍的是 $20$ 世纪最有价值的数值计算算法。
在用傅里叶矩阵 $F$ 和其逆矩阵 $F^{-1}$ 左乘其它矩阵或向量时，我们希望能更快的计算，快速傅里叶变换（FFT：Fast Fourier Transform） 可以实现这一点。通常情况下计算乘积 $F\mathbf{c}$ 需要 $n^2$ 次乘法，这是因为 $n$ 阶傅里叶矩阵 $F$ 有 $n^2$ 个元素，而 FFT 只需要 $\frac{1}{2}n{\log }_{2}n$ 次乘法。后面会详细介绍。
FFT 彻底改变了信号处理领域，整个产业的发展也被这一思想加速推进。电气工程师首先感受到了这种变换 —— 只要碰到函数，它们就会进行傅里叶变换。傅里叶变换是将 $f$ 表示为谐波函数 $c_k{e}^{ikx}$ 和的形式。函数 $f$ 通过系数 $c_k$ 被看成频率空间中的元素，而不是物理空间中的函数值 $f(x)$. 联系系数 $c$ 与函数 $f$ 之间的桥梁就是傅里叶变换，而 FFT 就是这种计算的快速算法。

二、单位根和傅里叶矩阵

二次方程有两个根（或是一个二重根），$n$ 次方程有 $n$ 个根（计算重数），这是代数的基本定理，为了使得该定理成立，我们允许有复根。这一节关注的是一个特殊的方程 $z^n=1$，它的解 $n$ 是 “$n$ 次单位根（$n$th roots of unity）”，在复平面上对应沿着单位圆均匀分布的 $n$ 个点。
Figure 9.4 显示了 $z^8=1$ 的八个根，它们的间隔是 $\frac{1}{8}(360°)=45°$. 第一个根位于 $45°$ 或 $\theta =\frac{2\pi }{8}$ 弧度处，它就是复数 $w=e^{i\theta }=e^{i2\pi /8}$. 为了强调这是八次方根，我们将其记为 $w_8$，也可以将其写成 $\cos \frac{2\pi }{8}$ 和 $\sin \frac{2\pi }{8}$ 组合的形式，但是我们不这样写。沿着圆周绕行一周，我们可以得到另外七个八次方根 $w^2,w^3,\cdots ,w^8$。这些 $w$ 的乘幂最好是用极形式表示，因为我们只会用到角度，这 $8$ 个角度的弧度表示是 $\frac{2\pi }{8},\frac{4\pi }{8},\cdots ,\frac{16\pi }{8}=2\pi$，度数表示就是 $45°,90°,135°,\cdots ,360°$.

$1$ 的四次方根也在图中，它们分别是 $i,-1,-i,1$，间隔的角度是 $\frac{2\pi }{4}$ 或 $90°$，第一个根 $w_4=e^{2\pi i/4}$ 就是 $i$；甚至 $1$ 的平方根也在图中，第一个根 $w_2=e^{i2\pi /2}=-1$，这里我们不要小瞧了平方根 $1$ 和 $-1$. FFT 背后的思想是将 $8\times 8$ 的傅里叶矩阵（包含八次方根 $w_8$ 的幂）变换成下面的 $4\times 4$ 傅里叶矩阵（包含四次方根 $w_4=i$ 的幂），同样的思想将四阶变换成二阶。通过探究 $F_8$ 降阶到 $F_4$ 以及升阶到 $F_{16}$（甚至更高阶）的规律，FFT 使得用 $F_{1024}$ 左乘非常迅速。
首先是阶数 $n=4$ 的傅里叶矩阵，它的行包含 $1,w,w^2$ 和 $w^3$ 以及它们的幂，它们是 $1$ 的四次方根，它们及其幂的排列很特别。$$\begin{array}{l}傅里叶矩阵\\ \text{Fourier matrix}\\ n=4,w=i\end{array}F=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& w& w^2& w^3\\ 1& w^2& w^4& w^6\\ 1& w^3& w^6& w^9\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& i& i^2& i^3\\ 1& i^2& i^4& i^6\\ 1& i^3& i^6& i^9\end{array}\right]$$这个矩阵是一个对称矩阵，$F=F^T$，它不是一个埃尔米特矩阵，因为它的对角线元素不全是实数。但是由于 $(\frac{1}{2}{F}^H)(\frac{1}{2}F)=I$，所以 $\frac{1}{2}F$ 是一个酉矩阵：

由 $F$ 的列向量可得 $F^{H}F=4I$，$F$ 的逆矩阵是 $\frac{1}{4}{F}^H$，即 $F^{-1}=\frac{1}{4}\overline{F}$.

逆矩阵中将 $w=i$ 变成 $\overline{w}=-i$，这就将 $F$ 变成了 $\overline{F}$，即可轻松求得逆矩阵 $F^{-1}$. 当 FFT 实现了用 $F$ 快速左乘时，也就同样快速的实现了用 $\overline{F}$ 和 $F^{-1}$ 的左乘。
$n$ 阶傅里叶矩阵的每列长度是 $\sqrt{n}$，所以酉矩阵是 $Q=\frac{F}{\sqrt{n}}$ 和 $Q^{-1}=\frac{\overline{F}}{\sqrt{n}}$，我们为了避免使用 $\sqrt{n}$，所以直接使用 $F$ 和 $F^{-1}=\frac{\overline{F}}{n}$. 重点是用 $F$ 左乘 $(c_0,c_1,c_2,c_3)$：$$\begin{array}{l}4点傅里叶级数\\ \text{4-point Fourier series}\end{array}\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\end{array}\right]=F\mathbf{c}=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& w& w^2& w^3\\ 1& w^2& w^4& w^6\\ 1& w^3& w^6& w^9\end{array}\right]\left[\begin{array}{c}{c}_0\\ c_1\\ c_2\\ c_3\end{array}\right](\mathrm{9.3.1})$$输入是四个复系数 $c_0,c_1,c_2,c_3$，输出是四个函数值 $y_0,y_1,y_2,y_3$. 第一个输出是 $y_0=c_0+c_1+c_2+c_3$ 是傅里叶级数 $\sum _{k=0}^3{c}_k{e}^{ikx}$ 在 $x=0$ 处的值，第二个输出是级数 $\sum _{k=0}^3{c}_k{e}^{ikx}$ 在 $x=\frac{2\pi }{4}$ 处的值：$$y_1=c_0+c_1{e}^{i2\pi /4}+c_2{e}^{i4\pi /4}+c_3{e}^{i6\pi /4}=c_0+c_{1}w+c_2{w}^2+c_3{w}^3$$第三个输出 $y_2$ 和第四个输出 $y_3$ 是 $\sum _{k=0}^3{c}_k{e}^{ikx}$ 分别在 $x=\frac{4\pi }{4}$ 和 $x=\frac{6\pi }{4}$ 处的值。这些是有限项傅里叶级数！它们包含 $n=4$ 项且为等距分布的 $n=4$ 个点，即点 $x=0,\frac{2\pi }{4},\frac{4\pi }{4},\frac{6\pi }{4}$ 处取值。

按照这个规律，下一个点是 $x=\frac{8\pi }{4}$，即是 $2\pi$，而由于 $e^{2\pi i}=e^0=1$，所以级数又回到了 $y_0$，每轮循环都以 $4$ 为周期，由于 $(w^2)(w^2)=w^0=1$，所以在这里指数 $2+2$ 与指数 $0$ 是相同的。依照惯例，规定行标 $j$ 和列标 $k$ 是从 $0$ 到 $n-1$（而不是 $1$ 到 $n$）. $F$ 的 “第零行” 和 “第零列” 都是全 $1$ 向量。
$n\times n$ 的傅里叶矩阵 $F_n$ 包含 $w=e^{2\pi i/n}$ 的幂：$$F_n\mathbf{c}=\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& w& w^2& \cdots & w^{n-1}\\ 1& w^2& w^4& \cdots & w^{2(n-1)}\\ ⋮& ⋮& ⋮& & ⋮\\ 1& w^{n-1}& w^{2(n-1)}& \cdots & w^{(n-1{)}^2}\end{array}\right]\left[\begin{array}{c}{c}_0\\ c_1\\ c_2\\ ⋮\\ c_{n-1}\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_{n-1}\end{array}\right]=y(\mathrm{9.3.2})$$$F_n$ 是对称矩阵，但不是埃尔米特矩阵。它的列是正交的，且 $F_n{\overline{F}}_n=nI$，其逆矩阵包含 ${\overline{w}}_n=e^{-2\pi i/n}$ 的幂。观察 $F$ 中元素的特征：

第 $j$ 行，第 $k$ 列的元素为 $w^{jk}$. 零行零列的元素均为 $w^0=1$.

当用 $F$ 左乘 $\mathbf{c}$ 时，相当于在 $n$ 个点处求级数的和。当用 $F_n^{-1}$ 左乘 $\mathbf{y}$ 时，即是由函数值 $\mathbf{y}$ 求系数 $\mathbf{c}$. 在 MATLAB 中，对应的命令是 $\mathbf{c}=\text{fft}(\mathbf{y})$. 矩阵 $F$ 能够将 “频率空间” 转换到 “物理空间” 即 “时域空间”。

重要注释： 很多作者更喜欢使用 $\omega =e^{-2\pi i/N}$，也就是这里 $w$ 的复共轭。（它们使用的是希腊字母 $\omega$，这里使用的是英文字母 $w$.）若使用这种选择，则 DFT 矩阵包含的是 $\omega$ 的幂而不是 $w$ 的幂，它就是 $\overline{F}$，$F$ 的共轭。$\overline{F}$ 将时域空间变换到频率空间。
$\overline{F}$ 是一个完全合理的选择！MATLAB 使用的是 $\omega =e^{-2\pi i/N}$，DFT 矩阵 fft(eye(N)) 包含数值 $w=\overline{w}$ 的幂。傅里叶矩阵 $F$ 的元素是 $w$ 的幂，它由 $\mathbf{c}$ 重构 $\mathbf{y}$，而 $\overline{F}$ 的元素是 $\omega$ 的幂，它用 fft(y) 来计算傅里叶系数。
另一个重要注释： 当函数 $f(x)$ 的周期是 $2\pi$ 时，我们将 $x$ 变为 $e^{i\theta }$，则这个函数就被定义在单位圆上（这里 $z=e^{i\theta }$）. 离散傅里叶变换（Discrete Fourier Transform）等同于插值（interpolation），即求出在 $n$ 个点 $z=1,w,\cdots ,w^{n-1}$ 处分别取 $n$ 个值 $f_0,f_1,\cdots ,f_{n-1}$ 的多项式 $p(z)=c_0+c_{1}z+\cdots +c_{n-1}{z}^{n-1}$：

插值 Interpolation$求出c_0,c_1,\cdots ,c_{n-1}，使得在n个点z=1,w,\cdots ,w^{n-1}处有p(w^j)=f_j$

傅里叶矩阵就是在 $n$ 个特殊点处插值时所得方程组的范德蒙德矩阵（Vandermonde matrix）.

三、快速傅里叶变换的一步

我们希望用 $F$ 左乘 $\mathbf{c}$ 时的计算越快越好。正常情况下，一个 $n$ 阶方阵有 $n^2$ 个元素，所以这个矩阵乘一个向量需要计算 $n^2$ 次乘法。如果矩阵中含有零元素，则对应的乘法可以被忽略，但是傅里叶矩阵里没有零元素！所以我们可能以为该算法不可能再得到改进了。而通过将其中的元素改写成特殊的形式 $w^{jk}$，将 $F$ 进行分解，从而产生很多零元素，这就是 FFT .
关键思想是将 $F_n$ 与阶数减半的傅里叶矩阵 $F_{n/2}$ 联系起来。假设 $n$ 为 $2$ 的幂（比如说 $n=2^{10}=1024$），下面我们将要把 $F_{1024}$ 与两个 $F_{512}$ 联系起来。
先从 $n=4$ 开始，关键是建立 $F_4$ 和两个 $F_2$ 之间的联系：$$F_4=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& i& i^2& i^3\\ 1& i^2& i^4& i^6\\ 1& i^3& i^6& i^9\end{array}\right],\left[\begin{array}{cc}{F}_2& \\ & F_2\end{array}\right]=\left[\begin{array}{cccc}1& 1& & \\ 1& i^2& & \\ & & 1& 1\\ & & 1& i^2\end{array}\right]$$左边的 $F_4$ 中没有零元素，右边的矩阵有一半的零元素，因此计算量会减半。但是，这两个矩阵还不相等，我们需要两个稀疏且简单的矩阵来完成 FFT 分解：$$\text{FFT}分解F_4=\left[\begin{array}{cccc}1& & 1& \\ & 1& & i\\ 1& & -1& \\ & 1& & -i\end{array}\right]\left[\begin{array}{cccc}1& 1& & \\ 1& i^2& & \\ & & 1& 1\\ & & 1& i^2\end{array}\right]\left[\begin{array}{cccc}1& & & \\ & & 1& \\ & 1& & \\ & & & 1\end{array}\right](\mathrm{9.3.3})$$最后一个矩阵是置换矩阵，它将偶数号系数（$c_0$ 和 $c_2$）放在了奇数号系数（$c_1$ 和 $c_3$）之前；中间的矩阵分别在偶数号和奇数号系数上执行规模减半的变换 $F_2$ 和 $F_2$；左边的矩阵将两个规模减半的输出合并起来，这样就可以得到正确的完整输出 $\mathbf{y}=F_4\mathbf{c}$.
同样的思想也可以应用在 $n=1024$ 和 $m=\frac{1}{2}n=512$ 的情形下。此时数 $w$ 是 $e^{2\pi i/1024}$，它在单位圆弧度 $\theta =\frac{2\pi }{1024}$ 处，傅里叶矩阵 $F_{1024}$ 中的元素全是 $w$ 的幂。FFT 的第一步就是这项伟大的分解，它是由 Cooley 和 Tukey 给出的（1805 年高斯曾经预言过）：

$$F_{1024}=\left[\begin{array}{cc}{I}_{512}& D_{512}\\ I_{512}& -D_{512}\end{array}\right]\left[\begin{array}{cc}{F}_{512}& \\ & F_{512}\end{array}\right]\left[\begin{array}{c}\text{even-odd}\\ \text{permutation}\end{array}\right](\mathrm{9.3.4})$$

$I_{512}$ 是单位矩阵，$D_{512}$ 是对角矩阵，其对角元素为 $1,w,\cdots ,w^{511}$；两个 $F_{512}$ 正是我们所期望的，注意这里使用的是 $512$ 次单位根（即 $w^2$ !!）；置换矩阵将输入向量 $\mathbf{c}$ 分成偶数号部分 ${\mathbf{c}}^{\prime }=(c_0,c_2,\cdots ,c_{1022})$和奇数号部分 ${\mathbf{c}}^{\prime \prime }=(c_1,c_3,\cdots ,c_{1023})$.
下面是代数公式，它与 $F_{1024}$ 的分解描述是一样的：

（FFT 的一步）令 $m=\frac{1}{2}n$，$\mathbf{y}=F_n\mathbf{c}$ 的前 $m$ 位分量和后 $m$ 位分量分别为规模减半的变换 ${\mathbf{y}}^{\prime }=F_m{\mathbf{c}}^{\prime }$ 和 ${\mathbf{y}}^{\prime \prime }=F_m{\mathbf{c}}^{\prime \prime }$. 这对应式（9.3.4）中的 $I{\mathbf{y}}^{\prime }+D{\mathbf{y}}^{\prime \prime }$ 和 $I{\mathbf{y}}^{\prime }-D{\mathbf{y}}^{\prime \prime }$，这个等式展示了从 $n$ 到 $m=\frac{n}{2}$ 的步骤：$$\begin{array}{r}{y}_j=y_j^{\prime }+(w_n{)}^j{y}_j^{\prime \prime },j=0,1,\cdots ,m-1\\ y_{j+m}=y_j^{\prime }-(w_n{)}^j{y}_j^{\prime \prime },j=0,1,\cdots ,m-1\end{array}(\mathrm{9.3.5})$$将 $\mathbf{c}$ 分为 ${\mathbf{c}}^{\prime }$ 和 ${\mathbf{c}}^{\prime \prime }$，用 $F_m$ 分别将它们变换为 ${\mathbf{y}}^{\prime }$ 和 ${\mathbf{y}}^{\prime \prime }$，从而（9.3.5）重构了 $\mathbf{y}$.

这些公式来自于将 $c_0,c_1,\cdots ,c_{n-1}$ 分成偶数号 $c_{2k}$ 和奇数号 $c_{2k+1}$，这里的 $w$ 是 $w_n$：$$\mathbf{y}=F\mathbf{c}{y}_j=\sum _{k=0}^{n-1}{w}^{jk}{c}_k=\sum _{k=0}^{m-1}{w}^{2jk}{c}_{2k}+\sum _{k=0}^{m-1}{w}^{j(2k+1)}{c}_{2k+1},其中m=\frac{1}{2}n(\mathrm{9.3.6})$$上式即（9.3.2）的展开形式，偶数号分量是 ${\mathbf{c}}^{\prime }=(c_0,c_2,\cdots ,c_{n-2})$，奇数号分量是 ${\mathbf{c}}^{\prime \prime }=(c_1,c_3,\cdots ,c_{n-1})$，然后进行变换 $F_m{\mathbf{c}}^{\prime }$ 和 $F_m{\mathbf{c}}^{\prime \prime }$，关键是 $w_n^2=w_m$，这得到 $w_n^{2jk}=w_m^{jk}$.$$重写(\mathrm{9.2.6})y_j=\sum _{k=0}^{m-1}(w_m{)}^{jk}{c}_k^{\prime }+(w_n{)}^j\sum _{k=0}^{m-1}(w_m{)}^{jk}{c}_k^{\prime \prime }=y_j^{\prime }+(w_n{)}^j{y}_j^{\prime \prime }(\mathrm{9.3.7})$$当 $j\ge m$ 时，（9.3.5）中的负号都来自于 $(w_n{)}^j$ 所分解出的 $(w_n{)}^m=-1$.
MATLAB 中可以使用 conj(F) 或傅里叶逆变换 ifft 轻易的将偶数号与奇数号分量分开，再乘上 $w_n^j$，因为 fft 是基于 $\omega =\overline{w}=e^{-2\pi i/n}$. $F$ 和 conj(F) 可以通过行置换联系起来。$$\begin{array}{l}\text{MATLAB}中从\\ n到n/2的\text{FFT}\\ 步骤\end{array}\begin{array}{l}{y}^{\prime }=\text{ifft}(c(0:2:n-2))\ast n/2;\\ y^{\prime \prime }=\text{ifft}(c(1:2:n-1))\ast n/2;\\ d=w.\text{^}(0:n/2-1{)}^{\prime };\\ y=[y^{\prime }+d.\ast y^{\prime \prime };y^{\prime }-d.\ast y^{\prime \prime }];\end{array}$$上述代码需要注意的地方：实际上 MATLAB 的标号是从 1 开始的，上述仅仅是伪代码；之所以后面会乘上 $n/2$，这是因为 MATLAB 中的 ifft 会自动归一化，即它会除以 $n/2$，这里再乘回来，不做归一化。
下面的流程图表明了用规模减半的矩阵 $F_2$ 左乘 ${\mathbf{c}}^{\prime }$ 和 ${\mathbf{c}}^{\prime \prime }$. 根据这个形状，我们称这个步骤称为 “蝴蝶形 butterflies”，然后将输出 ${\mathbf{y}}^{\prime }$ 和 ${\mathbf{y}}^{\prime \prime }$ 结合起来（用 $D$ 中的 $1,i$ 乘 ${\mathbf{y}}^{\prime \prime }$，再用 $-D$ 中的 $-1,-i$ 乘 ${\mathbf{y}}^{\prime \prime }$），从而得到 $\mathbf{y}=F_4\mathbf{c}$.

将 $F_n$ 化成两个 $F_m$ 几乎将计算量减半 —— 从分解的矩阵中出现的很多零元素可以看出。这种减半虽然很好但是不够彻底，FFT 的完整思想要更明显，它要节省出远超过一半的时间。

四、利用递归的完整版 FFT

到目前为止，我们只是将 $F_n$ 减少到了 $F_{n/2}$，我们还要继续将其减少至 $F_{n/4}$. 将每个 $F_{512}$ 都化成两个 $F_{256}$，然后将每个 $F_{256}$ 化成两个 $F_{128}$. 这就是递归（recursion）.
递归是很多快速算法的底层原理。我们下面进行第 $2$ 步，此时出现了四个 $F_{256}$ 和 $D$（包含 ${\omega }_{512}$ 的 $256$ 个幂）. 偶数号中的偶数号 $c_0,c_4,c_8,\cdots$ 会首先出现：$$\left[\begin{array}{cc}{F}_{512}& \\ & F_{512}\end{array}\right]=\left[\begin{array}{cccc}I& D& & \\ I& -D& & \\ & & I& D\\ & & I& -D\end{array}\right]\left[\begin{array}{cccc}F& & & \\ & F& & \\ & & F& \\ & & & F\end{array}\right]\left[\begin{array}{cc}\text{pick}& 0,4,8,\cdots \\ \text{pick}& 2,6,10,\cdots \\ \text{pick}& 1,5,9,\cdots \\ \text{pick}& 3,7,11,\cdots \end{array}\right]$$下面我们来计算单个乘法的次数，以搞清节省了多少计算量。在 FFT 发现之前，一般情况下我们需要 $n^2=(1024{)}^2$，这是有超过一百万次的乘法，可能这样计算并不会花费太长的时间，但是如果需要做很多次变换时（这是比较普遍的），时间的消耗就会非常巨大。而通过 FFT 我们可以节省相当可观的计算量：$$对于阶数n=2^l的矩阵，最终的乘法次数由n^2减少到\frac{1}{2}nl$$$1024$ 就是 $2^{10}$，所以 $l=10$. 那么乘法次数就由最初的 $(1024{)}^2$ 减少到了 $(5)(1024)$ ，即从一百万次减少到了五千多次，减少为原来的 $\frac{1}{200}$ 了。这也就是为什么说 FFT 在彻底改变了信号处理了。
下面解释乘法次数为什么 $\frac{1}{2}nl$。由 $n=2^l$ 降到 $n=1$ 的过程中，总共有 $l$ 个层级。每一级为了重新组合从更低层级而来的规模减半的输出，关于对角矩阵 $D$ 有 $\frac{2}{n}$ 次乘法，所以可以由此得到最终的 $\frac{1}{2}nl$ 次乘法，这个就是 $\frac{1}{2}n{\log }_{2}n$.
关于 FFT 这个意义非凡的算法最后一条注记，在对向量分量进行全部偶-奇置换后，这些分量进行 FFT 时的次序有一条令人惊叹的规则。将数字 $0$ 到 $n-1$ 写成二进制形式（如 $n=4$ 时的 $00,01,10,11$），然后反转每个数字的次序得到：$00,10,01,11$，这样得到了位反转排序（bit-reversed order）0,2,1,3，其中偶数排在奇数之前。向量分量按位反转排序，共有 $l={\log }_{2}n$ 次递归步骤，最终的输出 $y_0,y_1,\cdots ,y_{n-1}$ 就是 $F_n$ 左乘 $\mathbf{c}$.
这一节展示了矩阵左乘向量这一基本运算，也有可能对其计算进行优化。