文章目录

• 3. DFT与FFT
• • 3.1 单位复数根
  • 3.2 DFT
  • 3.3 FFT
  • 3.4 在单位复数根处插值
• 4. 高效FFT实现
• • 4.1 FFT的一种迭代实现
  • 4.2 并行FFT电路

3. DFT与FFT

（算导30.1节中断言）如果使用单位复数根，可以在 $\Theta (n\log n)$ 时间内完成求值与插值运算。在本节中给出单位复数根的定义，研究其性质，以及定义DFT，然后说明FFT如何仅用 $\Theta (n\log n)$ 时间、就可以计算出DFT和它的逆。

3.1 单位复数根

$n$ 次单位复数根 complex nth root of unity 是满足 ${\omega }^n=1$ 的复数 $\omega$（类比平方根、立方根，且是单位根、复数）。$n$ 次单位复数根恰好有 $n$ 个：对于 $k=0,1,\dots ,n-1$ 这些根是 $e^{2\pi ik/n}(i=\sqrt{-1})$ 。为了（理解）解释这个表达式，我们利用复数指数的定义 the definition of the exponential of a complex number ：$$e^{iu}=\cos (u)+i\sin (u)$$ 图30-2说明，$n$ 个单位复数根均匀地分布在「以复平面的原点为圆心的、单位半径的圆周」上。值 $${\omega }_n=e^{2\pi i/n}\text{\hspace{1em}(30.6)}$$ 称为主 $n$ 次单位根 the principal nth root of unity ，其余的 $n$ 次单位复数根都是 ${\omega }_n={\omega }_n^1$ 的幂次，即 $({\omega }_n^1{)}^k={\omega }_n^k$ 。注意，有时别人对 ${\omega }_n$ 有不同的定义：${\omega }_n=e^{-2\pi i/n}$ ，这个可替的定义一般用在信号处理应用中，这两个 ${\omega }_n$ 的定义、其背后的数学含义基本上是相同的。
例如（红字为主 $n$ 次单位根），$$\begin{array}{rlr}& {\omega }_8^0=e^{\frac{2\pi \times 0}{8}i}=\cos (0)+i\sin (0)& 复平面(1,0)\\ & {\omega }_8^1=e^{\frac{2\pi \times 1}{8}i}=\cos (\frac{\pi }{4})+i\sin (\frac{\pi }{4})& 复平面(\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2})\\ & {\omega }_8^2=e^{\frac{2\pi \times 2}{8}i}=\cos (\frac{\pi }{2})+i\sin (\frac{\pi }{2})& 复平面(0,1)\\ & {\omega }_8^3=e^{\frac{2\pi \times 3}{8}i}=\cos (\frac{3\pi }{4})+i\sin (\frac{3\pi }{4})& 复平面(-\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2})\\ & {\omega }_8^4=e^{\frac{2\pi \times 4}{8}i}=\cos (\pi )+i\sin (\pi )& 复平面(-1,0)\\ & {\omega }_8^5=e^{\frac{2\pi \times 5}{8}i}=\cos (\frac{5\pi }{4})+i\sin (\frac{5\pi }{4})& 复平面(-\frac{\sqrt{2}}{2},-\frac{\sqrt{2}}{2})\\ & {\omega }_8^6=e^{\frac{2\pi \times 6}{8}i}=\cos (\frac{3\pi }{2})+i\sin (\frac{3\pi }{2})& 复平面(0,-1)\\ & {\omega }_8^7=e^{\frac{2\pi \times 7}{8}i}=\cos (\frac{7\pi }{4})+i\sin (\frac{7\pi }{4})& 复平面(\frac{\sqrt{2}}{2},-\frac{\sqrt{2}}{2})\\ & {\omega }_8^8=e^{\frac{2\pi \times 8}{8}i}=\cos (2\pi )+i\sin (2\pi )& 复平面(1,0)\end{array}$$

$n$ 个 $n$ 次单位复数根 ${\omega }_n^0,{\omega }_n^1,\dots ,{\omega }_n^{n-1}$ 在乘法意义下形成一个群（见算导31.3节）。该群与加法群 $({\mathbb{Z}}_n,+)$ 具有相同的结构，因为 ${\omega }_n^n={\omega }_n^0=1$ 意味着 ${\omega }_n^j{\omega }_n^k={\omega }_n^{(j+k)\mathrm{m}\mathrm{o}\mathrm{d}n}$ 。类似地，${\omega }_n^{-1}={\omega }_n^{n-1}$ 。下面的引理给出了 $n$ 次单位复数根的一些基本性质。

引理30.3（消去引理 Cancellation lemma ）对任何整数 $n\ge 0,k\ge 0$ ，以及 $d>0$（$n$ 不应该 $>0$ 吗？），$${\omega }_{dn}^{dk}={\omega }_n^k\text{\hspace{1em}(30.7)}$$ 证明：由式30.6可以直接推出引理，因为：$${\omega }_{dn}^{dk}=(e^{\frac{2\pi i}{dn}}{)}^{dk}=(e^{\frac{2\pi i}{n}}{)}^k={\omega }_n^k\text{\hspace{1em}■}$$
推论30.4 对任意偶数 $n>0$ ，有 $${\omega }_n^{n/2}={\omega }_2=-1$$ 证明：证明见算导练习30.2-1。本人的证明如下：由于 $n$ 为大于零的偶数，所以 $n/2\times 2$ 仍等于 $n$ ，于是根据引理30.3有：$${\omega }_n^{n/2}={\omega }_{2n}^n=(e^{\frac{2\pi i}{2n}}{)}^n=e^{\frac{2\pi i}{2}}={\omega }_2$$ 又可知：$${\omega }_2=e^{\frac{2\pi i}{2}}=\cos (\pi )+i\sin (\pi )=-1$$ 所以推论30.4可证。$■$

引理30.5（折半引理 Halving lemma ）如果 $n>0$ 为偶数，那么 $n$ 个 $n$ 次单位复数根的平方的集合，就是 $n/2$ 个 $n/2$ 次单位复数根的集合 If n > 0 is even, then the squares of the n complex nth roots of unity are the n/2 complex (n/2)/th roots of unity 。
证明：根据消去引理，对任意非负整数 $k$ ，我们有 $({\omega }_n^k{)}^2={\omega }_n^{2k}={\omega }_{n/2}^k$（即对一个 $n$ 次单位复数根平方后得到一个 $n/2$ 次单位复数根）。注意，如果对所有 $n$ 次单位复数根进行平方，那么获得每个 $n/2$ 次单位根正好 $2$ 次，因为： $$({\omega }_n^{k+n/2}{)}^2={\omega }_n^{2k+n}={\omega }_n^{2k}{\omega }_n^n={\omega }_n^{2k}=({\omega }_n^k{)}^2$$ 因此可知，${\omega }_n^k$ 与 ${\omega }_n^{k+n/2}$ 平方相同。

我们也可以由推论30.4来证明该性质，因为 ${\omega }_n^{n/2}=-1$ 意味着 ${\omega }_n^{k+n/2}=-{\omega }_n^k$ ，从而 $({\omega }_n^{k+n/2}{)}^2=({\omega }_n^k{)}^2$ 。$■$

我们将会看到，对于用分治策略来对「多项式的系数与点值表达」进行相互转换，折半引理是非常重要的，因为它保证递归子问题的规模只是递归调用前的一半。

引理30.6（求和引理 Summation lemma ）对任意整数 $n\ge 1$ 和不能被 $n$ 整除的非负整数 $k$ ，有 $$\sum _{j=0}^{n-1}({\omega }_n^k{)}^j=0$$ 证明：对于实数 $x\ne 1$ ，和式：$$\sum _{k=0}^n{x}^k=1+x+x^2+\cdots +x^n=\frac{x^{n+1}-1}{x-1}$$ 这是几何级数/等比级数/指数级数求和公式（算导A.5），既适用于实数，也适用于复数。因此有：
$$\sum _{j=0}^{n-1}({\omega }_n^k{)}^j=\frac{({\omega }_n^k{)}^n-1}{{\omega }_n^k-1}=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}=\frac{(1{)}^k-1}{{\omega }_n^k-1}=0$$ 因为要求 $k$ 不能被 $n$ 整除、而仅当 $k$ 被 $n$ 整除时 ${\omega }_n^k=1$ 成立，我们确保分母不为 $0$ 。$■$

3.2 DFT

回顾一下，我们希望计算次数界为 $n$ 的多项式 $$A(x)=\sum _{j=0}^{n-1}{a}_j{x}^j$$ 在 ${\omega }_n^0,{\omega }_n^1,{\omega }_n^2,\dots ,{\omega }_n^{n-1}$ 处的值（即在 $n$ 个 $n$ 次复数根处）。假设 $A$ 以系数形式给出，即 $a=(a_0,a_1,\dots ,a_{n-1})$ 。接下来对 $k=0,1,\dots ,n-1$ ，定义结果 $y_k$ ：$$y_k=A({\omega }_n^k)=\sum _{j=0}^{n-1}{a}_j{\omega }_n^{kj}\text{\hspace{1em}(30.8)}$$ 向量 $y=(y_0,y_1,\dots ,y_{n-1})$ 就是系数向量 $a=(a_0,a_1,\dots ,a_{n-1})$ 的离散傅里叶变换 DFT 。我们也记为 $y=\mathrm{D}\mathrm{F}{\mathrm{T}}_{\mathrm{n}}(a)$ 。

3.3 FFT

通过使用一种称为快速傅里叶变换 FFT 的方法，利用复数单位根的特殊性质，我们就能在 $\Theta (n\log n)$ 时间内计算出 $\mathrm{D}\mathrm{F}{\mathrm{T}}_{\mathrm{n}}(a)$（即完成系数表达到点值表达的转换），而直接的方法所需时间为 $\Theta (n^2)$ 。这里全篇假设 $n$ 恰好是 $2$ 的整数幂，虽然处理非 $2$ 整数幂的策略已经存在，但算导正文没有提。

FFT利用了分治策略，采用 $A(x)$ 中偶数下标的系数与奇数下标的系数，分别定义两个新的、次数界为 $n/2$ 的多项式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ ： $$\begin{array}{rl}& A^{[0]}(x)=a_0+a_{2}x+a_4{x}^2+\cdots +a_{n-2}{x}^{n/2-1}\\ & A^{[1]}(x)=a_1+a_{3}x+a_5{x}^2+\cdots +a_{n-1}{x}^{n/2-1}\end{array}$$ 注意到，$A^{[0]}(x)$ 包含 $A$ 中所有偶数下标的系数（下标的相应二进制表示的最后一位为 $0$ ），$A^{[1]}(x)$ 包含 $A$ 中所有奇数下标的系数（下标的相应二进制表示的最后一位为 $1$ ）。于是有：$$A(x)=A^{[0]}(x^2)+x{A}^{[1]}(x^2)\text{\hspace{1em}(30.9)}$$

所以，求 $A(x)$ 在 ${\omega }_n^0,{\omega }_n^1,\dots ,{\omega }_n^{n-1}$ 处的值的问题，转换为：

1. 求次数界为 $n/2$ 的多项式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ 在点：$$({\omega }_n^0{)}^2,({\omega }_n^1{)}^2,\dots ,({\omega }_n^{n-1}{)}^2\text{\hspace{1em}(30.10)}$$ 的取值。
2. 根据式30.9综合上述结果。

根据折半引理（作用就在这里！），式30.10并不是由 $n$ 个不同值组成，而是仅由 $n/2$ 个 $n/2$ 次单位复数根组成，每个根正好出现 $2$ 次（$n$ 为 $2$ 的整数幂，也是偶数）。因此，我们递归地对次数界为 $n/2$ 的多项式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ 、在 $n/2$ 个 $n/2$ 次单位复数根处进行求值。这些子问题与原始问题形式相同、但规模变为一半。

现在，我们已经成功地把一个 $n$ 个元素的 $\mathrm{D}\mathrm{F}{\mathrm{T}}_{\mathrm{n}}$ 计算，划分为两个规模为 $n/2$ 个元素的 $\mathrm{D}\mathrm{F}{\mathrm{T}}_{\mathrm{n}/2}$ 计算。这一分解是下面递归FFT算法的基础，该算法计算出一个 $n$ 元素向量 $a=(a_0,a_1,\dots ,a_{n-1})$ 的DFT（其中 $n$ 是 $2$ 的幂）。 注意，这里的长度 $n$ 实际上是（算导30.1节）前面所指的 $2n$ ，因为我们需要在求值以前、加倍给定多项式的次数界，因此在多项式乘法的相关内容中，实际上处理的是 $2n$ 次单位根。

RECURSIVE-FFT(a)
	n = a.length // n is a power of 2
	if n == 1
		return a
	wn = e^(2πi/n)
	w = 1
	A0 = (a[0], a[2], ..., a[n-2]) // a中偶数下标的系数组成的多项式
	A1 = (a[1], a[3], ..., a[n-1]) // a中奇数下标的系数组成的多项式
	y0 = RECURSIVE-FFT(A0) // 递归FFT处理A0后得到的点值表达(偶数下标系数)
	y1 = RECURSIVE-FFT(A1) // 递归FFT处理A1后得到的点值表达(奇数下标系数)
	for k = 0 to (n/2-1)
		y[k] = y0[k] + w * y1[k]
		y[k + (n/2)] = y0[k] - w * y1[k]
		w = w * wn
	return y // y is assumed to be a column vector

RECURSIVE-FFT 的执行过程如下。前三行代表递归的基础：一个元素的DFT就是该元素自身，因为在这种情形下：$$y_0=a_0{\omega }_1^0=a_0\cdot 1=a_0$$ 定义 Y0, Y1 的两行，定义了多项式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ 的系数向量。定义 wn, w 的两行和 w = w * wn 这一行，则保证了 $\omega$ 可以正确更新，只要 y[k], y[k + (n/2)] 这两行被执行，就有 $\omega ={\omega }_n^k$（次次迭代中让 $\omega$ 的值改变，可以节约每次通过for循环重新计算 ${\omega }_n^k$ 的时间）。递归调用 RECURSIVE-FFT 的两行，执行递归计算 $\mathrm{D}\mathrm{F}{\mathrm{T}}_{\mathrm{n}/2}$ ，对于 $k=0,1,\dots ,n/2-1$ ：$$\begin{array}{rl}& y_k^{[0]}=A^{[0]}({\omega }_{n/2}^k)\\ & y_k^{[1]}=A^{[1]}({\omega }_{n/2}^k)\end{array}$$ 或者根据消去引理，有 ${\omega }_{n/2}^k={\omega }_n^{2k}=({\omega }_n^k{)}^2$ ，于是 $$\begin{array}{rl}& y_k^{[0]}=A^{[0]}({\omega }_n^{2k})\\ & y_k^{[1]}=A^{[1]}({\omega }_n^{2k})\end{array}$$ 更新 y[k], y[k + (n/2)] 的两行则综合了递归 $\mathrm{D}\mathrm{F}{\mathrm{T}}_{\mathrm{n}/2}$ 的计算结果。对 $y_0,y_1,\dots ,y_{n/2-1}$ ，前一行推出：$$\begin{array}{rlr}{y}_k& =y_k^{[0]}+{\omega }_n^k{y}_k^{[1]}& \\ & =A^{[0]}({\omega }_n^{2k})+{\omega }_n^k{A}^{[1]}({\omega }_n^{2k})& \\ & =A({\omega }_n^k)& (根据式30.9)\end{array}$$ 对 $y_{n/2},y_{n/2+1},\dots ,y_{n-1}$ ，设 $k=0,1,\dots ,n/2-1$ ，后一行推出：
$$\begin{array}{rlr}{y}_{k+(n/2)}& =y_k^{[0]}-{\omega }_n^k{y}_k^{[1]}& \\ & =y_k^{[0]}+{\omega }_k^{k+(n/2)}{y}_k^{[1]}& (因为{\omega }_n^{k+(n/2)}=-{\omega }_n^k)\\ & =A^{[0]}({\omega }_n^{2k})+{\omega }_n^{k+(n/2)}{A}^{[1]}({\omega }_n^{2k})& \\ & =A^{[0]}({\omega }_n^{2k+n})+{\omega }_n^{k+(n/2)}{A}^{[1]}({\omega }_n^{2k+n})& (因为{\omega }_n^{2k+n}={\omega }_n^{2k})\\ & =A({\omega }_n^{k+(n/2)})& (根据式30.9)\end{array}$$ 因此，由 RECURSIVE-FFT 返回的向量 y 确实是输入向量 a 的DFT。

更新 y[k], y[k + (n/2)] 的两行中，对 $k=0,1,\dots ,n/2-1$ ，每个值 $y_k^{[1]}$ 都乘了 ${\omega }_n^k$ 。前一行中这个乘积加到了 $y_k^{[0]}$ 上，后一行又减去它。因为应用了每个因子 ${\omega }_n^k$ 的正数形式和负数形式，我们把因子 ${\omega }_n^k$ 称为旋转因子 twiddle factor 。

为了确定过程 RECURSIVE-FFT 的运行时间，注意到除了递归调用外，每次调用所需的时间为 $\Theta (n)$ ，其中 $n$ 是输入向量的长度。因此，对运行时间有下列递归式：$$T(n)=2T(n/2)+\Theta (n)=\Theta (n\log n)$$ 因此，采用快速傅里叶变换，我们可以在 $\Theta (n\log n)$ 时间内，求出次数界为 $n$ 的多项式在 $n$ 次单位复数根处的值。

3.4 在单位复数根处插值

现在我们展示，如何在单位复数根处插值、来完成多项式乘法方案，使得我们把一个多项式从点值表达转换回系数表达。如下方法进行插值：把DFT写成一个矩阵方程，然后再观察其逆矩阵的形式。

根据等式30.4，我们可以把DFT写成矩阵乘积 $y=V_{n}a$ ，其中 $V_n$ 是一个由 ${\omega }_n$ 适当幂次填充成的范德蒙矩阵：
$$\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ ⋮\\ y_{n-1}\end{array}\right]=\left[\begin{array}{cccccc}1& 1& 1& 1& \dots & 1\\ 1& {\omega }_n^1& {\omega }_n^2& {\omega }_n^3& \dots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& {\omega }_n^6& \dots & {\omega }_n^{2(n-1)}\\ 1& {\omega }_n^3& {\omega }_n^6& {\omega }_n^9& \dots & {\omega }_n^{3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& {\omega }_n^{3(n-1)}& \dots & {\omega }_n^{(n-1)(n-1)}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ ⋮\\ a_{n-1}\end{array}\right]$$ 对 $j,k=0,1,\dots ,n-1$ ，$V_n$ 的 $(k,j)$ 处元素为 ${\omega }_n^{kj}$ 。$V_n$ 中元素的指数组成一张乘法表。对于逆运算 $a=\mathrm{D}\mathrm{F}{\mathrm{T}}_{\mathrm{n}}^{-1}(y)$ ，我们把 $y$ 乘以 $V_n$ 的逆矩阵 $V_n^{-1}$ 来进行处理。

定理30.7 对 $j,k=0,1,\dots ,n-1$ ，$V_n^{-1}$ 的 $(j,k)$ 处元素为 ${\omega }_n^{-kj}/n$ 。
证明：我们要做的只是证明，这种形式下 $V_n^{-1}{V}_n=I_n$ ，其中 $I_n$ 为 $n\times n$ 的单位矩阵。考虑 $V_n^{-1}{V}_n$ 中 $(j,j^{\prime })$ 处的元素：
$$[V_n^{-1}{V}_n{]}_{j{j}^{\prime }}=\sum _{k=0}^{n-1}({\omega }_n^{-jk}/n)({\omega }_n^{k{j}^{\prime }})=\sum _{k=0}^{n-1}{\omega }_n^{k(j^{\prime }-j)}/n$$ 如果 $j^{\prime }=j$ ，则此和为 $1$ ；否则，根据求和引理（引理30.6），此和为 $0$ 。注意，只有 $-(n-1)\le j^{\prime }-j\le n-1$ ，使得 $j^{\prime }-j$ 不能被 $n$ 整除，才能应用求和引理。$■$

给定逆矩阵 $V_n^{-1}$ ，可以推导出 $\mathrm{D}\mathrm{F}{\mathrm{T}}_{\mathrm{n}}^{-1}(y)$ ：$$a_j=\frac{1}{n}\sum _{k=0}^{n-1}{\omega }_n^{-jk}{y}_k(j=0,1,\dots ,n-1)\text{\hspace{1em}(30.11)}$$ 通过比较式30.8与式30.11，我们可以看到，对FFT算法进行如下修改就可以计算出逆DFT（算导练习30.2-4）：把 $a$ 与 $y$ 互换，用 ${\omega }_n^{-1}$ 替换 ${\omega }_n$ ，并将计算结果的每个元素除以 $n$ 。因此，我们也可以在 $\Theta (n\log n)$ 时间内计算出 $\mathrm{D}\mathrm{F}{\mathrm{T}}_{\mathrm{n}}^{-1}$ 。

我们可以看到，通过运用FFT与逆FFT，可以在 $\Theta (n\log n)$ 时间内、把次数界为 $n$ 的多项式在其系数表达与点值表达之间进行相互转换。在矩阵乘法的相关内容中，已经说明了下面的结论。

定理30.8（卷积定理 Convolution theorem）对任意两个长度为 $n$ 的向量 $a,b$ ，其中 $n$ 是 $2$ 的整数幂。$$a\otimes b=\mathrm{D}\mathrm{F}{\mathrm{T}}_{2\mathrm{n}}^{-1}(\mathrm{D}\mathrm{F}{\mathrm{T}}_{2\mathrm{n}}(\mathrm{a})\cdot \mathrm{D}\mathrm{F}{\mathrm{T}}_{2\mathrm{n}}(\mathrm{b}))$$ 其中，向量 $a,b$ 用 $0$ 填充、使其长度达到 $2n$ ，并用 $\cdot$ 表示两个 $2n$ 元向量的点乘。$■$

4. 高效FFT实现

因为DFT的实际应用（如信号处理）中需要尽可能快的速度，本节将探究两种高效的FFT实现方法。首先，讨论一种运行时间为 $\Theta (n\log n)$ 的FFT迭代实现方法，不过在此运行时间的 $\Theta$ 记号，隐含的常数要比（算导30.2节）前面递归实现方法的常数小（如果实现精确，这个递归方法可能会更加高效地应用硬件缓存）。然后，深入分析迭代实现方法，设计出一个高效的并行FFT电路。

4.1 FFT的一种迭代实现

首先我们注意到，在 RECURSIVE-FFT 中，第10~13行的for循环中，包含了 ${\omega }_n^k{y}_k^{[1]}$ 的 $2$ 次计算。在编译术语中，称该值为公用子表达式 common subexpression 。我们可以改变循环、使其仅计算一次，并将其存放在临时变量 $t$ 中。

for k = 0 to n/2-1
	t = w * Y1[k]
	y[k] = Y0[k] + t
	y[k + (n/2)] = Y0[k] - t
	w = w * wn

在这个循环中，把旋转因子 $\omega ={\omega }_n^k$ 乘以 $y_k^{[1]}$ ，把所得乘积存入 $t$ 中，然后从 $y_k^{[0]}$ 中增加及减去 $t$ ，这一系列操作称为一个蝴蝶操作 butterfly operation ，图30-3说明了执行步骤：

现在来说明，如何使FFT算法采用迭代结构、而非递归结构。在图30-4中，我们已把输入向量分配在「一次 RECURSIVE-FFT 调用中的各次递归调用中」 、形成一个树形结构 we have arranged the input vectors to the recursive calls in an invocation of RECURSIVE-FFT in a tree structure ，其中初始调用时有 $n=8$ 。树中的每个结点对应每次递归过程调用，由相应的输入向量标记。每次 RECURSIVE-FFT 调用产生两个递归调用，除非该调用已收到了 $1$ 个元素的向量。第一次调用作为左孩子，第二次调用作为右孩子。

观察此树，我们注意到，如果把初始向量 $a$ 中的元素按其在叶中出现次序进行安排，就可以对过程 RECURSIVE-FFT 的执行进行追踪，不过是自底向上、而非自顶向下。

• 首先，我们成对取出元素，利用一次蝴蝶操作计算出每对的DFT，然后用其DFT取代这对元素，这样向量中就包含了 $n/2$ 个二元素的DFT。
• 下一步，我们按对取出这 $n/2$ 个DFT，通过两次蝴蝶操作计算出含有四个元素向量的DFT，并用一个具有四个元素的DFT取代对应的两个二元素的DFT。于是，向量中包含 $n/4$ 个四元素的DFT。
• 继续进行这一过程，直到向量包含两个具有 $n/2$ 个元素的DFT，这时我们综合应用 $n/2$ 次蝴蝶操作，就可以合成最终的具有 $n$ 个元素的DFT。

为了把这个自底向上的方法变为代码，采用了一个数组 $A[0\dots n-1]$ ，初始时该数组包含输入向量 $a$ 中的元素，其顺序为它们在图30-4中树叶出现的顺序（在后面说明如何确定这个顺序，这也称为位逆序置换 bit-reversal permutation ）。因为需要在树的每一层进行组合，于是引入一个变量 $s$ 以计算树的层次，取值范围为从 $1$（在最底层，这时我们组合对来构成二元素的DFT）到 $\log n$（在最顶层，这里我们要对两个具有 $n/2$ 个元素的DFT进行组合，以产生最后结果）。因此，这个算法有如下结构：

我们可以用更精确的伪代码描述第3行中的循环主体部分。从子程序 RECURSIVE-FFT 中复制for循环，让 $y^{[0]}$ 与 $A[k\dots k+2^{s-1}-1]$ 一致、$y^{[1]}$ 与 $A[k+2^{s-1}\dots k+2^s-1]$ 一致。在每次蝴蝶操作中，使用的旋转因子依赖于 $s$ 的值，它是 ${\omega }_m$ 的幂，其中 $m=2^s$（引入变量 $m$ 仅为使代码易读）。又引入另一个临时变量 $u$ ，使得能恰当地执行蝴蝶操作。当用循环主体来取代第3行的整个结构时，就得到下面的伪代码，它是稍后将展示的、并行实现的基础。

这一代码首先调用辅助过程 BIT-REVERSE-COPY(a, A) ，把向量 a 按我们所需要的初始顺序复制到数组 A 中。BIT-REVERSE-COPY 如何做到这一点呢？在图30-4中，叶出现的顺序是一个位逆序置换。也就是说，如果让 $rev(k)$ 为「 $k$ 的二进制表示各位逆序」所形成的 $\log n$ 位的整数，那么我们希望把向量中的元素 $a_k$ 放在数组的 $A[rev(k)]$ 位置上。例如，图30-4中叶出现的次序为 $0,4,2,6,1,5,3,7$ ，这个序列用二进制表示为 $000,100,010,110,001,101,011,111$ ，当把二进制表示的各位逆序后，得到序列 $000,001,010,011,100,101,110,111$ 。为了获得一般情况下的位逆序置换，我们注意到，在树的最顶层，最低位为 $0$ 的下标在左子树中、最低位为 $1$ 的下标在右子树中；在每一层去掉最低位后，我们沿着树往下继续这一过程，直到在叶子得出「由位逆序置换给出的顺序」 at the top level of the tree, indices whose low-order bit is 0 go into the left subtree and indices whose low-order bit is 1 go into the right subtree. Stripping off the low-order bit at each level, we continue this process down the tree, until we get the order given by the bit-reversal permutation at the leaves 。由于很容易计算函数 $rev(k)$ ，因此过程 BIT-REVERSE-COPY 相对简单。

调用 BIT-REVERSE-COPY(a, A) 的运行时间当然是 $O(n\log n)$ ，因为迭代了 $n$ 次，并可以在 $O(\log n)$ 时间内、把一个 $0\sim n-1$ 之间的 $\log n$ 位整数逆序（实际上，通常事先知道 $n$ 的初始值，我们就可以编制出一张表、把 $k$ 映射为 $rev(k)$ ，使 BIT-REVERSE-COPY 的运行时间为 $\Theta (n)$ ，且该式中隐含的常数因子较小）。此外，可以采用算导思考题17-1中描述的「聪明的摊还逆序二进制计数器方案」。

这种迭代的FFT实现方法的运行时间为 $\Theta (n\log n)$ 。为了完成 ITERATIVE-FFT 的运行时间是 $\Theta (n\log n)$ 的证明，需要说明最内层循环体（第 $8\sim 13$ 行）执行次数 $L(n)$ 为 $\Theta (n\log n)$ 。对 $s$ 的每个值，第 $6\sim 13$ 行的for循环迭代了 $n/m=n/2^s$ 次、第 $8\sim 13$ 行的最内层循环迭代了 $m/2=2^{s-1}$ 次。因此：$$L(n)=\sum _{s=1}^{\log n}\frac{n}{2^s}\cdot 2^{s-1}=\sum _{s=1}^{\log n}\frac{n}{2}=\Theta (n\log n)$$

4.2 并行FFT电路

我们可以利用许多「允许我们实现一个高效的迭代FFT算法」的性质，来产生一个高效的并行FFT算法，并将其表示成一个电路。图30-5给出了 $n=8$ 时，已知 $n$ 个输入，一个并行FFT电路计算FFT。该电路开始时对输入进行位逆序置换，其后电路分为 $\log n$ 级，每一级由 $n/2$ 个并行执行的蝴蝶操作组成。电路的深度定义为：任意的输入和任意的输出之间、最大的可以达到的计算元素数目，因此上面电路的深度为 $\Theta (\log n)$ 。

并行FFT电路的最左边部分执行位逆序置换，其余部分模拟迭代的 ITERATIVE-FFT 过程。因为最外层for循环的每次迭代执行 $n/2$ 次独立的蝴蝶操作，于是电路并行地执行它们。

• 在 ITERATIVE-FFT 内，每次迭代的值 $s$ 对应于图30-5中的一个蝴蝶操作的阶段 a stage of butterflies 。
• 对于 $s=1,2,\dots ,\log n$ ，阶段 $s$ 有 $n/2^s$ 组蝴蝶操作（对应于 ITERATIVE-FFT 中每个 $k$ 值），每组中有 $2^{s-1}$ 个蝴蝶操作（对应于 ITERATIVE-FFT 中的每个 $j$ 值）。
• 图30-5所示的蝴蝶操作对应于最内层循环的蝴蝶操作（ ITERATIVE-FFT 的第 $9\sim 12$ 行）。此外还要注意，蝴蝶中用到的旋转因子对应于 ITERATIVE-FFT 中用到的那些旋转因子：在阶段 $s$ ，我们使用 ${\omega }_m^0,{\omega }_m^1,\dots ,{\omega }_m^{m/2-1}$ ，其中 $m=2^s$ 。