📘 第 6 章：傅立叶变换

这一篇在干嘛？

从 DFT 到 FFT：Cooley-Tukey、Good-Thomas、Winograd 三大算法与 Rader/Chirp-z 变换，频域处理的完整实现路径。 原书代码为 VHDL，本篇所有代码已改写为 Verilog。

• 第6章 傅立叶变换
• 6.1 傅立叶变换概述
• 6.2 离散傅立叶变换算法
• 例6.1 开窗操作
• 1. 实序列的 DFT
• 算法6.2 用一个 $N$ 点DFT计算长度为 $2N$ 的DFT变换
• 2. 利用 DFT 计算快速卷积
• 算法 6.4 Bluestein Chirp-z 算法
• 例6.5 $N=7$ 的Rader算法
• 例6.6 Rader算法的FPGA实现
• 例6.7 $N=5$ 的WinogradDFT算法
• 6.3 快速傅立叶变换算法
• 算法 6.8 Cooley-Tukey 算法
• 1. 基 $r$ Cooley-Tukey算法
• 2. 基 2 Cooley-Tukey 算法的实现
• 算法6.10 高效的复数乘法器
• 3. Cooley-Tukey 长度 256 FFT 算法的实现
• 例6.11 HDL中的256点FFT
• 定理 6.12 Good-Thomas 索引映射
• 算法 6.13 Good-Thomas FFT 算法
• 例 6.15 使用 Kronecker 乘积且 N=12 的 IDFT
• 例 6.20 生成长度为 256 的 FFT IP
• 6.4 与傅立叶相关的变换
• 6.5 练习

---

自测一下

FFT 为什么把 O(N^2) 降到 O(N logN)？

它利用旋转因子的对称性与周期性，把 N 点 DFT 分解成两个 N/2 点 DFT（奇偶分治），逐级减半，共 log2N 级、每级 N/2 个蝶形。

DIF 与 DIT 蝶形有何区别？

DIT（按时间抽取）输入是位反序、输出自然序；DIF（按频率抽取）相反。两者运算量相同，选哪个取决于希望哪端有序（方便串行读写）。

Goertzel 算法适合什么场景？

只需要少数几个频点的谱线时（如 DTMF 双音检测），Goertzel 相当于单点 DFT 的二阶 IIR，比整块 FFT 省资源。

实现 FFT 时位反序怎么处理？

三种方案：输入位反序预排序；输出位反序再用位反序计数器读出；或用”自然输入-自然输出”算法（Good-Thomas 素因子或双基 2

6.1 傅立叶变换概述

离散傅立叶变换（Discrete Fourier Transform, DFT）以及它的快速实现——快速傅立叶变换（Fast Fourier Transform, FFT），是数字信号处理的核心工具之一。简单说：DFT 把一段有限长的采样信号，从”时间的样子”翻译成”频率的样子”，告诉我们在哪些频率上有多少能量。没有它，频谱分析、滤波、图像压缩等任务在计算机和 FPGA 上都无从下手。

DFT 和 FFT 算法已经被”发明（和再发明）“了很多次。Heideman 等人${}^{[174]}$指出，高斯早在计算机出现之前就使用过一种与今天的 Cooley-Tukey FFT 同类的算法。本章讨论图 6-1 中总结的最重要的一些算法。

图 6-1 DFT 和 FFT 算法的分类

分类术语沿用 Burrus${}^{[175]}$的做法：按照输入序列与输出序列之间的（多维）索引映射关系来划分。凡是没有使用多维索引映射的算法统称为 DFT 算法（其中有些算法确实包含必不可少的少量计算，如 Winograd DFT 算法）。需要注意的是，DFT 算法和 FFT 算法并非彼此孤立：实际中的高效实现往往是二者组合的结果。例如 Rader 质数算法与 Good-Thomas FFT 的组合就产生了著名的 VLSI 实现。文献${}^{[176-181]}$给出了用 PDSP 和 ASIC 实现 FFT 的许多设计示例，而 FPGA 已经可以实现一维${}^{[182-184]}$和二维${}^{[52和185]}$的 FFT 变换。

本章将讨论 4 种最重要的 DFT 算法和 3 种最常用的 FFT 算法，并按计算量比较它们的实现问题；最后还会介绍与傅立叶相关的变换（如 DCT），它是 JPEG、MPEG 等图像压缩标准的重要工具。若想深入钻研，可参阅基础 DSP 书籍${}^{[5,88,186,187]}$和专门的 FFT 书籍${}^{[74,188-192]}$。

6.2 离散傅立叶变换算法

这一节先复习 DFT 的关键性质，再介绍 Bluestein、Goertzel、Rader 和 Winograd 提出的几种基本 DFT 算法。

6.2.1 用 DFT 近似傅立叶变换

连续时间的傅立叶变换对定义为：

$$X(f)={\int }_{-\infty }^{\infty }x(t)e^{-j2\pi ft}\mathrm{d}t\leftrightarrow x(t)={\int }_{-\infty }^{\infty }X(f)e^{j2\pi ft}\mathrm{d}f\text{\hspace{1em}(6-1)}$$

这个公式假设信号在时间和频率上都是连续、无限持续的。但真实系统只能处理有限个、量化过的采样点，于是就有了 DFT：在时域和频域各取 $N$ 次采样：

$$X[k]=\sum _{n=0}^{N-1}x[n]e^{-j2\pi kn/N}=\sum _{n=0}^{N-1}x[n]W_N^{kn}\text{\hspace{1em}(6-2)}$$

逆变换 IDFT 为：

$$x[n]=\frac{1}{N}\sum _{k=0}^{N-1}X[k]e^{j2\pi kn/N}=\frac{1}{N}\sum _{k=0}^{N-1}X[k]W_N^{-kn}\text{\hspace{1em}(6-3)}$$

写成向量/矩阵形式就是：

$$X=Wx\leftrightarrow x=\frac{1}{N}{W}^{\ast }X\text{\hspace{1em}(6-4)}$$

用 DFT 去近似傅立叶频谱时，必须记住”采样”这个动作本身带来的两个影响——否则结果会悄悄失真：

• 时域采样使频谱变成周期性的（周期为采样频率 $f_S$）。按”香农采样定理”，只有当 $x(t)$ 的频率分量全部集中在低于奈奎斯特频率 $f_S/2$ 的范围内时，DFT 才能合理地近似傅立叶变换；否则会出现混叠。
• 频域采样使时间函数变成周期性的，也就是说 DFT 默认时间序列是周期性的。如果信号在 $N$ 点窗口内没有完成整数个周期，就会产生”泄漏”（leakage）现象。因此如果 $x(t)$ 是周期信号，应尽量选择覆盖整数个周期的采样频率和分析窗口。

更实用的抗泄漏手段是使用窗函数：给数据乘上一个两边逐渐衰减为 0 的权函数。这类窗函数在第 3 章 FIR 滤波器设计上下文中已经讨论过（表 3-2）。图 6-2 给出了几种典型窗函数的时间和频率特性${}^{[128,193]}$。

典型窗函数的时域形状

图 6-2 时域和频域内的窗函数

例 6.1 开窗操作

图 6-3(a) 给出了一个在采样窗口内不能完成整数个周期的正弦信号。它的理想傅立叶变换应该只在 $\pm {\omega }_0$ 处有两个单位脉冲函数，如图 6-3(b) 所示。图 6-3(c) 和 (d) 分别是采用框函数（即不加窗）和汉宁窗做 DFT 分析的结果。可以直观看到：框函数的旁瓣纹波明显更多；而汉宁窗虽然把泄漏压下去了，但主瓣宽度变宽了——频率分辨率略有牺牲。这就是加窗的基本折衷：用主瓣变宽换旁瓣降低。

图 6-3(a) DFT 窗口

图 6-3(b) 理想傅立叶变换

图 6-3(c) 采用框函数的 DFT

图 6-3(d) 采用汉宁窗的 DFT

6.2.2 DFT 的性质

表 6-1 总结了 DFT 最重要的性质。许多性质与连续傅立叶变换一致：变换是唯一（双射）的、满足叠加原理、实部与虚部通过希尔伯特变换相联系。

表 6-1 DFT 定理

定理	$x(n)$	$X(k)$
变换	$x(n)$	${\sum }_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}$
逆变换	$\frac{1}{N}{\sum }_{n=0}^{N-1}X[k]e^{j2\pi nk/N}$	$X[k]$
叠加	$s_1{x}_1[n]+s_2{x}_2[n]$	$s_1{X}_1[k]+s_2{X}_2[k]$
时间反向	$x[-n]$	$X[-k]$
共轭复数	$x^{\ast }[n]$	$X^{\ast }[-k]$
实部	$\Re (x[n])$	$(X[k]+X^{\ast }[-k])/2$
虚部	$\Im (x[n])$	$(X[k]+X^{\ast }[-k])/(2j)$
实偶数部分	$x_e[n]=(x[n]+x[-n])/2$	$\Re (X[k])$
实奇数部分	$x_o[n]=(x[n]-x[-n])/2$	$j\Im (X[k])$
对称性	$X[n]$	$Nx[-k]$
循环卷积	$x[n]\oplus f[n]$	$X[k]F[k]$
乘法	$x[n]\times f[n]$	$\frac{1}{N}X[k]\oplus F[k]$
周期平移	$x[n-d\mathrm{mod}N]$	$X[k]e^{-j2\pi dk/N}$
帕斯瓦尔定理	${\sum }_{n=0}^{N-1}|x[n]{|}^2$	$\frac{1}{N}{\sum }_{k=0}^{N-1}|X[k]{|}^2$

正、逆变换结构的相似性还提供了一个巧妙的”反演技巧”。由向量/矩阵表示（式 6-5）：

$$X=Wx\leftrightarrow x=\frac{1}{N}{W}^{\ast }X\text{\hspace{1em}(6-5)}$$

两边取共轭可得：

$${\mathbf{x}}^{\ast }=\frac{1}{N}({\mathbf{W}}^{\ast }\mathbf{X}{)}^{\ast }=\frac{1}{N}\mathbf{W}{\mathbf{X}}^{\ast }\text{\hspace{1em}(6-6)}$$

也就是说：只要对 $X^{\ast }$ 做一次正向 DFT 再除以 $N$，就得到了逆 DFT。这样硬件上只需实现一个方向的变换核。

1. 实序列的 DFT

工程中输入往往是实数序列，而 DFT 本身按复数设计，这就有可省的计算。利用”实序列的频谱具有偶对称的实部和奇对称的虚部”这一性质（表 6-1），可以合成两个经典技巧${}^{[188]}$：要么用一个 $N$ 点 DFT 同时算两个 $N$ 点实序列，要么用一个 $N$ 点 DFT 算一个长度 $2N$ 的实序列。

算法 6.2　用一个 $N$ 点 DFT 计算长度为 $2N$ 的 DFT

由时间序列 $x[n]$ 计算 $2N$ 点 DFT $X[k]=X_{\mathrm{r}}[k]+j{X}_{\mathrm{i}}[k]$：

1. 构造 $N$ 点复序列 $y[n]=x[2n]+\mathrm{j}x[2n+1]$，$n=0,1,\dots ,N-1$（偶数位放实部、奇数位放虚部）。
2. 计算 $y[n]\circ -\bullet Y[k]=Y_{\mathrm{r}}[k]+j{Y}_{\mathrm{i}}[k]$。
3. 最后组合：

$$X_{\mathrm{r}}[k]=\frac{Y_{\mathrm{r}}[k]+Y_{\mathrm{r}}[-k]}{2}+\cos (\pi k/N)\frac{Y_{\mathrm{i}}[k]+Y_{\mathrm{i}}[-k]}{2}-\sin (\pi k/N)\frac{Y_{\mathrm{r}}[k]-Y_{\mathrm{r}}[-k]}{2}$$ $$X_{\mathrm{i}}[k]=\frac{Y_{\mathrm{i}}[k]-Y_{\mathrm{i}}[-k]}{2}-\sin (\pi k/N)\frac{Y_{\mathrm{i}}[k]+Y_{\mathrm{i}}[-k]}{2}-\cos (\pi k/N)\frac{Y_{\mathrm{r}}[k]-Y_{\mathrm{r}}[-k]}{2}$$

其中 $k=0,1,\dots ,N-1$。除了一次 $N$ 点 DFT（或 FFT）之外，额外开销只有来自旋转因子 $\pm \exp (j\pi k/N)$ 的 $4N$ 次实数加法和乘法——用一半规模的变换算出了两倍长度的谱。

算法 6.3　用一个 $N$ 点 DFT 计算两个长度为 $N$ 的 DFT

要同时计算 $g[n]\circ -\bullet G(k)$ 和 $h[n]\circ -\bullet H(k)$：

1. 构造 $N$ 点复序列 $y[n]=h[n]+\mathrm{j}g[n]$，$n=0,1,\dots ,N-1$。
2. 计算 $y[n]\circ -\bullet Y[k]=Y_{\mathrm{r}}[k]+j{Y}_{\mathrm{i}}[k]$。
3. 最后分离：

$$H[k]=\frac{Y_{\mathrm{r}}[k]+Y_{\mathrm{r}}[-k]}{2}+\mathrm{j}\frac{Y_{\mathrm{i}}[k]-Y_{\mathrm{i}}[-k]}{2},G[k]=\frac{Y_{\mathrm{i}}[k]+Y_{\mathrm{i}}[-k]}{2}-\mathrm{j}\frac{Y_{\mathrm{r}}[k]-Y_{\mathrm{r}}[-k]}{2}$$

其中 $k=0,1,\dots ,N-1$。原理是：实序列的频谱共轭对称，纯虚序列的频谱反对称，两者叠加后仍可无损拆开。除一次 $N$ 点 DFT 外，额外开销只有 $2N$ 次实数加法。

2. 利用 DFT 计算快速卷积

DFT（或 FFT）最常见的用途之一是计算卷积：时域卷积对应频域相乘。先把两个序列变换到频域，做标量点积，再变换回时域。但要注意关键区别：DFT 算的是循环卷积，而不是线性卷积。用 FFT 做快速卷积时必须处理这一点，于是产生了两种经典方法：

• 重叠保留（overlap-save）：直接做循环卷积，然后简单放弃边界处被”绕回”污染的采样点；
• 重叠相加（overlap-add）：对信号和滤波器补 0，再把各段部分结果直接相加拼接成完整输出。

快速卷积的输入通常是实序列，因此可以用实变换（如练习 6.15 讨论的 Hartley 变换）高效完成。Hartley 变换也能构造类 FFT 算法，与复变换${}^{[194]}$相比性能可提高两倍。如果手头只有复数 FFT 程序，可以套用算法 6.2 或 6.3。图 6-4 给出一种与算法 6.3 类似的替代方案：用一个 $N$ 点 DFT 实现两个 $N$ 点变换——实部用作 DFT，虚部用作 IDFT，根据卷积理论，逆变换时需要用到虚部。

图 6-4 采用复数 FFT 的实数卷积${}^{[74]}$

此外，若实值滤波器的 DFT（满足 $F[k]=F[-k{]}^{\ast }$）已经离线算好，频域内只需要 $N/2$ 次乘法即可完成 $X[k]F[k]$。

6.2.3 Goertzel 算法

如果只需要频谱中的某一个分量 $X[k]$，用整套 DFT 就太浪费了。把 DFT 的单项展开：

$$X[k]=x[0]+x[1]W_N^k+x[2]W_N^{2k}+\cdots +x[N-1]W_N^{(N-1)k}$$

将所有 $x[n]$ 用公共因子 $W_N^k$ 逐步括起来，就得到 Horner 嵌套形式：

$$X[k]=x[0]+W_N^k\left(x[1]+W_N^k\left(x[2]+\cdots +W_N^{k}x[N-1]\right)\dots \right)$$

这是一个可行的递归计算——这就是 Goertzel 算法，图 6-5 给出了图形解释。递归从输入序列的最后一个值 $x[N-1]$ 开始计算，经过 $N$ 步递推后，输出端就得到 $X[k]$ 这一个频谱值。

图 6-5 长度为 4 的 Goertzel 算法

如果需要计算好几个相邻的频谱分量，还可以把 $e^{\pm j2\pi n/N}$ 类型的因子两两组合，得到一个带分母的二阶系统：

$$z^2-2z\cos \left(\frac{2\pi n}{N}\right)+1$$

这样一来，所有的复数乘法都简化成了实数乘法，硬件代价大幅下降。

什么时候用 Goertzel？判断标准很简单：若只需要少量频谱分量，Goertzel 非常划算；若要算整个 DFT，计算量仍是 $N^2$ 数量级，与直接计算相比就没有任何优势了——此时该换 FFT。

6.2.4 Bluestein Chirp-z 变换

Bluestein Chirp-z 变换（CZT）的核心技巧是把 DFT 指数 $nk$ 二次展开：

$$nk=-(k-n{)}^2/2+n^2/2+k^2/2\text{\hspace{1em}(6-7)}$$

代入 DFT 定义后，原本”每个输出都要乘遍所有输入”的 $nk$ 交叉项消失了，DFT 变成：

$$X[k]=W_N^{k^2/2}\sum _{n=0}^{N-1}\left(x[n]W_N^{n^2/2}\right)W_N^{-(k-n{)}^2/2}\text{\hspace{1em}(6-8)}$$

这个求和正是一个线性卷积。图 6-6 给出了算法的图形解释，于是得到：

图 6-6 Bluestein Chirp-z 算法

算法 6.4　Bluestein Chirp-z 算法，DFT 的计算分为 3 步：

1. $x[n]$ 与 $W_N^{n^2/2}$ 的 $N$ 次乘法（预调制）；
2. $x[n]W_N^{n^2/2}$ 与 $W_N^{n^2/2}$ 的线性卷积（核心步骤，可用 FIR 滤波器实现）；
3. 结果再与 $W_N^{k^2/2}$ 的 $N$ 次乘法（后调制）。

一次完整变换需要一个长度为 $N$ 的卷积加 $2N$ 次复数乘法。与后面要讲的 Rader 算法相比，CZT 最大的优点是：变换长度 $N$ 不必是质数，可以定义成任意长度。

Narasimha${}^{[195]}$等人进一步注意到：CZT 中 FIR 滤波器部分的许多系数其实是无关紧要或彼此相同的。例如长度为 8 的 CZT，其 FIR 滤波器长度为 14，但图 6-7 中只有 4 个不同的复系数——分别是 $1$、$\mathrm{j}$ 和 $\pm e^{\mathrm{j}22.5^{\circ }}$，也就是说真正需要实现的不可或缺的实系数只有 2 个。

图 6-7 CZT 系数 $C(n)=e^{\mathrm{j}2\pi \frac{n^2/2\mathrm{mod}8}{8}}$，$n=1,2,\dots ,14$

那么固定系数位宽 $C_N$ 下，DFT 能做到多长？表 6-2 给出了按复系数总数 $C_N$ 统计的数据。

表 6-2 DFT 的最大长度

DFT 的长度	8	12	16	24	40	48	72	80	120	144	168	180	240	360	504
$C_N$	4	6	7	8	12	14	16	21	24	28	32	36	42	48	64

不过，复系数的数量与实际计算量之间并没有直接对应关系，因为其中一些系数可能是平凡的（如 $\pm 1$ 或 $\pm \mathrm{j}$）或对称的。特别是 2 的幂对应的长度具有大量对称性，如图 6-8 所示。若按”不可或缺的重要实系数”个数来统计，就会得到表 6-3 的最大长度变换。

图 6-8 CZT 的复系数和重要的实数乘法的数量

表 6-3 具体数量的重要实系数的 DFT 的最大长度

DFT 的长度	10	16	20	32	40	48	50	80	96	160	192
sin/cos	2	3	5	6	8	9	10	11	14	20	25

因此，长度 16 和 32 是分别只需要 3 个和 6 个实数乘法器的最大长度——这也是它们在硬件实现中格外受欢迎的原因。

一般情况下，2 的幂长度是受欢迎的 FFT 构造模块。表 6-4 给出了以转置形式实现长度 $N=2^n$ 的 CZT 滤波器时的工作量。

表 6-4 在转置形式中实现长度为 $N=2^n$ 的 CZT 滤波器的工作量

$N$	$C_N$	sin/cos	CSD 加法器	RAG 加法器	14 位系数的 NOF
8	4	2	23	7	3,5,33,49,59
16	7	3	91	8	3,25,59,63,387
32	12	6	183	13	2,25,49,73,121,375
64	23	11	431	18	5,25,27,93,181,251,7393
128	44	22	879	31	5,15,25,175,199,319,403,499,1567
256	87	42	1911	49	5,25,765,1443,1737,2837,4637

逐列解读：第 1 列是 DFT 长度 $N$；第 2 列是复指数系数总数 $C_N$，最坏情况下需要实现 $2C_N$ 个实系数；第 3 列给出实际不同重要实系数的数量——对比第 2、3 列可见，对称系数和平凡系数显著减少了工作量。最后三列是用 15 位（14 个无符号位加 1 个符号位）系数精度、分别采用（第 2 章讨论过的）CSD 算法和 RAG 算法实现时所需的加法器数量，以及 RAG 使用的辅助系数（NOF）。CSD 编码无法利用系数对称性，加法器较多；而 RAG 算法能从根本上减少工作量——对本表长度至少减少 16 倍以上。

6.2.5 Rader 算法

Rader 算法${}^{[196,197]}$只能定义质数长度 $N$ 的 DFT：

$$X[k]=\sum _{n=0}^{N-1}x[n]W_N^{nk}k,n\in Z_N;\mathrm{ord}(W_N)=N\text{\hspace{1em}(6-9)}$$

它的思路分两步。第一步单独算出 DC 分量：

$$X[0]=\sum _{n=0}^{N-1}x[n]\text{\hspace{1em}(6-10)}$$

第二步处理剩下的部分。由于 $N=p$ 是质数，根据第 2 章的讨论，存在本原元素（生成器）$g$，它的幂可以取遍 $Z_p$ 中除 0 之外的所有元素。用 $g^n\mathrm{mod}N$ 代替 $n$、$g^k\mathrm{mod}N$ 代替 $k$，就得到索引变换：

$$X[g^k\mathrm{mod}N]-x[0]=\sum _{n=0}^{N-2}x[g^n\mathrm{mod}N]W_N^{g^{n+k}\mathrm{mod}N}\text{\hspace{1em}(6-11)}$$

其中 $k\in \{1,2,3,\dots ,N-1\}$。注意式 (6-11) 右端正是一个循环卷积：

$$\left[x[g^0\mathrm{mod}N],x[g^1\mathrm{mod}N],\dots ,x[g^{N-2}\mathrm{mod}N]\right]\ast \left[W_N,W_N^g,\dots ,W_N^{g^{N-2}\mathrm{mod}N}\right]\text{\hspace{1em}(6-12)}$$

也就是说：Rader 把质数长度的 DFT 变成了一个循环卷积，而卷积可以用 FIR 滤波器硬件高效实现。下面通过 $N=7$ 的例子把全过程算一遍。

例 6.5 $N=7$ 的 Rader 算法

对 $N=7$，取生成器 $g=3$（可参阅文献${}^{[5]}$表 B-7），其幂对 7 取模构成索引变换：

$$[3^0,3^1,3^2,3^3,3^4,3^5]\mathrm{mod}7\equiv [1,3,2,6,4,5]\text{\hspace{1em}(6-13)}$$

首先计算 DC 分量：

$$X[0]=\sum _{n=0}^{6}x[n]=x[0]+x[1]+x[2]+x[3]+x[4]+x[5]+x[6]$$

然后计算 $X[k]-x[0]$ 的循环卷积：

$$\left[x[1],x[3],x[2],x[6],x[4],x[5]\right]\odot \left[W_7,W_7^3,W_7^2,W_7^6,W_7^4,W_7^5\right]$$

写成矩阵形式：

$$\left[\begin{array}{l}X[1]\\ X[3]\\ X[2]\\ X[6]\\ X[4]\\ X[5]\end{array}\right]=\left[\begin{array}{cccccc}{W}_7^1& W_7^3& W_7^2& W_7^6& W_7^4& W_7^5\\ W_7^3& W_7^2& W_7^6& W_7^4& W_7^5& W_7^1\\ W_7^2& W_7^6& W_7^4& W_7^5& W_7^1& W_7^3\\ W_7^6& W_7^4& W_7^5& W_7^1& W_7^3& W_7^2\\ W_7^4& W_7^5& W_7^1& W_7^3& W_7^2& W_7^6\\ W_7^5& W_7^1& W_7^3& W_7^2& W_7^6& W_7^4\end{array}\right]\left[\begin{array}{l}x[1]\\ x[3]\\ x[2]\\ x[6]\\ x[4]\\ x[5]\end{array}\right]+\left[\begin{array}{l}x[0]\\ x[0]\\ x[0]\\ x[0]\\ x[0]\\ x[0]\end{array}\right]\text{\hspace{1em}(6-14)}$$

注意这个矩阵的每一行都是上一行的循环移位——这正是”循环卷积”的结构特征。图 6-9 给出了用 FIR 滤波器实现的图形解释。

图 6-9 长度为 $p=7$ 的 Rader 质数因子 DFT 实现

用一个三角形信号 $x[n]=10n$（步长为 10）来验证上式。直接代入式 (6-14)：

$$\left[\begin{array}{l}X[1]\\ X[3]\\ X[2]\\ X[6]\\ X[4]\\ X[5]\end{array}\right]=\left[\begin{array}{cccccc}{W}_7^1& W_7^3& W_7^2& W_7^6& W_7^4& W_7^5\\ W_7^3& W_7^2& W_7^6& W_7^4& W_7^5& W_7^1\\ W_7^2& W_7^6& W_7^4& W_7^5& W_7^1& W_7^3\\ W_7^6& W_7^4& W_7^5& W_7^1& W_7^3& W_7^2\\ W_7^4& W_7^5& W_7^1& W_7^3& W_7^2& W_7^6\\ W_7^5& W_7^1& W_7^3& W_7^2& W_7^6& W_7^4\end{array}\right]\left[\begin{array}{l}20\\ 40\\ 30\\ 70\\ 50\\ 60\end{array}\right]+\left[\begin{array}{l}10\\ 10\\ 10\\ 10\\ 10\\ 10\end{array}\right]=\left[\begin{array}{l}-35+j72\\ -35+j8\\ -35+j28\\ -35-j72\\ -35-j8\\ -35-j28\end{array}\right]$$

而 $X[0]$ 是时间序列之和：$10+20+\dots +70=280$。可以检查：6 个非 DC 输出的实部都是 $-35$，虚部按 $\pm$ 成对出现——这正是实输入序列频谱共轭对称的体现。

在 Rader 算法中还可以利用复数对 $e^{\pm \mathrm{j}2k\pi /N}$，$k\in [0,N/2]$ 的对称性，构造更高效的 FIR 实现（练习 6.5）。既然 Rader 质数因子 DFT 等价于一个 FIR 滤波器，那么第 3 章的所有快速 FIR 技术——完全流水线 DA、RAG 转置结构——都可以直接搬过来用。下面给出一个 RAG 的 FPGA 实现示例。

例 6.6 Rader 算法的 FPGA 实现

长度为 7 的 Rader 算法的 RAG 实现过程如下。第一步是系数量化：假定输入和系数都用有符号 8 位字表示（缩放因子 256），量化后的系数如表 6-5 所示。

表 6-5 量化后的系数

$k=$	0	1	2	3	4	5	6
实部 $\{256\times W_7^k\}$	256	160	−57	−231	−231	−57	160
虚部 $\{256\times W_7^k\}$	0	−200	−250	−111	111	250	200

如果所有系数都用直接形式的常系数乘法器实现（参见表 2-3），共需要 24 个加法器。改用转置结构后，利用若干系数仅符号不同这一事实，每个系数的工作量可降到 11 个加法器。进一步做 RAG 优化（参阅图 2-4），加法器数量可以压缩到最小值 7——相比直接 FIR 体系结构有 3 倍以上的提高。

原书 VHDL 改写为 Verilog。下面的 Verilog 代码给出了采用转置 FIR 结构、长度为 7 的 Rader DFT 的一个可行实现：

// 原书 VHDL 改写为 Verilog：7 点 Rader 质数因子 DFT（转置 FIR + RAG）
module rader7 (
    input  wire               clk,     // 系统时钟
    input  wire               reset,   // 异步复位
    input  wire signed [7:0]  x_in,    // 实数输入
    output reg  signed [10:0] y_real,  // 实部输出
    output reg  signed [10:0] y_imag   // 虚部输出
);
 
    // ---------- 状态机状态定义 ----------
    localparam [1:0] START = 2'd0,   // 初始化
                     LOAD  = 2'd1,   // 装载阶段：算 X[0]
                     RUN   = 2'd2;   // 运行阶段：FIR 输出各 X[k]
 
    reg [1:0]  state;                 // 状态变量
    reg [3:0]  count;                 // 时钟周期计数器
    reg signed [17:0] accu;           // X[0] 累加器
    reg signed [7:0]  x, x_0;         // 延迟输入与 x[0] 暂存
 
    // ---------- 转置 FIR 抽头延迟线（实部/虚部各 6 级） ----------
    reg signed [17:0] real_part [0:5];
    reg signed [17:0] imag_part [0:5];
 
    // ---------- RAG 系数与辅助因子 ----------
    reg  signed [17:0] x57, x111, x160, x200, x231, x250;
    reg  signed [17:0] x5, x25, x110, x125;
    wire signed [17:0] xs   = {{10{x_in[7]}}, x_in}; // 输入扩展到 18 位
    wire signed [17:0] x256 = xs <<< 8;              // 256*x，仅移位
 
    // ---------- 状态机：区分 Start / Load / Run 三个处理阶段 ----------
    always @(posedge clk or posedge reset) begin
        if (reset) begin                 // 异步复位
            state <= START; accu <= 0; count <= 0;
            y_real <= 0;    y_imag <= 0;
        end else begin
            case (state)
                START: begin             // 初始化
                    state  <= LOAD;
                    count  <= 1;
                    x_0    <= x_in;      // 保存 x[0]
                    accu   <= 0;         // X[0] 累加器清零
                    y_real <= 0;
                    y_imag <= 0;
                end
                LOAD: begin              // 依次输入 x[5],x[4],x[6],x[2],x[3],x[1]
                    if (count == 8)      // 装载阶段结束？
                        state <= RUN;
                    else begin
                        state <= LOAD;
                        accu  <= accu + x; // 累加求 X[0]
                    end
                    count <= count + 1;
                end
                RUN: begin               // 再次输入 x[5],x[4],x[6],x[2],x[3]
                    if (count == 15) begin // 运行阶段结束？
                        y_real <= accu;    // 输出 X[0]
                        y_imag <= 0;       // 只有实输入，Im(X[0]) = 0
                        state  <= START;   // 输出结果，重新开始
                        count  <= 0;
                    end else begin
                        y_real <= (real_part[0] >>> 8) + {{3{x_0[7]}}, x_0};
                        y_imag <=  imag_part[0] >>> 8;
                        state  <= RUN;
                        count  <= count + 1;
                    end
                end
            endcase
        end
    end
 
    // ---------- 转置结构的两条 FIR 滤波器通路 ----------
    integer k;
    always @(posedge clk or posedge reset) begin
        if (reset) begin                 // 异步清零
            for (k = 0; k < 6; k = k + 1) begin
                real_part[k] <= 0;
                imag_part[k] <= 0;
            end
            x <= 0;
        end else begin
            x <= x_in;
            // 实部通路（转置 FIR，系数为 256*W7^k 的实部）
            real_part[0] <= real_part[1] + x160;  // W^1
            real_part[1] <= real_part[2] - x231;  // W^3
            real_part[2] <= real_part[3] - x57;   // W^2
            real_part[3] <= real_part[4] + x160;  // W^6
            real_part[4] <= real_part[5] - x231;  // W^4
            real_part[5] <= -x57;                 // W^5
            // 虚部通路（转置 FIR，系数为 256*W7^k 的虚部）
            imag_part[0] <= imag_part[1] - x200;  // W^1
            imag_part[1] <= imag_part[2] - x111;  // W^3
            imag_part[2] <= imag_part[3] - x250;  // W^2
            imag_part[3] <= imag_part[4] + x200;  // W^6
            imag_part[4] <= imag_part[5] + x111;  // W^4
            imag_part[5] <= x250;                 // W^5
        end
    end
 
    // ---------- RAG 系数合成（寄存一拍，保证时序） ----------
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            x160 <= 0; x200 <= 0; x250 <= 0;
            x57  <= 0; x111 <= 0; x231 <= 0;
        end else begin
            x160 <= x5  <<< 5;         // 32 * 5x   = 160x
            x200 <= x25 <<< 3;         // 8  * 25x  = 200x
            x250 <= x125 <<< 1;        // 2  * 125x = 250x
            x57  <= x25 + (xs <<< 5);  // 25x + 32x = 57x
            x111 <= x110 + xs;         // 110x + x  = 111x
            x231 <= x256 - x25;        // 256x - 25x = 231x
        end
    end
 
    // ---------- RAG 辅助因子（组合逻辑，不占用触发器） ----------
    always @* begin
        x5   = (xs  <<< 2) + xs;             // 5x
        x25  = (x5  <<< 2) + x5;             // 25x
        x110 = (x25 <<< 2) + (x5  <<< 1);    // 110x
        x125 = (x25 <<< 2) + x25;            // 125x
    end
 
endmodule

这个设计由 4 个 always 块组成。第一个是状态机，区分 Start、Load、Run 三个处理阶段：Start 阶段保存 $x[0]$ 并清零累加器；Load 阶段把置换顺序的输入依次送入，同时累加出 $X[0]$；Run 阶段让 FIR 滤波器输出各频谱分量。第二个 always 块定义实部和虚部两条转置 FIR 滤波器通路。第三个块用 RAG 方式合成系数。第四个块计算 RAG 算法的未寄存辅助因子——可以看到，所有系数最终只由 6 个加法器和 1 个减法器（再加移位）实现，完全不需要乘法器。本设计使用了 443 个 LE，未使用嵌入式乘法器；采用 TimeQuest 缓慢 85C 模型时，时序性能 $F_{\max }=138.45\mathrm{MHz}$。图 6-10 给出了 Quartus II 对三角形输入序列 $x[n]=\{10,20,30,40,50,60,70\}$ 的仿真结果。注意：输入和输出序列从 1 μs 开始、按置换后的顺序出现（对应索引变换 $[1,3,2,6,4,5]$）；如果仿真器使用有符号数据类型，负的结果会带负号显示；最后在 1.7 μs 处 $X[0]=280$ 被送到输出端，电路随即准备好处理下一个输入帧——数值与例 6.5 的手算结果完全一致。

图 6-10 7 点 Rader 算法的仿真结果

由于 Rader 算法受限于质数长度，与 CZT 相比系数缺乏对称性。表 6-6 给出了质数长度 $2^n\pm 1$ 时，以转置形式实现循环滤波器所需的工作量。

表 6-6 实现转置形式的循环滤波器所需的工作量

DFT 的长度	sin/cos	CSD 加法器	RAG 加法器	14 位系数的 NOF
7	6	52	13	7, 11, 31, 59, 101, 177, 319
17	16	138	23	3, 35, 103, 415, 1153, 1249, 8051
31	30	244	38	3, 9, 133, 797, 877, 975, 1179, 3235
61	60	496	66	5, 39, 51, 205, 265, 3211
127	124	1060	126	5

第 1 列是循环卷积长度 $N$（即复系数数量）。将第 2 列与最坏情况 $2N$ 个实 sin/cos 系数相比，可见对称性和平凡系数已经把重要系数减半。接下来两列分别是用 CSD 和 RAG 算法、14 位（加符号位）系数精度实现时所需的加法器数量，最后一列是 RAG 所用的辅助系数 NOF。规律很清晰：CSD 类型滤波器的工作量约为 $BN/2$（$B$ 为系数位宽，此处 14）；而 RAG 的工作量仅为 $N$ 量级——对长滤波器提高了约 $B/2=7$ 倍。原因是 RAG 已经合成的系数生成了一个”密集的”小系数网格，每多一个系数只需再增加约一个加法器。

6.2.6 Winograd DFT 算法

第一种以”精简必要乘法数量”为目标的算法是 Winograd DFT 算法。它是 Rader 算法（把 DFT 转换成循环卷积）与前面实现快速 FIR 时用过的 Winograd${}^{[124]}$短卷积算法（参阅 5.3.2 节）的结合。因此其长度被限制在质数或质数幂范围内。表 6-7 简要列出了所需运算量。

表 6-7 带有实输入的 Winograd DFT 的工作量（无关紧要的乘法指乘以 $\pm 1$ 或 $\pm \mathrm{j}$；复数输入时运算量翻倍）

块长度	实数乘法总数	重要乘法总数	实数加法总数
2	2	0	2
3	3	2	6
4	4	0	8
5	6	5	17
7	9	8	36
8	8	2	26
9	11	10	44
11	21	20	84
13	21	20	94
16	18	10	74
17	36	35	157
19	39	38	186

下面用 $N=5$ 的例子详细演示 Winograd DFT 算法的构造步骤。

例 6.7 $N=5$ 的 Winograd DFT 算法

采用文献${}^{[5]}$给出的 Rader 算法的另一种表示形式（用 $X[0]$ 代替 $x[0]$ 的位置）：

$$X[0]=\sum _{n=0}^{4}x[n]=x[0]+x[1]+x[2]+x[3]+x[4]$$ $$X[k]-X[0]=\left[x[1],x[2],x[4],x[3]\right]\odot \left[W_5-1,W_5^2-1,W_5^4-1,W_5^3-1\right]k=1,2,3,4$$

若用 Winograd 算法实现这个长度为 4 的循环卷积，只需 5 次重要乘法，得到下面的完整算法：

$$X[k]=\sum _{n=0}^{4}x[n]e^{-j2\pi kn/5}k=0,1,\dots ,4$$ $$\left[\begin{array}{l}X[0]\\ X[4]\\ X[3]\\ X[2]\\ X[1]\end{array}\right]=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 1& 1& 1& 1& 0& -1\\ 1& 1& -1& 1& 1& 0\\ 1& 1& -1& -1& -1& 0\\ 1& 1& 1& -1& 0& 1\end{array}\right]$$ $$\times \mathrm{diag}\left(1,\frac{1}{2}(\cos (2\pi /5)+\cos (4\pi /5))-1,\frac{1}{2}(\cos (2\pi /5)-\cos (4\pi /5)),\mathrm{j}\sin (2\pi /5),\mathrm{j}(-\sin (2\pi /5)+\sin (4\pi /5)),\mathrm{j}(\sin (2\pi /5)+\sin (4\pi /5))\right)$$ $$\times \left[\begin{array}{ccccc}1& 1& 1& 1& 1\\ 0& 1& 1& 1& 1\\ 0& 1& -1& -1& 1\\ 0& 1& -1& 1& -1\\ 0& 1& 0& 0& -1\\ 0& 0& -1& 1& 0\end{array}\right]\left[\begin{array}{l}x[0]\\ x[1]\\ x[2]\\ x[3]\\ x[4]\end{array}\right]$$

读法：最右侧矩阵做输入加法，中间对角阵放全部乘法，最左侧矩阵做输出加法。对实数（或虚数）输入序列 $x[n]$，总计算量分别只有 5 次（或 10 次）重要实数乘法。图 6-11 的信号流程图还展示了如何用有效的方式安排加法。

图 6-11 Winograd 5 点 DFT 信号流程图

用矩阵表示 Winograd DFT 算法非常方便：

$${\mathbf{W}}_{N_l}={\mathbf{C}}_l{\mathbf{B}}_l{\mathbf{A}}_l\text{\hspace{1em}(6-15)}$$

其中 $A_l$ 合并输入加法，$B_l$ 是傅立叶系数的对角矩阵，$C_l$ 包括输出加法。唯一的缺点是：短卷积算法中输入加法、输出加法的精确执行顺序在矩阵表达式中消失了，不容易直接排布时序。

这一”Rader 算法 + 短 Winograd 卷积”的组合就是 Winograd DFT 算法；后面它将与索引映射结合，构成 Winograd FFT 算法。它是目前所有已知 FFT 算法中实数乘法次数最少的算法。

6.3 快速傅立叶变换算法

如 6.1 节所述，这里沿用 Burrus${}^{[175]}$的术语：所有 FFT 算法都按照输入序列与输出序列之间不同的（多维）索引映射来分类。这建立在把长度为 $N$ 的 DFT 变换到多维表示 $N={\prod }_l{N}_l$ 的基础上：

$$X[k]=\sum _{n=0}^{N-1}x[n]W_N^{nk}\text{\hspace{1em}(6-16)}$$

一般情况下只讨论两个因子的情形就足够了，因为更高维数可以通过反复迭代替换其中一个因子来实现。为了简化表示，下面只在二维索引变换内讨论 3 种 FFT 算法。

先把（时域）索引 $n$ 变换为：

$$n=A{n}_1+B{n}_2\mathrm{mod}N\{\begin{array}{l}0⩽n_1⩽N_1-1\\ 0⩽n_2⩽N_2-1\end{array}\text{\hspace{1em}(6-17)}$$

其中 $N=N_1{N}_2$，$A,B\in Z$ 是待定义的常数。利用这种索引变换，可以按下面的公式：

$$\left[x[0]x[1]x[2]\dots x[N-1]\right]=\left[\begin{array}{cccc}x[0,0]& x[0,1]& \dots & x[0,N_2-1]\\ x[1,0]& x[1,1]& \dots & x[1,N_2-1]\\ ⋮& ⋮& ⋮& ⋮\\ x[N_1-1,0]& x[N_1-1,1]& \dots & x[N_1-1,N_2-1]\end{array}\right]\text{\hspace{1em}(6-18)}$$

把一维数据映射成二维数组 $f:C^N\to C^{N_1\times N_2}$——一长串数据被”折”成了一个 $N_1\times N_2$ 的矩阵。再对输出（频）域的索引 $k$ 应用另一个映射：

$$k=C{k}_1+D{k}_2\mathrm{mod}N\{\begin{array}{l}0⩽k_1⩽N_1-1\\ 0⩽k_2⩽N_2-1\end{array}\text{\hspace{1em}(6-19)}$$

其中 $C,D\in Z$ 也是待定义的常数。由于 DFT 是双射，必须选择 $A,B,C,D$ 使得变换后的表示仍然唯一（保持双射）。Burrus${}^{[175]}$已经给出了针对具体 $N_1,N_2$ 选择这四个常数的一般方法（参阅练习 6.7 和 6.8）；本章给出的变换都是唯一的。

区分不同 FFT 算法的关键一点是：是否允许 $N_1$ 和 $N_2$ 有公因数，即 $\gcd (N_1,N_2)>1$（gcd 为最大公约数），还是二者必须互质。含 $\gcd (N_1,N_2)>1$ 的算法称为公因数算法（Common Factor Algorithm, CFA），而 $\gcd (N_1,N_2)=1$ 的称为质因数算法（Prime Factor Algorithm, PFA）。接下来讨论的 Cooley-Tukey FFT 属于 CFA 类型，而 Good-Thomas FFT 和 Winograd FFT 属于 PFA 类型。需要强调的是：Cooley-Tukey 算法的两个因子 $N=N_1{N}_2$ 彼此可以是互质的；反过来，PFA 只要求因子 $N_1,N_2$ 互质，它们本身并不一定是质数。例如长度 $N=12$ 的变换可分解为 $N_1=4$、$N_2=3$（互质），因此既可以用于 CFA FFT，也可以用于 PFA FFT。

6.3.1 Cooley-Tukey FFT 算法

Cooley-Tukey FFT 是所有 FFT 算法中最通用的，因为 $N$ 可以任意分解。最流行的形式是变换长度 $N$ 为基 $r$ 的幂，即 $N=r^v$，这类算法通常称为基 $r$ 算法。

Cooley 和 Tukey（更早是高斯）提出的索引映射也是最简单的。在式 (6-17) 中取 $A=N_2$、$B=1$，得到输入映射：

$$n=N_2{n}_1+n_2\{\begin{array}{l}0⩽n_1⩽N_1-1\\ 0⩽n_2⩽N_2-1\end{array}\text{\hspace{1em}(6-20)}$$

从 $n_1,n_2$ 的取值范围可以看出，式 (6-17) 中的取模运算根本不需要显式计算——两项之和天然落在 $[0,N-1]$ 内。对输出映射（式 6-19），Cooley 和 Tukey 取 $C=1$、$D=N_1$，得到：

$$k=k_1+N_1{k}_2\{\begin{array}{l}0⩽k_1⩽N_1-1\\ 0⩽k_2⩽N_2-1\end{array}\text{\hspace{1em}(6-21)}$$

这里同样可以省略取模计算。把式 (6-20) 和 (6-21) 代入 $W_N^{nk}$，得到：

$$W_N^{nk}=W_N^{N_2{n}_1{k}_1+N_1{N}_2{n}_1{k}_2+n_2{k}_1+N_1{n}_2{k}_2}\text{\hspace{1em}(6-22)}$$

由于 $W$ 的阶数为 $N=N_1{N}_2$，有 $W_N^{N_1}=W_{N_2}$ 和 $W_N^{N_2}=W_{N_1}$，于是式 (6-22) 化简为：

$$W_N^{nk}=W_{N_1}^{n_1{k}_1}{W}_N^{n_2{k}_1}{W}_{N_2}^{n_2{k}_2}\text{\hspace{1em}(6-23)}$$

注意关键结果：四项交叉项中 $N_1{N}_2{n}_1{k}_2$ 一项被完全消掉，剩下的三项分别只依赖 $\{n_1,k_1\}$、$\{n_2,k_1\}$、$\{n_2,k_2\}$。把式 (6-23) 代入式 (6-16) 的 DFT，求和就可以分解为内外两层：

$$\begin{array}{rl}& X[k_1,k_2]={\sum }_{n_2=0}^{N_2-1}{W}_{N_2}^{n_2{k}_2}\underset{\overline{x}[n_2,k_1]}{\underbrace{\left(W_N^{n_2{k}_1}\underset{N_1\text{点变换}}{\underbrace{\sum _{n_1=0}^{N_1-1}x[n_1,n_2]W_{N_1}^{n_1{k}_1}}}\right)}}\\ & =\underset{N_2\text{点变换}}{\underbrace{\sum _{n_2=0}^{N_2-1}{W}_{N_2}^{n_2{k}_2}\overline{x}[n_2,k_1]}}\end{array}\text{\hspace{1em}(6-24)}$$

(6-25)

这就是完整的 Cooley-Tukey 算法：先对所有 $n_2$ 行做 $N_1$ 点 DFT，乘上旋转因子 $W_N^{n_2{k}_1}$，再对所有 $k_1$ 列做 $N_2$ 点 DFT。原本 $N^2$ 量级的计算量被拆成两组小变换——这正是”快速”二字的来源。

算法 6.8 Cooley-Tukey 算法

前面我们已经知道，直接计算一个 $N$ 点 DFT 大约需要 $N^2$ 次复数乘法。当 $N$ 达到几百上千时，这样的计算量在实时系统里是难以接受的。Cooley-Tukey 算法的核心思想非常朴素：把一个大 DFT 拆成许多小 DFT。具体地，若 $N=N_1{N}_2$，一个 $N$ 点 DFT 就可以通过下列五个步骤完成：

(1) 根据式(6-20)计算输入序列的索引变换，把一维数据”折叠”成二维排布；

(2) 计算长度为 $N_1$ 的 $N_2$ 个小 DFT；

(3) 在第一个变换级的输出上应用旋转因子 $W_N^{n_2{k}_1}$；

(4) 计算长度为 $N_2$ 的 $N_1$ 个小 DFT；

(5) 根据式(6-21)计算输出序列的索引变换，把二维结果”展开”回一维。

为什么这样就快了？因为计算量从 $N^2$ 变成了大约 $N_1{N}_2(N_1+N_2)$ 的量级——两段小规模的代价之和，远小于一次大规模的代价。接下来用长度为 12 的变换把这个过程完整走一遍。

例 6.9 N = 12 的 Cooley-Tukey FFT

取 $N_1=4$、$N_2=3$，则索引映射为 $n=3n_1+n_2$ 和 $k=k_1+4k_2$。按这两个式子，把输入 $x[0]\sim x[11]$ 和输出 $X[0]\sim X[11]$ 分别排成表格（表 6-8、表 6-9）：

表 6-8 $N_1=4,N_2=3,n=3n_1+n_2$ 时的索引映射表

$n_2\backslash n_1$	0	1	2	3
0	x[0]	x[3]	x[6]	x[9]
1	x[1]	x[4]	x[7]	x[10]
2	x[2]	x[5]	x[8]	x[11]

表 6-9 $N_1=4,N_2=3,k=k_1+4k_2$ 时的索引映射表

$k_2\backslash k_1$	0	1	2	3
0	X[0]	X[1]	X[2]	X[3]
1	X[4]	X[5]	X[6]	X[7]
2	X[8]	X[9]	X[10]	X[11]

有了这两张表，就可以画出如图 6-12 所示的信号流程图：先用 4 点 DFT 计算 3 个变换，然后乘以旋转因子，最后计算 4 个 3 点 DFT。

图 6-12 N=12 的 Cooley-Tukey FFT

我们来算一算省了多少工作量。直接计算 12 点 DFT 需要 $12^2=144$ 次复数乘法和 $11^2=121$ 次复数加法。而 Cooley-Tukey FFT 中，旋转因子总共只需 12 次复数乘法（其中 8 次是 $\pm 1$ 或 $\pm j$ 这类”无关紧要”的乘法，不用真正去乘）。根据表 6-7，4 点 DFT 只需 8 次实数加法、不需要乘法；3 点 DFT 需要 4 次乘法和 6 次加法。若按算法 6.10 用 3 次加法、3 次乘法实现固定系数复数乘法，则 12 点 Cooley-Tukey FFT 的总工作量为：

$$\begin{array}{rl}& 3\times 16+4\times 3+4\times 12=108\text{ 次实数加法}\\ & 4\times 3+4\times 4=28\text{ 次实数乘法}\end{array}$$

而直接实现需要 $2\times 11^2+12^2\times 3=674$ 次实数加法和 $12^2\times 3=432$ 次实数乘法。加法减少到约 1/6，乘法减少到约 1/15——这就是”快速傅立叶变换”（Fast Fourier Transform, FFT）中”快速”二字的由来。

1. 基 r Cooley-Tukey 算法

Cooley-Tukey 算法区别于其他 FFT 算法的一个重要特点是：因子 $N_1$、$N_2$ 可以任意选取。这样就可以使用 $N=r^S$ 的基 r（radix-r）算法。最流行的是基 2 和基 4 算法，因为由表 6-7 可知，2 点、4 点 DFT 不需要任何乘法就能实现。

例如在 $S$ 级、$r=2$ 的情形下，索引映射为：

$$n=2^{S-1}{n}_1+\cdots +2n_{S-1}+n_S\text{\hspace{1em}(6-26)}$$ $$k=k_1+2k_2+\cdots +2^{S-1}{k}_S\text{\hspace{1em}(6-27)}$$

当 $S>2$ 时的一个一般惯例是：在信号流程图中，2 点 DFT 以蝶形图（butterfly）的形式表示。这种简化的画法基于两个约定：汇入同一节点的所有箭头相加，常系数乘法则标在箭头上。基 r 算法共有 ${\log }_r(N)$ 级，且每级使用同一类型的旋转因子。

图6-13 基2且长度为8的频域抽取算法

从图 6-13 的信号流程图可以看出一个对硬件实现极其重要的性质：计算可以**“就地”（in-place）完成**。也就是说，蝶形计算完就可以把结果写回原来两个输入所在的存储位置，因为下一步计算不再需要旧数据了——不需要额外的缓冲区。基 2 变换的旋转因子乘法总数是：

$${\log }_2(N)\cdot N/2\text{\hspace{1em}(6-28)}$$

原因是每两个箭头只共用一个旋转因子。

由于图 6-13 的算法是从频域一侧开始把原始 DFT 分成更短的 DFT（第一步就把输出分成奇偶两半），因此称为频域抽取（Decimation-In-Frequency, DIF）算法。它的典型特征是：输入按自然顺序出现，而输出频率值的索引是”位逆序”的。表 6-10 汇总了 DIF 基 2 算法各级的特征。

表 6-10 频率抽取的基 2 FFT

不同的指标	第1级	第2级	第3级	…	第 ${\log }_2(N)$ 级
组数	1	2	4	…	N/2
每组的蝶形数量	N/2	N/4	N/8	…	1
增量指数旋转因子	1	2	4	…	N/2

对称地，也可以从时域一侧拆分，构造时域抽取（Decimation In Time, DIT）算法。DIT 的特点是先把输入（时间）序列分开，结果所有频率值都按自然顺序出现（练习 6.10）。DIF 与 DIT 本质上是同一算法的转置，计算量完全相同。

那么位逆序到底是什么意思？图 6-14 给出了第 41 个索引在基 2 和基 4 算法下的索引变换。基 2 算法需要把二进制位的顺序整个翻转过来，这就是位逆序（bitreverse）；基 4 算法则先把每两位二进制组成一个”数字”，再逆转这些数字的顺序，称为数字逆序（digitreverse）。初学者只需记住：FFT 的每一级计算本身很规则，“不规则”全被推到了输入或输出的排序上，而排序在硬件里只是一次查表或计数器的翻转。

图 6-14 位逆序和数字逆序

2. 基 2 Cooley-Tukey 算法的实现

基 2 FFT 可以用蝶形处理器高效实现。一个蝶形处理器除了蝶形本身（一个复数加法器 + 一个复数减法器）之外，还包含旋转因子的复数乘法器。

普通的复数乘法 $(X+jY)\times (C+jS)$ 需要 4 次实数乘法和 2 次加/减法。但既然 $C$、$S$ 是旋转因子，可以预先计算并存在查找表中，我们不妨多存两个组合系数，把乘法次数降到 3 次——这就是下面的技巧。

算法 6.10 高效的复数乘法器

复数旋转因子乘法 $R+jI=(X+jY)\times (C+jS)$ 中，$C$ 和 $S$ 可预先算好存入表里。表中额外存储 3 个系数：

$$C、C+S\text{ 和 }C-S\text{\hspace{1em}(6-29)}$$

首先计算：

$$E=X-Y,Z=C\times E=C\times (X-Y)\text{\hspace{1em}(6-30)}$$

然后用：

$$R=(C-S)\times Y+Z\text{\hspace{1em}(6-31)}$$ $$I=(C+S)\times X-Z\text{\hspace{1em}(6-32)}$$

检验一下它确实等价于普通复数乘法：

$$\begin{array}{rl}R& =(C-S)Y+C(X-Y)=CY-SY+CX-CY=CX-SY\\ I& =(C+S)X-C(X-Y)=CX+SX-CX+CY=CY+SX\end{array}$$

正好是 $(X+jY)(C+jS)$ 的实部和虚部。这个算法共用了 3 次乘法、1 次加法和 2 次减法，代价是系数表要多存一列，并且关键路径上多了一级加法——在时序很紧的设计里，这一点需要在”省乘法”和”省延迟”之间权衡。

3. Cooley-Tukey 长度 256 FFT 算法的实现

掌握了蝶形和高效复数乘法器之后，就可以把所有部件拼成一个全尺寸 FFT。设计遵循表 6-10 的方案，为每一级维护蝶形数据索引、旋转因子增量、双节点间距和组大小的更新逻辑，最终得到图 6-15 所示的仿真结果：先把数据加载进 FFT 机器，然后逐级计算；结束时数据经过位逆序，并由 fft_valid 标志指示 FFT 结果可以送往输出端口。从仿真中还能直观验证表 6-10 的规律：旋转因子角 dw 的增量在每级以因子 2 递增，蝶形对索引间距 k2 每级减半。

图 6-15 256 点 FFT 模块的总体 MODELSIM 仿真

例 6.11 HDL 中的 256 点 FFT

原书 VHDL 改写为 Verilog：下面的可综合 Verilog 代码实现了 256 点 DIF（频域抽取）FFT。设计用一个状态机依次经历 load（载入数据）→ calc/update（计算蝶形并更新索引）→ reverse（位逆序输出）→ done 各状态；正余弦系数存放在两个 $128\times 16$ 位的查找表中（用 $readmemh 初始化）；复数乘法采用普通的 4 乘 2 加结构，以缩短最坏情况延迟路径。

//=============================================================
// 256 点 DIF FFT：寄存器阵列存放数据，ROM 存放旋转因子
//=============================================================
module fft256 #(
    parameter N    = 256,   // 变换点数
    parameter LdN  = 8,     // log2(N)
    parameter W    = 16     // 数据位宽
)(
    input  wire clk,        // 时钟
    input  wire reset,      // 异步复位
    input  wire signed [W-1:0] xr_in, xi_in, // 实部/虚部输入
    output reg         fft_valid,            // 输出有效标志
    output reg  signed [W-1:0] fftr, ffti,   // 实部/虚部输出
    output reg  [LdN:0] rcount_o,            // 位逆序索引计数器
    output wire signed [W-1:0] xr_out [0:7], // 前 8 个寄存器(调试)
    output wire signed [W-1:0] xi_out [0:7],
    output reg  [LdN:0] stage_o, gcount_o,   // 级计数/组计数
    output reg  [LdN:0] i1_o, i2_o,          // 蝶形双节点索引
    output reg  [LdN:0] k1_o, k2_o,          // 索引增量
    output reg  [LdN+1:0] w_o, dw_o,         // 旋转因子角度/增量
    output reg  [2:0]   wo                   // 状态机位置指示
);
 
    //---------- 状态定义 ----------
    localparam START  = 3'd0,
               LOAD   = 3'd1,
               CALC   = 3'd2,
               UPDATE = 3'd3,
               REVERSE= 3'd4,
               DONE   = 3'd5;
 
    reg [2:0] s;                              // 状态机变量
 
    //---------- 数据寄存器阵列 ----------
    reg signed [W-1:0] xr [0:N-1];
    reg signed [W-1:0] xi [0:N-1];
 
    //---------- 正余弦系数 ROM（Q14 定点，128 点一个象限）----
    reg signed [W-1:0] cos_rom [0:127];
    reg signed [W-1:0] sin_rom [0:127];
    initial $readmemh("cos_rom.txt", cos_rom); // 16384,16379,...,-16379
    initial $readmemh("sin_rom.txt", sin_rom); // 0,402,...,402
 
    reg [LdN-1:0] w;                          // 旋转因子角度
    reg signed [W-1:0] sin_r, cos_r;
 
    // 下降沿读 ROM，给出一个周期的读延迟
    always @(negedge clk) begin
        sin_r <= sin_rom[w];
        cos_r <= cos_rom[w];
    end
 
    //---------- 主状态机 ----------
    integer count, i1, i2, gcount, k1, k2;
    integer stage, dw, rcount, k;
    reg signed [2*W-1:0] tr, ti, pr, pi;
    reg [LdN-1:0] slv, rslv;
 
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            s <= START;
        end else begin
            case (s)
            //------------------------------------------------
            START: begin                       // 初始化
                count  <= 0;  gcount <= 0;  stage <= 1;
                i1     <= 0;  i2 <= N/2;   k1 <= N;  k2 <= N/2;
                dw     <= 1;  fft_valid <= 0;  w <= 0;
                s      <= LOAD;
            end
            //------------------------------------------------
            LOAD: begin                        // 读入 256 点数据
                xr[count] <= xr_in;
                xi[count] <= xi_in;
                count <= count + 1;
                s     <= (count == N-1) ? CALC : LOAD;
            end
            //------------------------------------------------
            CALC: begin                        // 计算一个蝶形
                tr = xr[i1] - xr[i2];          // 下支路差
                ti = xi[i1] - xi[i2];
                xr[i1] <= xr[i1] + xr[i2];     // 上支路和(就地写回)
                xi[i1] <= xi[i1] + xi[i2];
                pr = cos_r*tr + sin_r*ti;      // W = cos - j*sin
                pi = cos_r*ti - sin_r*tr;
                xr[i2] <= pr >>> 14;           // Q14 定点缩放
                xi[i2] <= pi >>> 14;
                s <= UPDATE;
            end
            //------------------------------------------------
            UPDATE: begin                      // 更新各层计数器
                s <= CALC;                     // 默认继续下一个蝶形
                i1 <= i1 + k1;  i2 <= i1 + k1 + k2;
                wo <= 1;
                if (i1 + k1 >= N-1) begin      // 本组蝶形算完？
                    gcount <= gcount + 1;
                    i1 <= gcount + 1;  i2 <= gcount + 1 + k2;
                    wo <= 2;
                    if (gcount + 1 >= k2) begin // 本级所有组算完？
                        gcount <= 0;  i1 <= 0;  i2 <= k2/2;
                        dw <= dw * 2;           // 旋转因子增量加倍
                        stage <= stage + 1;
                        wo <= 3;
                        if (stage + 1 > LdN) begin // 全部级算完
                            s <= REVERSE;  count <= 0;  wo <= 4;
                        end else begin             // 开始新的一级
                            k1 <= k2;  k2 <= k2/2;
                            i1 <= 0;  i2 <= k2/2;
                            w <= 0;  wo <= 5;
                        end
                    end else begin              // 开始新的一组
                        i1 <= gcount + 1;  i2 <= gcount + 1 + k2;
                        w <= w + dw;  wo <= 6;
                    end
                end
            end
            //------------------------------------------------
            REVERSE: begin                     // 位逆序并输出
                fft_valid <= 1;
                slv  = count[LdN-1:0];
                for (k = 0; k < LdN; k = k + 1)
                    rslv[k] = slv[LdN-1-k];    // 位序翻转
                rcount = rslv;
                fftr <= xr[rcount];
                ffti <= xi[rcount];
                count <= count + 1;
                s <= (count == N-1) ? DONE : REVERSE;
            end
            //------------------------------------------------
            DONE: s <= START;                  // 回到起点做下一帧
            default: s <= START;
            endcase
        end
    end
 
    //---------- 调试信号引出 ----------
    always @(*) begin
        i1_o = i1;  i2_o = i2;  stage_o = stage; gcount_o = gcount;
        k1_o = k1;  k2_o = k2;  w_o = {2'b0, w}; dw_o = dw;
        rcount_o = rcount;
    end
    genvar g;
    generate for (g = 0; g < 8; g = g + 1) begin : dbg
        assign xr_out[g] = xr[g];
        assign xi_out[g] = xi[g];
    end endgenerate
 
endmodule

这段代码的结构值得仔细读一遍。状态机从 START 状态初始化主要变量和循环计数器；LOAD 状态把 256 个实部、虚部数据依次读入内部的 $256\times 16$ 寄存器文件；随后 CALC 与 UPDATE 两个状态交替运行多个时钟周期——CALC 中完成一次蝶形（注意这里刻意使用常规的 4 乘 2 加复数乘法，以缩短最坏情况延迟路径），UPDATE 中更新循环计数器和数据索引。可以看到三个嵌套的循环层次：蝶形 → 组 → 级，正好对应表 6-10 中的三行。所有级算完后进入 REVERSE 状态，用位翻转计数器逐个把 FFT 值送到输出端口，同时把 fft_valid 置 1，通知外部”输出端口上现在是有效数据”。

原设计在 FPGA 上的资源消耗为 34340 个 LE、8 个嵌入式乘法器，用 TimeQuest 缓慢 85C 模型测得时序性能 Fmax = 31.12 MHz。一个有意思的细节是优化目标的影响：若以面积（Area）为目标综合，正弦、余弦 LUT 会被合并进嵌入式 M9K 存储块，LE 数量下降；若以速度（Speed）为目标，LUT 会被综合进 LE。而在 Verilog 版本中使用 $readmemh() 初始化 ROM，则无论速度还是面积优化，每次都有两个嵌入式 M9K 被综合出来。

仿真验证

图 6-16 显示了用 MODELSIM 对三角形输入序列 $x[n]=\{20,40,60,80,100,120,140,160,\dots \}$（只有前八个值非零）的仿真起始波形。测试序列可以在 MATLAB 中生成。

图 6-16 FFT 模拟输入数据。输入帧的开始

用以下指令量化并列出前 10 个样本，与 MODELSIM 仿真结果对照：

sprintf('%d', real(round(Y(1:9))))
sprintf('%d', imag(round(Y(1:9))))

得到的（预期）测试数据为：

720 714 698 671 634 587 532 471 403 ... (实数部分)
0 -82 -163 -240 -313 -380 -439 -490 -532 ... (虚数部分)

最后，在 $57.7\mu s$ 处，从图 6-17 可以看到 DC 值即 $\sum x_i=720$——这正是 DFT 在 $k=0$ 处等于输入之和的性质的直接体现。除了一些量化误差外，仿真结果与 MATLAB 测试数据的其他 FFT 值都吻合。

图 6-17 FFT 核心模拟输出结果。输出帧的开始

还可以怎样改进？在 LE 资源方面最有效的一步是用嵌入式存储器块替换三个大型寄存器文件——但代价是 M9K 只能作同步存储器使用，FFT 机器里要增加额外的等待状态。

6.3.2 Good-Thomas FFT 算法

Good 和 Thomas 提出的索引变换能把长度为 $N=N_1{N}_2$ 的 DFT 变成”实际的”二维 DFT——也就是说，没有旋转因子。天下没有免费的午餐：它要求各因子互质（即 $\gcd (N_k,N_l)=1,k\ne l$），而且索引映射比 Cooley-Tukey 复杂，在线计算索引（而不是查预存的表）时开销明显。

回忆 Cooley-Tukey 的做法：通过 $n$ 和 $k$ 的索引映射，在 $W_N^{nk}$ 中引入了交叉项旋转因子。Good-Thomas 的目标是把交叉项彻底消掉。把 $n=A{n}_1+B{n}_2$、$k=C{k}_1+D{k}_2$ 代入 $W_N^{nk}$：

$$\begin{array}{rl}{W}_N^{nk}& =W_N^{(A{n}_1+B{n}_2)(C{k}_1+D{k}_2)}\\ & =W_N^{AC{n}_1{k}_1+AD{n}_1{k}_2+BC{k}_1{n}_2+BD{n}_2{k}_2}\\ & =W_N^{N_2{n}_1{k}_1}{W}_N^{N_1{k}_2{n}_2}=W_{N_1}^{n_1{k}_1}{W}_{N_2}^{k_2{n}_2}\end{array}\text{\hspace{1em}(6-33)}$$

要让交叉项 $AD{n}_1{k}_2$ 和 $BC{k}_1{n}_2$ 消失、剩余项恰好变成小 DFT 的旋转因子，必须同时满足：

$$⟨AD{⟩}_N=⟨BC{⟩}_N=0\text{\hspace{1em}(6-34)}$$ $$⟨AC{⟩}_N=N_2\text{\hspace{1em}(6-35)}$$ $$⟨BD{⟩}_N=N_1\text{\hspace{1em}(6-36)}$$

Good 和 Thomas 给出的一组解是：

$$A=N_{2}B=N_{1}C=N_2⟨N_2^{-1}{⟩}_{N_1}D=N_1⟨N_1^{-1}{⟩}_{N_2}\text{\hspace{1em}(6-37)}$$

检验：$AD$ 与 $BC$ 都含有因子 $N_1{N}_2=N$，所以模 $N$ 后为 0，式(6-34)成立。再由 $\gcd (N_1,N_2)=1$ 和欧拉定理，逆元可写成 $N_2^{-1}\mathrm{mod}{N}_1=N_2^{\varphi (N_1)-1}\mathrm{mod}{N}_1$（$\varphi$ 是欧拉函数）。于是式(6-35)可以改写为：

$$⟨AC{⟩}_N={⟨N_2{N}_2{⟨N_2^{\varphi (N_1)-1}⟩}_{N_1}⟩}_N\text{\hspace{1em}(6-38)}$$

把内部的模化简展开（$v\in Z$，且 $v{N}_1{N}_2\mathrm{mod}N=0$），得到最终形式：

$$⟨AC{⟩}_N={⟨N_2{N}_2(N_2^{\varphi (N_1)-1}+v{N}_1)⟩}_N=N_2\text{\hspace{1em}(6-39)}$$

同样的证明也适用于式(6-36)，因此采用 Good-Thomas 映射后，式(6-34)~式(6-36)三个条件全部满足。于是有如下定理。

定理 6.12 Good-Thomas 索引映射

对 $n$ 的索引映射是：

$$n=(N_2{n}_1+N_1{n}_2)\mathrm{mod}N,\{\begin{array}{l}0⩽n_1⩽N_1-1\\ 0⩽n_2⩽N_2-1\end{array}\text{\hspace{1em}(6-40)}$$

对 $k$ 的索引映射是：

$$k=(N_2⟨N_2^{-1}{⟩}_{N_1}{k}_1+N_1⟨N_1^{-1}{⟩}_{N_2}{k}_2)\mathrm{mod}N,\{\begin{array}{l}0⩽k_1⩽N_1-1\\ 0⩽k_2⩽N_2-1\end{array}\text{\hspace{1em}(6-41)}$$

式(6-41)与定理(2-13)的中国余数定理（CRT）形式完全一致，所以 $k_1$、$k_2$ 可以简单地用模化简得到：$k_l=k\mathrm{mod}{N}_l$。这也是”输出映射复杂、但反推容易”的来源。

把 Good-Thomas 索引映射代入 DFT 矩阵（式(6-16)），就得到：

$$X[k_1,k_2]=\sum _{n_2=0}^{N_2-1}{W}_{N_2}^{n_2{k}_2}\underset{N_1\text{点变换}:\overline{x}[n_2,k_1]}{\underbrace{\left(\sum _{n_1=0}^{N_1-1}x[n_1,n_2]W_{N_1}^{n_1{k}_1}\right)}}=\underset{N_2\text{点变换}}{\underbrace{\sum _{n_2=0}^{N_2-1}{W}_{N_2}^{n_2{k}_2}\overline{x}[n_2,k_1]}}\text{\hspace{1em}(6-42)}$$

这就是”真正的”二维 DFT：先做 $N_1$ 点变换，再做 $N_2$ 点变换，两级之间不需要乘任何旋转因子。完整算法如下，形式上与 Cooley-Tukey 算法 6.8 相似，区别只在索引映射和没有旋转因子一步。

算法 6.13 Good-Thomas FFT 算法

$N=N_1{N}_2$ 点的 DFT 可按下列步骤计算：

(1) 根据式(6-40)进行输入序列的索引变换。

(2) 计算长度为 $N_1$ 的 $N_2$ 次 DFT。

(3) 计算长度为 $N_2$ 的 $N_1$ 次 DFT。

(4) 根据式(6-41)进行输出序列的索引变换。

例 6.14 N = 12 的 Good-Thomas FFT

取 $N_1=4$、$N_2=3$（互质，满足前提）。输入索引按 $n=3n_1+4n_2\mathrm{mod}12$ 映射，输出索引按 $k=9k_1+4k_2\mathrm{mod}12$ 映射（$9=3⟨3^{-1}{⟩}_4$，$4=4⟨4^{-1}{⟩}_3$），得到表 6-11 和表 6-12：

表 6-11 输入索引的映射

$n_2\backslash n_1$	0	1	2	3
0	x[0]	x[3]	x[6]	x[9]
1	x[4]	x[7]	x[10]	x[1]
2	x[8]	x[11]	x[2]	x[5]

表 6-12 输出索引结果的映射

$k_2\backslash k_1$	0	1	2	3
0	X[0]	X[9]	X[6]	X[3]
1	X[4]	X[1]	X[10]	X[7]
2	X[8]	X[5]	X[2]	X[11]

对照例 6.9 的表 6-8 可以看到：输入映射 $n=3n_1+4n_2$ 与 $n=3n_1+n_2$ 在第一行相同，但后续行元素开始”绕回”（如 $x[1]$ 出现在 $n_1=3,n_2=1$ 处）——这就是 $\mathrm{mod}12$ 的效果。利用这些索引变换可以构造图 6-18 的信号流程图：第一级是 3 个 4 点 DFT，第二级是 4 个 3 点 DFT，各级之间不需要乘旋转因子，这正是 Good-Thomas 的卖点。

图6-18 $N=12$ 的 Good-Thomas FFT

6.3.3 Winograd FFT 算法

Winograd FFT 算法建立在对 $N_1\times N_2$ 维逆 DFT 矩阵（不含前因子 $N^{-1}$）的观察之上，要求 $\gcd (N_1,N_2)=1$：

$$x[n]=\sum _{k=0}^{N-1}X[k]W_N^{-nk}\text{\hspace{1em}(6-44)}$$ $$\mathbf{x}={\mathbf{W}}_N^{\ast }\mathbf{X}\text{\hspace{1em}(6-45)}$$

关键工具是 Kronecker 乘积：当 $N_1$、$N_2$ 互质时，上述矩阵可以写成两个分别为 $N_1$ 维和 $N_2$ 维的二次 IDFT 矩阵的 Kronecker 乘积。与 Good-Thomas 一样，$X[k]$ 和 $x[n]$ 的索引先按二维模式排列，再逐行读出。看一个 $N=12$ 的例子。

例 6.15 使用 Kronecker 乘积且 N=12 的 IDFT

取 $N_1=4$、$N_2=3$，输出索引按 Good-Thomas 映射 $k=9k_1+4k_2\mathrm{mod}12$ 排列。其中 Kronecker 乘积的定义为 $\mathbf{A}\otimes \mathbf{B}=[a[i,j]]\mathbf{B}$，即用 $a[i,j]$ 逐个缩放整个矩阵 $\mathbf{B}$ 拼成分块大矩阵。排列过程如下：

$$\left[\begin{array}{c}X[0]\\ X[1]\\ X[2]\\ X[3]\\ X[4]\\ X[5]\\ X[6]\\ X[7]\\ X[8]\\ X[9]\\ X[10]\\ X[11]\end{array}\right]\to \begin{array}{ccccc}{k}_2\backslash k_1& 0& 1& 2& 3\\ 0& X[0]& X[9]& X[6]& X[3]\\ 1& X[4]& X[1]& X[10]& X[7]\\ 2& X[8]& X[5]& X[2]& X[11]\end{array}\to \left[\begin{array}{c}X[0]\\ X[9]\\ X[6]\\ X[3]\\ X[4]\\ X[1]\\ X[10]\\ X[7]\\ X[8]\\ X[5]\\ X[2]\\ X[11]\end{array}\right]$$

据此构造长度为 12 的 IDFT 矩阵：

$$\left[\begin{array}{c}x[0]\\ x[9]\\ x[6]\\ x[3]\\ x[4]\\ x[1]\\ x[10]\\ x[7]\\ x[8]\\ x[5]\\ x[2]\\ x[11]\end{array}\right]=\left(\left[\begin{array}{ccc}{W}_{12}^0& W_{12}^0& W_{12}^0\\ W_{12}^0& W_{12}^{-4}& W_{12}^{-8}\\ W_{12}^0& W_{12}^{-8}& W_{12}^{-4}\end{array}\right]\otimes \left[\begin{array}{cccc}{W}_{12}^0& W_{12}^0& W_{12}^0& W_{12}^0\\ W_{12}^0& W_{12}^{-3}& W_{12}^{-6}& W_{12}^{-9}\\ W_{12}^0& W_{12}^{-6}& W_{12}^0& W_{12}^{-6}\\ W_{12}^0& W_{12}^{-9}& W_{12}^{-6}& W_{12}^{-3}\end{array}\right]\right)\left[\begin{array}{c}X[0]\\ X[9]\\ X[6]\\ X[3]\\ X[4]\\ X[1]\\ X[10]\\ X[7]\\ X[8]\\ X[5]\\ X[2]\\ X[11]\end{array}\right]$$

用速记符号 $\tilde{x}$ 表示重排后的序列 $x$，上式就是矩阵/向量形式：

$$\tilde{\mathbf{x}}=({\mathbf{W}}_{N_1}\otimes {\mathbf{W}}_{N_2})\tilde{\mathbf{X}}\text{\hspace{1em}(6-46)}$$

现在轮到 Winograd DFT 算法（算法 6.7）登场了。对每个短 DFT，有：

$${\mathbf{W}}_{N_l}={\mathbf{C}}_l\times {\mathbf{B}}_l\times {\mathbf{A}}_l\text{\hspace{1em}(6-47)}$$

其中 ${\mathbf{A}}_l$ 实现输入加法，${\mathbf{B}}_l$ 是傅立叶系数构成的对角矩阵，${\mathbf{C}}_l$ 实现输出加法。把式(6-47)代入式(6-46)，并利用”矩阵乘法与 Kronecker 乘积可交换结合顺序”的性质，得到：

$${\mathbf{W}}_{N_1}\otimes {\mathbf{W}}_{N_2}=({\mathbf{C}}_1{\mathbf{B}}_1{\mathbf{A}}_1)\otimes ({\mathbf{C}}_2{\mathbf{B}}_2{\mathbf{A}}_2)=({\mathbf{C}}_1\otimes {\mathbf{C}}_2)({\mathbf{B}}_1\otimes {\mathbf{B}}_2)({\mathbf{A}}_1\otimes {\mathbf{A}}_2)\text{\hspace{1em}(6-48)}$$

这一步是整个 Winograd FFT 的灵魂，值得停下来想想它为什么省乘法：

• ${\mathbf{A}}_l$、${\mathbf{C}}_l$ 都是简单加法矩阵，它们的 Kronecker 乘积 ${\mathbf{A}}_1\otimes {\mathbf{A}}_2$、${\mathbf{C}}_1\otimes {\mathbf{C}}_2$ 仍然只是加法矩阵——加法不需要乘法器；
• 两个对角矩阵的 Kronecker 乘积仍是对角矩阵。若按表 6-7，两个小 Winograd DFT 各需 $M_1$、$M_2$ 次乘法，那么总乘法次数就是对角阵 ${\mathbf{B}}_1\otimes {\mathbf{B}}_2$ 的对角元个数 $M_1{M}_2$——乘法次数相乘，而不是相加。

把上述步骤组装起来，就得到 Winograd FFT 的构造方法和执行流程。

定理 6.16 Winograd FFT 的设计

$N=N_1{N}_2$ 点且 $N_1$、$N_2$ 互质的变换按下列步骤构造：

(1) 根据式(6-40)中的 Good-Thomas 映射（按索引行读取）对输入序列进行索引变换；

(2) 利用 Kronecker 乘积对 DFT 矩阵进行因数分解；

(3) 用 Winograd DFT 算法替换长度为 $N_1$ 和 $N_2$ 的 DFT 矩阵；

(4) 集中乘法。

定理 6.17 Winograd FFT 算法

构造完成后，计算只需 3 步：

(1) 计算前置加法 ${\mathbf{A}}_1$ 和 ${\mathbf{A}}_2$；

(2) 根据矩阵 ${\mathbf{B}}_1\otimes {\mathbf{B}}_2$ 计算 $M_1{M}_2$ 次乘法；

(3) 根据 ${\mathbf{C}}_1$ 和 ${\mathbf{C}}_2$ 计算后置加法。

例 6.18 长度为 12 的 Winograd FFT

按定理 6.16 构造 $N_1=3$、$N_2=4$ 所需的矩阵：

$${\mathbf{A}}_1\otimes {\mathbf{A}}_2=\left[\begin{array}{ccc}1& 1& 1\\ 0& 1& 1\\ 0& 1& -1\end{array}\right]\otimes \left[\begin{array}{cccc}1& 1& 1& 1\\ 1& -1& 1& -1\\ 1& 0& -1& 0\\ 0& 1& 0& -1\end{array}\right]\text{\hspace{1em}(6-49)}$$ $${\mathbf{B}}_1\otimes {\mathbf{B}}_2=\mathrm{diag}(1,-3/2,\sqrt{3}/2)\otimes \mathrm{diag}(1,1,1,-i)\text{\hspace{1em}(6-50)}$$ $${\mathbf{C}}_1\otimes {\mathbf{C}}_2=\left[\begin{array}{ccc}1& 0& 0\\ 1& 1& i\\ 1& 1& -i\end{array}\right]\otimes \left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 1& 0& -1& 0\\ 0& 0& 1& -1\end{array}\right]\text{\hspace{1em}(6-51)}$$

三个矩阵中，输入加法阵和输出加法阵都不需要乘法器；对角阵 ${\mathbf{B}}_1\otimes {\mathbf{B}}_2$ 共有 $3\times 4=12$ 个对角元，但其中系数 $1$ 的乘法可以省去，实数乘法总数为 $2\times 3\times 4=24$。对比直接算法的 432 次实数乘法，节省是惊人的。

到目前为止我们用 Winograd FFT 算的是 IDFT。要借助 IDFT 计算 DFT，可以复用式(6-6)的共轭技巧，用矩阵/向量表示即：

$$x^{\ast }=(W_N^{\ast }X{)}^{\ast }\text{\hspace{1em}(6-52)}$$ $$x^{\ast }=W_N{X}^{\ast }\text{\hspace{1em}(6-53)}$$

也就是说，若 $W_N=[e^{2\pi jnk/N}]$ 是 DFT 矩阵，则计算 DFT 的流程为：先取输入序列的共轭复数，再用 IDFT 算法变换它，最后取输出序列的共轭复数——三次”取共轭”换一次 DFT。

也可以用 Kronecker 乘积公式（即 Winograd FFT 本身）直接计算 DFT，只是输出索引映射会略有不同，如下例。

例 6.19 用 Kronecker 乘积公式计算 12 点 DFT

$$\left[\begin{array}{c}X[0]\\ X[3]\\ X[6]\\ X[9]\\ X[4]\\ X[7]\\ X[10]\\ X[1]\\ X[8]\\ X[11]\\ X[2]\\ X[5]\end{array}\right]=\left(\left[\begin{array}{ccc}{W}_{12}^0& W_{12}^0& W_{12}^0\\ W_{12}^0& W_{12}^4& W_{12}^8\\ W_{12}^0& W_{12}^8& W_{12}^4\end{array}\right]\otimes \left[\begin{array}{cccc}{W}_{12}^0& W_{12}^0& W_{12}^0& W_{12}^0\\ W_{12}^0& W_{12}^3& W_{12}^6& W_{12}^9\\ W_{12}^0& W_{12}^6& W_{12}^0& W_{12}^6\\ W_{12}^0& W_{12}^9& W_{12}^6& W_{12}^3\end{array}\right]\right)\left[\begin{array}{c}x[0]\\ x[9]\\ x[6]\\ x[3]\\ x[4]\\ x[1]\\ x[10]\\ x[7]\\ x[8]\\ x[5]\\ x[2]\\ x[11]\end{array}\right]\text{\hspace{1em}(6-54)}$$

输入序列 $x[n]$ 仍按 Good-Thomas 映射的顺序排列；对比例 6.15 的输出排列可以发现，频率输出索引映射中第一个和第三个元素交换了位置（$X[0],X[3],X[6],X[9]$ 出现在开头而非 $X[0],X[9],X[6],X[3]$）。使用时只需记住这一平滑修正即可。

6.3.4 DFT 和 FFT 算法的比较

到目前为止，实现 DFT 的途径已经相当多了：先选定一种短 DFT 算法（直接法或 Rader/Winograd 短 DFT），再通过 Cooley-Tukey、Good-Thomas 或 Winograd 的索引映射组装成长 DFT。选择时的共同目标通常是把乘法复杂度降到最低——这是一个合理的准则，因为乘法的实现成本（面积、功耗、延迟）远高于加法、数据访问或索引计算。

图 6-19 给出了各种 FFT 长度所需乘法的次数。单从乘法复杂度看，Winograd FFT 最有吸引力。本章已经用 $N=4\times 3=12$ 点 FFT 的几种设计做了示范，表 6-13 汇总了直接算法、Rader 质因子算法（RPFA）以及三种索引映射（Good-Thomas、Cooley-Tukey、Winograd）下的 Winograd FFT 的乘法数量。

图 6-19 基于所需的实数乘法次数的不同 FFT 算法的比较

表 6-13 长度为 12 的复数输入 FFT 算法的实数乘法数量（不考虑旋转因子乘以 $W^0$，假设一次复数乘法使用 4 次实数乘法）

DFT 方法	Good-Thomas（图 6-18）	Cooley-Tukey（图 6-2）	Winograd（例 6.15）
直接算法	$4\times 12^2=576$	$4\times 12^2=576$	$4\times 12^2=576$
RPFA	$4(3(4-1{)}^2+4(3-1{)}^2)=172$	$4(43+6)=196$	—
WFTA	$3\times 0\times 2+4\times 2\times 2=16$	$16+4\times 6=40$	$2\times 3\times 4=24$

但乘法次数并不是唯一标准，还要考虑：可实现的变换长度、加法次数、索引计算的系统开销、系数与数据存储器大小、运行时代码长度等。综合权衡后，Cooley-Tukey 方法往往提供最佳总体解决方案，见表 6-14。

表 6-14 长度 $N=\prod N_k$ 的 FFT 算法的重要性质

性质	Cooley-Tukey	Good-Thomas	Winograd
任意变换长度	是	否，$\gcd (N_k,N_l)=1$	否，$\gcd (N_k,N_l)=1$
$W$ 的最大阶数	$N$	$\max (N_k)$	$\max (N_k)$
是否需要旋转因子	是	否	否
乘法	差	中	最佳
加法	中	中	中
索引计算量	最佳	中	差
原地数据	是	是	否
实现的优点	小规模蝶形处理器	可使用 RPFA，快速、简单的 FIR 阵列	小规模的完全并行、中等规模 FFT（<50）

怎么理解这张表？Cooley-Tukey 乘法最多，但索引最简单、支持任意长度、支持就地计算，工程上最省心；Winograd 乘法最少，但索引计算复杂、不能原地计算、长度受限（因子必须互质），适合小规模完全并行结构；Good-Thomas 介于两者之间，免去旋转因子又保留原地计算。

表 6-15 一些 FPGA FFT 实现的比较

名称	数据类型	FFT 类型	N 点 FFT 的时间	时钟速率	内部 RAM/ROM	设计目标/源
Xilinx FPGA	8 位	基 2 FFT	N=256, 102.4 μs	70 MHz, 4.8 W@3.3V	否	573 个 CLB [184]
Xilinx FPGA	16 位	AFT	N=256, 82.48/42.08 μs	50 MHz, 15.6 W@3.3V / 29.5 W	否	2602/4922 个 CLB [200]
Xilinx FPGA ERNS-NTT	12.7 位	使用 FFT 的 NTT	N=97, 9.24 μs	26 MHz, 3.5 W@3.3V	否	1178 个 CLB [201]

表中 Goslin 的设计基于基 2 FFT，蝶形部分采用第 2 章讨论的分布式算法实现；Dandalis 等人的设计基于所谓算术傅立叶变换（AFT）的 DFT 近似；Meyer-Bäse 等人的 ERNS FFT 则把 Rader 算法与数论变换（NTT）结合了起来。

当 FPGA 的规模超过 1M 门以后，FFT 完全可以集成在单片 FPGA 上。不过由于 FFT 模块的设计是劳动密集型工作，实践中通常直接购买商用”知识产权”（Intellectual Property, IP）模块（有时也称虚拟元件 Virtual Component, VC）更划算，例如访问 www.xilinx.com 或 www.altera.com 参阅其 IP 合作程序。多数商用设计基于基 2 或基 4 FFT。

6.3.5 IP 内核 FFT 设计

Altera 和 Xilinx 都提供 FFT 生成器——这是除 FIR 滤波器之外最重要的 IP 模块类别（IP 模块的介绍见 1.4.4 节）。Xilinx 虽然提供一些免费的固定长度、固定位宽的硬内核，但通常从 FPGA 供应商处采购参数化 FFT 内核需要支付一笔（合理的）许可费用。

我们再以例 6.11 讨论过的 256 点 FFT 为例，看看用 Altera FFT 编译器生成同样功能的设计是什么样——这次连控制进程的 FSM 都包含在内。Altera 优化的 FFT 宏函数内核是高性能、高度参数化的 FFT 处理器，专用支持 Stratix 和 Cyclone II 器件（不支持已停产的 APEX 和 Flex 器件）。它的主要特性包括：

• 实现基 2/4 频域抽取（DIF）FFT 算法，变换长度为 $2^S$，其中 $6⩽S⩽16$（即 64 点到 65536 点）；
• 内部的浮点体系结构模块用于在变换计算中最大化信号的动态范围——定点 FFT 每级都可能引入位增长，块浮点是一种折中方案；
• 可用 IP 工具 MegaWizard 设计环境设定各种 FFT 体系结构，包括 $4\times 2+$ 和 $3\times 5+$ 蝶形体系结构以及不同的并行体系结构；
• 内含旋转因子的系数生成器，旋转因子存储在 M9K 存储模块中。

对初学者而言，这一节的启示是：理解了前面手写的 256 点 FFT 状态机，你就具备了评估和使用这类商用 IP 内核的能力——它们做的正是同一件事，只是在长度、位宽、动态范围和流水线结构上做了更全面的参数化封装。

例 6.20 用 IP 核生成长度为 256 的 FFT

前面几节我们一直自己动手搭 FFT 蝶形运算单元，这一节换个思路：让 FPGA 厂商的工具直接”生成”一个现成的 FFT。这种现成的、参数可配置的电路模块叫做 IP 核（Intellectual Property Core）。它就像集成电路设计里的”标准件”——旋转因子表、流水线调度、存储器组织这些繁琐又容易出错的细节都由工具替你做好，你只需要回答几个问题：“变换长度是多少？数据多少位？用什么架构？”

生成流程是这样的：在 Quartus II 的 Tools 菜单下打开 MegaWizard Plug-In Manager（向导式生成器），在库目录 DSP | Transform 下找到 FFT 生成器。先给内核起一个设计名称，然后进入 ToolBench 配置界面，如图 6-20(a) 所示。

(a) IP 工具测试台

(b) 系数的指定

图 6-20 FFT 的 IP 设计

配置分几步走：

1. 在 Parameters 模块中把 FFT 长度设为 256，数据与系数精度都设为 16 位。
2. 选择体系结构。向导提供三种选择：Streaming（流式）、Buffered Burst（缓冲突发） 和 Burst（突发）。区别在于内部缓冲器的多少：Streaming 完全流水化，每 256 个时钟周期送入一组新数据，处理后连续吐出一组新结果，中间不需要额外的等待周期；另外两种则要靠额外的缓冲器逐块处理，吞吐率依次下降，但占用的存储器资源更少。如果不做会怎样？如果你选了省资源的 Burst 模式却按流式节奏喂数据，数据就会丢失或错位——架构选择必须与系统数据流速匹配。
3. 在 Implementation Options 选项卡中选择 “3×5+” 乘法器结构。

以 Cyclone II 器件为目标时，三种架构的逻辑资源估算如表 6-16 所示。注意看最后一行”模块吞吐量周期”：Streaming 只要 256 个周期就能完成一个模块的吞吐，而 Burst 需要 775 个周期，差了 3 倍。

表 6-16 Cyclone II 系列的逻辑资源估算

资源	Streaming	Buffered Burst	Burst
LE	4581	4638	4318
M9K	11	9	5
M114K RAM	0	0	0
MLAB	0	0	0
DSP模块9位	18	18	18
变换计算周期	256	258	262
模块吞吐量周期	256	331	775

工具测试平台的第 2 步会生成 ModelSim 仿真所需的仿真模型；第 3 步生成 HDL 代码和所有辅助文件。几秒钟之内，旋转因子系数文件、MATLAB 测试平台和 ModelTech 仿真脚本就全部就绪了。表 6-17 列出了这些文件——值得注意的是，工具不仅给出综合用的设计文件，还同时给出 MATLAB（位精度）和 ModelSim（周期精度）两套测试向量，这给验证工作省了大量时间。

表 6-17 为 FFT 内核生成的 IP 文件

文件	说明
fft256.vhd	定义定制的宏函数的顶层说明的宏函数变量文件
fft256.	例化文件
fft256.cmp	MegaCore 函数变量的元件声明
fft256.bsf	用在 Quartus II 模块图编辑器中的符号文件
fft256.vho	ModelSim 仿真所用的功能模型
fft_tb.vhd	ModelSim 仿真所用的测试平台
fft256_model.m	为定制的 FFT 提供 MATLAB 仿真模型
fft256_tb.m	为定制的 FFT 提供 MATLAB 测试平台
*.txt	具备随机实部和虚部输入数据的两个文本文件
*.hex	6 个 sin/cos 旋转因子表
fft256_nativelink.tcl	用于设置 Quartus II 软件中 NativeLink 的一个 Tcl 脚本
fft256.qip	包含 MegaCore 函数变量的 Quartus II 项目信息
fft256_syn.v	用于某些第三方综合工具的时间和资源估计网表
fft256.html	MegaCore 函数报告文件

为了与 DIF 示例 6.11 的仿真结果对照，我们要换成自己的测试数据：把工具生成的 fft256_real_input.txt 和 fft256_imag_input.txt 替换掉。测试序列用一个短的三角形序列，虚部置零。在 MATLAB 中这样生成数据文件（共 1024 个样本，多出来的部分补零）：

x=[(1:8)*20,zeros(1,248+3*256)]; % Total 1024 samples
fid=fopen('fft256_real_umb.txt','w');
fprintf(fid,'%d\r\n',x);fclose(fid);
Y=fft(x);

将零向量作为虚数存储。FFT 内核内部采用模块浮点（块浮点）格式，输出自带一个块指数，这里对应缩放 $2^{-3}$。用下面的指令可以在 ModelSim 仿真中量化并列出缩放后的前 10 个样本，作为期望值：

sprintf('%d', real(round(Y(1:10)*2^-3)))
sprintf('%d', imag(round(Y(1:10)*2^-3)))

（期望的）测试数据如下：

90 89 87 84 79 73 67 59 50 41... (实部)
0 -10 -20 -30 -39 -47 -55 -61 -66 -70... (虚部)

工具生成的是 VHDL 例化封装，我们可以把它改写成等价的 Verilog HDL。下面的顶层模块展示了如何例化这个 256 点 FFT 内核，端口采用 Altera FFT MegaCore 的标准 Avalion-ST 流接口（原书 VHDL 改写为 Verilog）：

// 原书 VHDL 例化封装改写为 Verilog（可综合风格）
module fft256_wrap (
    input  wire        clk,           // 工作时钟
    input  wire        reset_n,       // 低有效复位
    // 输入数据流接口（sink 端）
    input  wire        sink_valid,    // 数据有效标志
    input  wire        sink_sop,      // 输入帧起始（start of packet）
    input  wire        sink_eop,      // 输入帧结束（end of packet）
    output wire        sink_ready,    // 内核准备好接收
    input  wire [15:0] sink_real,     // 输入实部
    input  wire [15:0] sink_imag,     // 输入虚部
    // 输出数据流接口（source 端）
    output wire        source_valid,
    output wire        source_sop,
    output wire        source_eop,
    input  wire        source_ready,
    output wire [15:0] source_real,   // 输出实部
    output wire [15:0] source_imag,   // 输出虚部
    output wire [5:0]  source_exp     // 块指数（如 -3 表示缩放 1/8）
);
    fft256 u_fft256 (
        .clk          (clk),
        .reset_n      (reset_n),
        .sink_valid   (sink_valid),
        .sink_sop     (sink_sop),
        .sink_eop     (sink_eop),
        .sink_ready   (sink_ready),
        .sink_real    (sink_real),
        .sink_imag    (sink_imag),
        .source_valid (source_valid),
        .source_sop   (source_sop),
        .source_eop   (source_eop),
        .source_ready (source_ready),
        .source_real  (source_real),
        .source_imag  (source_imag),
        .source_exp   (source_exp)
    );
endmodule

拿到输出后，还要根据块指数把数值”扶正”。用筒式移位器（barrel shifter，可在 1 个时钟内完成任意位数的移位）把实部和虚部同时按指数移位，即可去掉块缩放（原书 VHDL 缩放恢复逻辑改写为 Verilog）：

// 原书 VHDL 缩放恢复逻辑改写为 Verilog
// source_exp 为补码表示的块指数，左移即放大 2^exp 倍
wire signed [15:0] real_fix = source_real <<< source_exp;
wire signed [15:0] imag_fix = source_imag <<< source_exp;

图 6-21 和图 6-22 给出了 FFT 模块的仿真结果，可以一步步读懂握手过程。复位 reset 变低后，数据源置起数据可用信号 sink_valid；一个时钟周期后，帧起始标识 sink_sop 有效，同时第一个输入数据（测试中的值 20）随时钟送入内核。256 个时钟周期后，全部输入数据被接收完毕。由于选用 Streaming 模式没有额外等待，第一个输出数据出现在 $256\times 2\times 10\mathrm{ns}\approx 5\mu \mathrm{s}$ 处，对应 source_sop 脉冲，见图 6-22(a)。再经过 256 个周期，全部输出传完，source_eop 上出现脉冲标志帧结束（见图 6-22(b)），随后内核立刻开始处理下一个数据块。

图 6-21 IP FFT 模块初始化步骤的 Quartus II 仿真

(a) 输出帧的开始

(b) 输出帧的结束
图 6-22 FFT 内核仿真输出结果

还有一个细节值得注意：仿真得到的输出值看似几乎没有量化损失，但整体被缩小了 $2^{-3}$（即 1/8）。这不是误差，而是块浮点格式的结果——多级 FFT 在模块内部按块调整指数，把量化噪声降到最低。拿到结果后按指数值用筒式移位器把所有实部、虚部移回去，就恢复了无缩放的真实数值。

最后看实测资源：该设计运行速度 213.31 MHz（TimeQuest 缓慢 85C 模型），用了 4811 个 LE、18 个 9×9 位嵌入式乘法器（等价于 9 个 18×18 位乘法器）和 20 个 M9K 存储块。与工具测试台给出的估算（11 个 M9K、18 个乘法器、4581 个 LE）对照：乘法器估算完全准确，LE 偏差仅 5%，而 M9K 估算偏差高达 81%。这提醒我们：IP 工具的资源估算只能当参考，存储器类资源的估算尤其不可靠，最终一定要以实际编译报告为准。

6.4 与傅立叶相关的变换

学完了 DFT 和 FFT，为什么还要讲 DCT 和 DST？先说它们是什么：离散余弦变换（DCT）和离散正弦变换（DST）是与 DFT 平行的另外两族变换，区别只在于”核函数”用的是纯余弦或纯正弦，而不是 DFT 中的复指数 $e^{-j2\pi nk/N}$。

再说为什么值得关注。DCT 有一个非常宝贵的性质：对自然图像这类强相关信号，DCT 的能量集中能力几乎逼近理论上最优的 Kahunen-Loève 变换（KLT）——信号能量会高度集中到极少数低频系数上，丢弃高频系数造成的失真很小。这就是 JPEG、H.261、H.263、MPEG 等图像/视频压缩标准全都采用 $8\times 8$ 的二维 DCT-II 的原因。反过来讲，如果不用 DCT 而直接用 DFT 做压缩，会碰到边界不连续带来的高频泄漏问题，压缩效率明显变差。

不过 DCT 也有短板：卷积定理对它不成立——不能像 DFT 那样”乘频谱再逆变换”来实现快速卷积。所以 DCT 不像 FFT 那样无处不在，它的主战场是变换域压缩与分析。

所有 DCT 都遵循同一个变换模式（Wang 观察到的）：

$$X[k]=\sum _{n}x[n]C_N^{n,k}\leftrightarrow x[n]=\sum _{k}X[k]C_N^{n,k}\text{\hspace{1em}(6-55)}$$

四种 DCT 的核函数 $C_N^{n,k}$ 定义如下（其中除 $c[0]=1/\sqrt{2}$ 外，$c[m]=1$）：

DCT-I:

$$C_N^{n,k}=\sqrt{2/N}c[n]c[k]\cos (nk\frac{\pi }{N})n,k=0,1,\dots ,N$$

DCT-II:

$$C_N^{n,k}=\sqrt{2/N}c[k]\cos (k(n+\frac{1}{2})\frac{\pi }{N})n,k=0,1,\dots ,N-1$$

DCT-III: $C_N^{n,k}=\sqrt{2/N}c[n]\cos (n(k+\frac{1}{2})\frac{\pi }{N})$，$n,k=0,1,\dots ,N-1$

$$\text{DCT-IV:}{C}_N^{n,k}=\sqrt{2/N}\cos ((k+\frac{1}{2})(n+\frac{1}{2})\frac{\pi }{N})n,k=0,1,\dots ,N-1$$

DST 的结构与 DCT 完全相同，只是把余弦项全部换成正弦项。DCT 的重要性质归纳如下：

1. DCT 用余弦基实现函数展开。
2. 所有变换均是正交的，即 $C\times C^t=k[n]I$。
3. 与 DFT 不同，DCT 是实变换——实数输入得到实数输出，不需要复数运算，硬件实现代价天然减半。
4. DCT-I 是其本身的逆矩阵。
5. DCT-II 与 DCT-III 互为逆变换。
6. DCT-IV 是其本身的逆矩阵，且对称：$C=C^t$。
7. DCT 的卷积性质与 DFT 的卷积乘法关系不一样（这正是前面说的短板）。
8. DCT 是 KLT 的一种近似。

实际做图像压缩时用的是二维 $8\times 8$ 的 DCT-II。好消息是二维变换是可分离的：可以先把每一行做一维 DCT，再对每一列做一维 DCT（反过来也行）。这样二维问题就降解为 $2N$ 次一维变换，硬件上只需集中实现一个一维 DCT 电路，行列复用即可。

6.4.1 利用 DFT 计算 DCT

DCT 既然与 DFT 同源，能否”借 FFT 的东风”来算 DCT？答案是肯定的。Narasimha 和 Peterson 给出了 DCT 到 DFT 的映射模式。这个思路非常具有吸引力：一旦映射到 DFT，各种现成的 FFT 算法（基 2、分裂基、IP 核）都可以直接拿来用。

由于 DCT-II 最常用（图像压缩用的就是它），下面专门推导 DFT 与 DCT-II 的关系。为了记号简洁，先省略缩放因子——它可以在 DFT/FFT 计算的末尾一并补上。设变换长度 $N$ 为偶数，先做一个”奇偶重排”置换：

$$y[n]=x[2n]\text{和}y[N-n-1]=x[2n+1],\text{其中}n=0,1,\dots ,N/2-1$$

直观理解：把原序列的偶数下标样本依次放到 $y$ 的前半段，奇数下标样本倒序放到后半段。以 $N=8$ 为例，$x=\{x_0,x_1,\dots ,x_7\}$ 重排后得到 $y=\{x_0,x_2,x_4,x_6,x_7,x_5,x_3,x_1\}$。这个”先排偶、再倒排奇”的动作，让余弦相位里的 $(n+\frac{1}{2})$ 项正好凑成 $2n+\frac{1}{2}$ 与 $2n+\frac{3}{2}$ 两段，于是 DCT-II：

$$X[k]=\sum _{n=0}^{N-1}x[n]\cos \left(k(n+\frac{1}{2})\frac{\pi }{N}\right)\text{\hspace{1em}(6-56)}$$

被拆成两部分之和：

$$\begin{array}{l}X[k]={\sum }_{n=0}^{N/2-1}y[n]\cos \left(k(2n+\frac{1}{2})\frac{\pi }{N}\right)+{\sum }_{n=0}^{N/2-1}y[N-n-1]\cos \left(k(2n+\frac{3}{2})\frac{\pi }{N}\right)\\ X[k]={\sum }_{n}y[n]\cos \left(k(2n+\frac{1}{2})\frac{\pi }{N}\right)\end{array}\text{\hspace{1em}(6-57)}$$

现在计算 $y[n]$ 的 DFT，记为 $Y[k]$。把余弦写成复指数的实部，对比后发现一个干净利落的结果：

$$\begin{array}{rl}X[k]& =\Re (W_{4N}Y[k])\\ & =\cos \left(\frac{\pi k}{2N}\right)\Re (Y[k])-\sin \left(\frac{\pi k}{2N}\right)\Im (Y[k])\end{array}\text{\hspace{1em}(6-58)}$$

也就是说：DCT-II = 一个 $N$ 点 FFT 的输出，再乘一个复旋转因子 $W_{4N}$，取实部。计算量从 $O(N^2)$ 降到 $O(N\log N)$，只多出 $N$ 次额外乘法。式(6-58)很容易转写成程序（参阅练习 6.17），借助于 DFT 或 FFT 就可以计算 DCT。原书给出的 MATLAB 参考实现如下（先重排、再 FFT、最后逐点乘权重）：

function X = DCTII(x)
    N = length(x); % get length
    Y = [x(1:2:N), x(N:-2:2)]; % re-order elements
    Y = fft(y); % Compute the FFT
    w = 2 * exp(-i * (0:N-1)' * pi / (2*N)) / sqrt(2*N); % get weights
    w(1) = w(1) / sqrt(2); % make it unitary
    X = real(w.*Y); % compute pointwise product

逐行解读：第 2 行完成前面说的奇偶重排（偶数下标正序 + 奇数下标倒序）；第 3 行做 $N$ 点 FFT；第 4 行构造 $N$ 个复数权重 $w[k]$（对应 $W_{4N}$ 并带上归一化因子）；第 5 行单独处理 $k=0$ 的直流分量（补 $1/\sqrt{2}$，使变换成为幺正变换）；最后一行逐点相乘并取实部，即完成式(6-58)。

6.4.2 快速直接 DCT 实现

借 FFT 算 DCT 固然方便，但能否像 Cooley-Tukey 那样为 DCT 量身定做一个快速算法？Lee 给出了肯定的答案。由于该算法与基 2 Cooley-Tukey FFT 在结构上高度相似，被称为快速 DCT（Fast DCT, FCT），同样可以用矩阵结构开发。

推导从逆变换入手。因为 DCT 是正交变换，所以可以通过转置逆 DCT（IDCT）得到 DCT。式(6-55)引入的 IDCT-II 为：

$$x[n]=\sum _{k=0}^{N-1}\hat{X}[k]C_N^{n,k},n=0,1,\dots ,N-1\text{\hspace{1em}(6-59)}$$

注意其中 $\hat{X}[k]=c[k]X[k]$（直流分量要带上那个特殊的 $1/\sqrt{2}$ 权重）。把 $x[n]$ 分解成偶数部分和奇数部分后可以发现，$x[n]$ 可以由两个 $N/2$ 点 DCT 重构。先在频域定义两个半长序列：

$$G[k]=\hat{X}[2k]\text{\hspace{1em}(6-60)}$$ $$H[k]=\hat{X}[2k+1]+\hat{X}[2k-1],k=0,1,\dots ,N/2-1\text{\hspace{1em}(6-61)}$$

即 $G[k]$ 取偶数下标的频谱，$H[k]$ 取相邻奇数下标频谱之和。对应的时域半长变换为：

$$g[n]=\sum _{k=0}^{N/2-1}G[k]C_{N/2}^{n,k}\text{\hspace{1em}(6-62)}$$ $$h[n]=\sum _{k=0}^{N/2-1}H[k]C_{N/2}^{n,k},k=0,1,\dots ,N/2-1\text{\hspace{1em}(6-63)}$$

最后用”加减重构”把两个半长结果拼回全长序列：

$$x[n]=g[n]+1/(2C_N^{n,k})h[n]\text{\hspace{1em}(6-64)}$$ $$x[N-1-n]=g[n]-1/(2C_N^{n,k})h[n],n=0,1,\dots ,N/2-1\text{\hspace{1em}(6-65)}$$

一眼就能认出这副”骨架”：一个 $N$ 点问题拆成两个 $N/2$ 点子问题，再加减重组——与 DIT 基 2 FFT 的分解思想如出一辙。区别在于 FFT 里蝶形输出只需乘/加减旋转因子 $W_N^{nk}$，而这里出现的是除法 $1/(2C_N^{n,k})$。把图 6-13 的基 2 FFT 旋转因子与式(6-62)对比即可看出这一点。除法在硬件里代价高昂，但好在这些系数是固定常数：应该预先计算并存储在查找表中。这种制表方法对 Cooley-Tukey FFT 同样适用，因为在线计算三角函数一般非常耗时。

接下来用 8 点的例子把整个算法走一遍。

例 6.21：8 点 FCT。 对于 8 点 FCT，式(6-60)~(6-65)具体化为：

$$G[k]=\hat{X}[2k]\text{\hspace{1em}(6-66)}$$ $$H[k]=\hat{X}[2k+1]+\hat{X}[2k-1],k=0,1,2,3\text{\hspace{1em}(6-67)}$$

在时域中就有两个 4 点 DCT：

$$g[n]=\sum _{k=0}^{3}G[k]C_4^{n,k}\text{\hspace{1em}(6-68)}$$ $$h[n]=\sum _{k=0}^{3}H[k]C_4^{n,k},n=0,1,2,3\text{\hspace{1em}(6-69)}$$

这样，重构就变成：

$$x[n]=g[n]+1/(2C_8^{n,k})h[n]\text{\hspace{1em}(6-70)}$$ $$x[N-1-n]=g[n]-1/(2C_8^{n,k})h[n],n=0,1,2,3\text{\hspace{1em}(6-71)}$$

式(6-66)和(6-67)构成了流程图的第一级（频域的偶/奇拆分与相邻相加），而式(6-70)和(6-71)构成了流程图的最后一级（乘系数 $c[p]$ 后的加减重组）。由于 DCT 是正交变换，把上述流程图的所有箭头反向、信号反向流动，就得到正向 DCT 的流程。完整结构如图 6-23 所示，图中使用速记符号 $c[p]=1/(2\cos (p\pi /16))$。

图 6-23 采用速记符号 $c[p]=1/(2\cos (p\pi /16))$ 的 8 点快速 DCT 流程图

流程图里还藏着两个”排序机关”，初学者容易在这里栽跟头：

输入侧：进入流程图的序列 $\hat{X}[k]$ 是位逆序排列的——这与基 2 FFT 的输出排序规律一致，即 8 个下标按二进制位颠倒的顺序（0,4,2,6,1,5,3,7）排列。

输出侧：$x[n]$ 的顺序按”前缀生成法”排列：从集合 $(0,1)$ 出发，每轮通过给现有元素增加前缀 0 和 1 来生成新集合；前缀为 1 时，还要把前面模式中所有的位都颠倒。例如从序列 10 出发，得到两个子序列 010 和 $1\overline{1}\overline{0}=101$（前缀 1，后两位 10 取反成 01，合并即 101）。图 6-24 给出了这种生成模式的图解。写成硬件时，这些置换都是固定的连线关系，不消耗任何逻辑资源，只是布线时要照表接线。

图 6-24 8 点快速 DCT 的输入和输出置换

与 8 点直接 DCT 需要 $8\times 8=64$ 次乘法相比，图 6-23 的 FCT 流程只需约 13 次乘法——这正是”快速”二字的分量。练习 6.31 和 6.32 将引导你把这个流程写成软件和 HDL 代码，并在 FPGA 上实现一个带串行 I/O 的 8 点 IDCT。

6.5 练习

本节练习覆盖本章全部内容：从基础 DFT 计算、各种快速算法（Goertzel、CZT、Rader、Winograd、CFA、PFA、Good-Thomas），到 IP 核应用、DCT/DHT 变换以及 FPGA 设计实践。做 HDL 类练习时请使用 Cyclone IVE 系列 EP4CE115F29C7 器件，并用 TimeQuest 缓慢 85C 模型报告 $F_{\max }$ 和资源占用（LE、嵌入式乘法器、M9K）。没有 Quartus II 经验的读者可先复习 1.4.3 节的案例研究。

基础与卷积类

• 6.1 用傅立叶变换计算矩形窗和三角窗的 3dB 带宽、第一个零点、最大旁瓣和每倍频程衰减。
• 6.2 与 6.24：计算 $x[n]=\{3,1,-1\}$、$f[n]=\{2,1,5\}$（或 5 点序列 $\{1,2,3,4,5\}$ 与 $\{-1,0,-2,0,4\}$）的循环卷积；写出 $W_3$（或 $W_5$）DFT 矩阵；分别求 DFT，再按 $y=W^{-1}Y$ 验证 $Y[k]=X[k]F[k]$。建议用 C 或 MATLAB。

快速算法类

• 6.3 与 6.19：Goertzel 算法把单边频谱分量改写为嵌套形式 $X[k]=x[0]+W_N^k(x[1]+W_N^k(\cdots ))$，构成一阶递归结构，只需少数频点时非常划算。要求画 $N=5$ 的递归信号流程图（含输入/输出寄存器），并用三组输入 $\{20,40,60,80,100\}$、$\{j20,\dots ,j100\}$、$\{20+j20,\dots \}$ 在 $k=1$ 时计算各寄存器内容；6.19 要求用 8 位系数量化在 Quartus II 上实现并测资源与时序。
• 6.4：确定 $N=4$ 的 Bluestein Chirp-z（CZT）算法，用 MATLAB 计算三角序列 $\{0,1,2,3\}$ 的 CZT，再扩展到 $N=256$ 并与同长 FFT 交叉检验。
• 6.5：对 $N=7$ 的 Rader 算法设计非递归滤波器的直接实现，找出可组合的系数并比较工作量。
• 6.6：设计长度 $N=3$ 的 Winograd DFT 算法并画出信号流程图。
• 6.7 与 6.8：练习索引映射双射性的判定。例如 $n=3n_1+2n_2\mathrm{mod}6$（$N_1=2,N_2=3$）是否双射？Burrus 条件 $\gcd (A,N_1)=\gcd (B,N_2)=1$ 如何用于 $N_1=5,N_2=8$ 时 $A=8,B=5$、$A=8,B=10$、$A=24,B=15$ 等情形？并求 $A=7$ 或 $A=8$ 时所有有效的 $B\in Z_{20}$。
• 6.9/6.10：绘制 $N=16$ 的 DIF、$N=8$ 的 DIT 基 2 FFT 信号流程图，编写 C/MATLAB 程序，用复三角输入 $x[n]=n+jn$ 验证。
• 6.11 与 6.25：公共因子算法（CFA）练习。按索引映射 $n=4n_1+n_2$、$k=k_1+4k_2$（$N=16$）或 $n=3n_1+n_2$、$k=k_1+5k_2$（$N=15$，$N_1=5,N_2=3$）编制映射表，补全流程图，并分四步（映射输入并算第一级 DFT → 乘旋转因子 → 算第二级 DFT → 排序输出）为给定稀疏序列手算 16 点/15 点 FFT。

图 6-25 未完成的 16 点基 4 FFT 的信号流程图

图 6-26 未完成的 15 点 CFA FFT 信号流程图

• 6.12：绘制 $N=12$ 的 Good-Thomas FFT 信号流程图，要求图中不出现交叉（提示：用行/列 DFT 的三维表示）。
• 6.13：对 Burrus-Eschenbacher 索引变换 $n=N_2{n}_1+N_1{n}_2\mathrm{mod}N$、$k=N_2{k}_1+N_1{k}_2\mathrm{mod}N$（$N_1=3,N_2=4$），算出映射表与 $W^{nk}$，代入 DFT 矩阵，判断算法类型，并讨论 Rader 算法的适用性。
• 6.14：计算 $W_2$、$W_3$ 矩阵及 Kronecker 积 $W_6^{\prime }=W_2\otimes W_3$，确定使 $X=W_6^{\prime }x$ 成为 6 点 DFT 的输入/输出索引，再对 IDFT 情形重复。

DCT/DHT 类

• 6.15：离散 Hartley 变换 $H[n]={\sum }_k\mathrm{cas}(2\pi nk/N)h[k]$（$\mathrm{cas}(x)=\sin x+\cos x$）与 DFT 通过 $H[n]=\Re \{F[n]\}-\Im \{F[n]\}$、$F[n]=E[n]-jO[n]$、$E/O$ 取 $H$ 的偶/奇部分相联系。要求给出逆 DHT、用 DHT 求卷积的步骤，以及输入为偶序列时的简化形式。
• 6.16：对 DCT-II 标准形式（式(6-81)(6-82)）给出逆变换方程、$N=4$ 的 DCT 矩阵，手算 $\{1,2,2,1\}$ 与 $\{1,1,-1,-1\}$ 的变换，并总结奇/偶对称序列的 DCT 特点。
• 6.17 与 6.18：用 6.4.1 节的 DCTII 程序计算 $\{1,2,2,1\}$ 与 $\{1,1,-1,-1\}$；再利用 DCT 的可分离性，对 $x[k,l]=k+l$ 的 $8\times 8$ 数据分别按”先行后列""先列后行”及式(6-83)直接实现三种方式计算二维 DCT 并比较结果。
• 6.31：快速 IDCT 软件设计。按图 6-23 编写 8 点 IDCT 代码（注意 $X[0]$ 的 $\sqrt{1/2}$ 缩放与式(6-55)的 $\sqrt{2/N}$ 缩放未在图中画出），用 MATLAB 的 idct 函数以 $X=10,20,\dots ,80$ 验证；用 abs、max、dctmtx、sum 估计每个频谱分量的最大位增长率；用 csd.exe 求各系数 $c[p]=0.5/\cos (p\pi /16)$ 的最少 8 位 CSD 表示；分别以浮点与整型（含 4 个保护位，缩放 $2^4=16$）计算输出并列出三级乘法后的中间值。
• 6.32：8 点 IDCT 的 HDL 设计。按图 6-23 编写 HDL 代码，带异步复位、ena 使能与串行 I/O；输入 8 位，输出与内部数据用 14 位整型格式（4 个小数位，即输入乘 16、输出除 16）。用 6.31(f) 的数据调试并与图 6-30 的仿真结果比对，报告 $F_{\max }$ 与资源，确定 HDL 与软件结果的最大相对输出误差百分比。

图 6-30 8 点 IDCT 的 VHDL 仿真结果

HDL/FPGA 设计类

• 6.20/6.21：用 Quartus II 分别设计实数输入 4 点 Winograd DFT（输入/输出精度 8/10 位）和复数输入 3 点 Winograd DFT（10/12 位）的 Component，各带输入/输出寄存器与同步使能信号；分别用 $\{40,70,100,10\}$ 等三组、$\{180,220,260\}$ 输入仿真。
• 6.22：把 6.20 与 6.21 的 3 点、4 点 Component 用元件例化组成类似图 6-18 的完全并行 12 点 Good-Thomas FFT（8/12 位精度），用 $x[n]=10n$（$0\le n\le 11$）仿真。
• 6.23：设计与算法 6.10 类似的旋转因子复数乘法器 ccmulp：乘法器 3 级流水、输入减法 $X-Y$ 1 级流水、8 位精度；验证 $(70+j50)(121+j39)$ 计算正确；再实现完整流水线蝶形处理器，用例 6.11 数据仿真并报告时序与资源。
• 6.28/6.29：用 HDL 分别设计 5 点实数输入 Winograd DFT（8/11 位，系数用 csd.exe 量化为最少 8 位 CSD 编码）和 2 点复数输入 Winograd DFT（11/12 位），均带寄存器与同步使能；用 $\{10,30,50,70,90\}$、$\{60,80,100,20,40\}$ 和 $\{250,300\}$、$\{-50+j67,-j85\}$ 仿真，并与图 6-27、图 6-28 的结果匹配。

图 6-27 5 点实数输入的 Winograd DFT 的仿真结果

图 6-28 两点复数输入的 Winograd DFT 的仿真结果

• 6.30：综合练习。将 6.28、6.29 的 5 点与 2 点元件按元件例化组成完全并行的 10 点 Good-Thomas FFT（8/12 位），I/O 状态机和寄存器用异步复位，一组 I/O 值传完时由 ENA 信号指示；用 $x[n]=10n$（$0\le n\le 10$）仿真并与图 6-29 的开始帧/结束帧比对，报告时序与资源。

(a) 开始帧

(b) 结束帧
图 6-29 10 点 Good-Thomas FFT 的仿真结果