📘 第 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 个。