注：本文为 “快速傅里叶变换（FFT）” 相关合辑。
图片清晰度受引文原图所限。
略作重排，未整理去重。
如有内容异常，请看原文。

---

快速理解 FFT 算法

cronus-裁 发布于 2021-10-23 21:19

0. 引言

学习 FFT 算法应优先聚焦其傅里叶变换本质，而非高关注度资料中提及的多项式系数表示法与点值表示法。这类表示法仅为 FFT 算法的应用场景之一，若未先明晰傅里叶变换的物理意义便学习 FFT，易导致认知混淆，无法达成对算法的完整理解。

逻辑梳理：FFT 算法是 DFT 算法的优化实现形式，DFT 算法是傅里叶变换的离散化表达形式。遵循"傅里叶变换 $\to$ DFT $\to$ FFT"的推导路径，可系统理解 FFT 算法的设计原理。

1. 傅里叶变换的物理意义

为简化表述，本节省略傅里叶变换的完整数学推导过程，仅阐释其物理意义。

若已掌握傅里叶变换的物理意义，可跳过本节直接阅读 2. DFT 部分；若尚未理解傅里叶变换的基本概念，建议参考傅里叶变换专项讲解资料或高等数学教材。

1.1 周期函数与非周期函数的傅里叶分析

周期函数的傅里叶级数本质是将函数 $f(t)$ 分解为无穷多个不同频率、不同幅值的正弦与余弦信号的叠加，且这些信号的频率均为基频 ${\omega }_0$ 的整数倍（即 $n{\omega }_0$）。通过无穷多个不同频率、不同幅值的正余弦信号线性组合，可逼近周期函数 $f(t)$；分解过程中，每个频率分量 $n{\omega }_0$ 对应一个幅值，该幅值与频率 $n{\omega }_0$ 构成的函数关系（自变量为 $n{\omega }_0$，因变量为对应频率信号的幅值）被定义为频谱函数。

非周期函数的傅里叶变换以频谱密度函数为求解目标，该函数以角频率 $\omega$ 为自变量，因变量为信号幅值在频域中的分布密度（即单位频率区间内的信号强度）。（若具备概率论基础，可将频谱函数与频谱密度函数分别类比为离散型随机变量的概率质量函数与连续型随机变量的概率密度函数。）

若上述物理意义未完全理解，可继续阅读下文；若已理解，可跳转至 2. DFT 部分。

1.2 傅里叶变换的数学推导

将周期函数的傅里叶级数转化为指数形式，表达式如下：

$$\begin{array}{rl}{c}_n& =\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t){\mathrm{e}}^{-\mathrm{j}n{\omega }_{0}t}\mathrm{d}t\\ f(t)& =\sum _{n=-\infty }^{+\infty }{c}_n{\mathrm{e}}^{\mathrm{j}n{\omega }_{0}t}\end{array}$$

傅里叶级数的周期为 $T$，无法直接用于表示非周期函数。针对非周期函数，将其视为周期 $T\to \infty$ 的极限情况，此时基频 ${\omega }_0=\frac{2\pi }{T}\to 0$；相邻频率分量的差值 $(n+1){\omega }_0-n{\omega }_0={\omega }_0$ 记为 $\Delta \omega$，在极限情况下成为微分元 $\mathrm{d}\omega$，离散频率 $n{\omega }_0$ 对应连续变量 $\omega$。由此推导非周期函数的频谱密度函数：

令 $c_n=\frac{\Delta \omega }{2\pi }F(\omega )$，其中 $F(\omega )={\int }_{-\infty }^{+\infty }f(\tau ){\mathrm{e}}^{-\mathrm{j}\omega \tau }\mathrm{d}\tau$（为避免混淆，积分变量改用 $\tau$）。

$\tau$ 是积分哑变量（dummy variable / bound variable），可任意替换，仅存在于积分号内，积分完成后"消失"；$f(\tau )$ 中的 $\tau$ 为被积函数的求值点，由积分遍历积分区间所有值。

当 $T\to \infty$ 时，$c_n\to 0$，将其代入 $f(t)={\sum }_{n=-\infty }^{+\infty }{c}_n{\mathrm{e}}^{\mathrm{j}n{\omega }_{0}t}$，可得：

$$\begin{array}{rl}f(t)& =\sum _{n=-\infty }^{+\infty }\left(\frac{\Delta \omega }{2\pi }{\int }_{-\infty }^{+\infty }f(\tau ){\mathrm{e}}^{-\mathrm{j}n{\omega }_0\tau }\mathrm{d}\tau \right){\mathrm{e}}^{\mathrm{j}n{\omega }_{0}t}\\ & \to \frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\left({\int }_{-\infty }^{+\infty }f(\tau ){\mathrm{e}}^{-\mathrm{j}\omega \tau }\mathrm{d}\tau \right){\mathrm{e}}^{\mathrm{j}\omega t}\mathrm{d}\omega \end{array}$$

据此定义非周期函数的傅里叶变换：

$$F(\omega )={\int }_{-\infty }^{+\infty }f(t){\mathrm{e}}^{-\mathrm{j}\omega t}\mathrm{d}t$$

由推导过程可知 $c_n=\frac{\mathrm{d}\omega }{2\pi }F(\omega )$，即傅里叶变换选取信号幅值在频域中的分布密度 $F(\omega )$ 作为表征量（而非对应频率的信号幅值 $c_n$）。若直接选取 $c_n$，在 $T\to \infty$ 时函数值均趋于 0，无实际应用价值。

傅里叶变换也可采用频率 $f$（单位：Hz）作为自变量，表达式为：

$$F(f)={\int }_{-\infty }^{+\infty }f(t){\mathrm{e}}^{-\mathrm{j}2\pi ft}\mathrm{d}t$$

1.3 结论

傅里叶变换的目标是将时域信号分解为无穷多个不同频率正弦信号的叠加，并揭示各正弦信号的强度与频率之间的对应关系。

2. DFT（离散傅里叶变换）

若已掌握 DFT 的基本原理，可跳过本节直接阅读 3. FFT 部分。

傅里叶变换的目标是求解信号的频谱密度函数（自变量为 $\omega$，因变量为信号幅值在频域中的分布密度），该函数反映信号强度与频率的对应关系。其基本表达式为：

$$F(\omega )={\int }_{-\infty }^{+\infty }f(t){\mathrm{e}}^{-\mathrm{j}\omega t}\mathrm{d}t$$

理论分析中 $f(t)$ 多为时域连续函数，但计算机采集的信号经采样处理后为时域离散信号，且计算机仅能计算有限个频率点的频谱值（频率需离散化）。DFT 是时域与频域均离散化的傅里叶变换形式。

2.1 采样的数学原理

若已掌握采样的基本原理，可跳过本节直接阅读 2.2 部分。

采样的数学实现依赖冲激函数（Dirac 函数），其定义为：

• 当 $t\ne 0$ 时，$\delta (t)=0$；
• 积分满足 ${\int }_{-\infty }^{+\infty }\delta (t)\mathrm{d}t=1$。

冲激函数的筛选性质：${\int }_{-\infty }^{+\infty }\delta (t-t_0)f(t)\mathrm{d}t=f(t_0)$，可提取函数 $f(t)$ 在 $t_0$ 处的函数值，实现单点采样。

单点采样无法满足时域离散化需求，需将多个不同时刻的冲激函数叠加构成梳状函数（理想采样函数）：

$${\delta }_T(t)=\sum _{n=-\infty }^{+\infty }\delta (t-n{T}_s)$$

其中 $T_s$ 为采样周期，。

将时域连续信号与梳状函数相乘，得到时域离散的采样信号：

$$f_s(t)=f(t)\cdot {\delta }_T(t)=\sum _{n=-\infty }^{+\infty }f(n{T}_s)\delta (t-n{T}_s)$$

2.2 时域离散化计算

对采样得到的离散信号 $x_s(t)={\sum }_{n=-\infty }^{+\infty }x(n{T}_s)\delta (t-n{T}_s)$ 进行傅里叶变换。

根据傅里叶变换定义：

$$X(\omega )={\int }_{-\infty }^{+\infty }x(t){\mathrm{e}}^{-\mathrm{j}\omega t}\mathrm{d}t$$

采样信号的傅里叶变换为：

$$X_s(\omega )={\int }_{-\infty }^{+\infty }\left(\sum _{n=-\infty }^{+\infty }x(n{T}_s)\delta (t-n{T}_s)\right){\mathrm{e}}^{-\mathrm{j}\omega t}\mathrm{d}t$$

交换积分与求和顺序：

$$X_s(\omega )=\sum _{n=-\infty }^{+\infty }x(n{T}_s){\int }_{-\infty }^{+\infty }\delta (t-n{T}_s){\mathrm{e}}^{-\mathrm{j}\omega t}\mathrm{d}t$$

根据冲激函数筛选性质，最终得到时域离散化结果（即离散时间傅里叶变换，DTFT）：

$$X_s(\omega )=\sum _{n=-\infty }^{+\infty }x(n{T}_s){\mathrm{e}}^{-\mathrm{j}\omega n{T}_s}$$

或简记为（令 $x[n]=x(n{T}_s)$，$\Omega =\omega T_s$）：

$$X({\mathrm{e}}^{\mathrm{j}\Omega })=\sum _{n=-\infty }^{+\infty }x[n]{\mathrm{e}}^{-\mathrm{j}\Omega n}$$

2.3 频域离散化计算

时域离散化得到的 $X_s(\omega )$ 仍为频域连续函数（周期为 ${\omega }_s=2\pi /T_s$），且对应无限长序列 $x[n]$，而计算机只能存储有限长序列并计算有限个频率点，不符合工程需求。

2.3.1 处理思路

对信号 $x(t)$ 进行 $N$ 次（有限值）采样，得到有限长序列 $x[n]$，$n=0,1,\dots ,N-1$。DFT 隐含假设该有限长序列是某个周期序列的一个周期（时域周期延拓），周期为 $T_0=N{T}_s$。周期序列的频谱为离散谱（傅里叶级数），由此实现频域离散化。

注：时域周期延拓是 DFT 的数学构造，并非实际物理操作。

2.3.2 数学推导

一个周期内离散信号的表达式：

$$\tilde{x}(t)=\sum _{n=0}^{N-1}x(n{T}_s)\delta (t-n{T}_s)$$

周期信号 $\tilde{x}(t)$ 的傅里叶级数系数为：

$$\tilde{X}(k{\omega }_0)=\frac{1}{T_0}{\int }_0^{T_0}\tilde{x}(t){\mathrm{e}}^{-\mathrm{j}k{\omega }_{0}t}\mathrm{d}t$$

其中基频 ${\omega }_0=\frac{2\pi }{T_0}$。

代入 $\tilde{x}(t)$ 并交换积分与求和顺序：

$$\tilde{X}(k{\omega }_0)=\frac{1}{T_0}\sum _{n=0}^{N-1}x(n{T}_s){\int }_0^{T_0}\delta (t-n{T}_s){\mathrm{e}}^{-\mathrm{j}k{\omega }_{0}t}\mathrm{d}t$$

根据冲激函数筛选性质：

$$\tilde{X}(k{\omega }_0)=\frac{1}{T_0}\sum _{n=0}^{N-1}x(n{T}_s){\mathrm{e}}^{-\mathrm{j}k{\omega }_{0}n{T}_s}$$

代入 $T_0=N{T}_s$ 与 ${\omega }_0=\frac{2\pi }{T_0}$ 化简：

$$\begin{array}{rl}\tilde{X}(k{\omega }_0)& =\frac{1}{N{T}_s}\sum _{n=0}^{N-1}x(n{T}_s){\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N{T}_s}kn{T}_s}\\ & =\frac{1}{N{T}_s}\sum _{n=0}^{N-1}x(n{T}_s){\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}kn}\end{array}$$

定义简化符号：

• $X[k]=N{T}_s\cdot \tilde{X}(k{\omega }_0)$（或根据归一化方式不同，定义为 $X[k]=\tilde{X}(k{\omega }_0)\cdot T_s$ 等，需保持前后一致）
• $x[n]=x(n{T}_s)$

最终得到 DFT 表达式（采用无系数形式，与 FFT 文献常用形式一致）：

$$X[k]=\sum _{n=0}^{N-1}x[n]{\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}kn},k=0,1,\dots ,N-1$$

注：DFT 定义中系数 $1/N$ 可置于正变换或逆变换，或两边均分 $1/\sqrt{N}$，不同领域惯例不同。本文采用信号处理领域常见形式（正变换无系数，逆变换含 $1/N$）。

2.4 DFT 的复杂度

对于给定的 $k$ 值，求解 $X[k]$ 需完成 $N$ 次复数乘法（$x[n]$ 与 ${\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}kn}$）和 $N-1$ 次复数加法。

由于 ${\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}kn}$ 对 $k$ 具有周期性（周期为 $N$），实际计算仅需 $k=0,1,\dots ,N-1$。因此 DFT 算法的时间复杂度为 $\mathcal{O}(N^2)$。

3. FFT（快速傅里叶变换）

3.1 数学原理

DFT 表达式为 $X[k]={\sum }_{n=0}^{N-1}x[n]{\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}kn}$，时间复杂度为 $\mathcal{O}(N^2)$。其中：

• $n$ 对应时域采样序号；
• $k$ 对应频率序号，实际频率为 $f_k=\frac{k}{N{T}_s}$ 或角频率 ${\omega }_k=\frac{2\pi k}{N{T}_s}$；
• $X[k]$ 对应频率 ${\omega }_k$ 处的频谱值。

3.1.1 旋转因子的周期性

降低 DFT 复杂度的思路是利用 ${\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}kn}$ 的周期性。定义旋转因子 $W_N^{nk}={\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}nk}$。

当 $n$ 为偶数（令 $n=2m$）时：

$$W_N^{2mk}={\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}2mk}={\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N/2}mk}=W_{N/2}^{mk}$$

且满足 $W_{N/2}^{mk}=W_{N/2}^{m(k+N/2)}$，即 $W_{N/2}^{mk}$ 对 $k$ 的周期为 $N/2$。

3.1.2 奇偶项拆分

将 $X[k]$ 按 $n$ 的奇偶性拆分为两部分：

$$X[k]=\sum _{m=0}^{N/2-1}x[2m]W_N^{2mk}+\sum _{m=0}^{N/2-1}x[2m+1]W_N^{(2m+1)k}$$

化简过程：

$$\begin{array}{rl}X[k]& =\sum _{m=0}^{N/2-1}x[2m]W_{N/2}^{mk}+W_N^k\sum _{m=0}^{N/2-1}x[2m+1]W_{N/2}^{mk}\end{array}$$

定义：

• 偶数项部分：$E[k]={\sum }_{m=0}^{N/2-1}x[2m]W_{N/2}^{mk}$；
• 奇数项部分：$O[k]={\sum }_{m=0}^{N/2-1}x[2m+1]W_{N/2}^{mk}$。

则 $X[k]$ 简化为：

$$X[k]=E[k]+W_N^{k}O[k]$$

3.1.3 周期性推导

由于 $W_{N/2}^{mk}$ 对 $k$ 的周期为 $N/2$，因此 $E[k]=E[k+N/2]$、$O[k]=O[k+N/2]$。

且旋转因子满足 $W_N^{k+N/2}=-W_N^k$，推导如下：

$$\begin{array}{rl}{W}_N^{k+N/2}& ={\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}(k+N/2)}={\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}k}{\mathrm{e}}^{-\mathrm{j}\pi }=-{\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}k}=-W_N^k\end{array}$$

由此可得：

$$X[k+N/2]=E[k]-W_N^{k}O[k]$$

以 $N=8$ 为例：

$$X[0]=E[0]+W_8^{0}O[0],X[4]=E[0]-W_8^{0}O[0]$$

即求解得到 $E[0]$ 与 $O[0]$ 后，可同时得到 $X[0]$ 与 $X[4]$，无需重复计算。

3.2 算法流程

以 $N=8$ 为例，求解 $X[0]$ 与 $X[4]$ 仅需计算 $E[0]$ 与 $O[0]$。
对任意 $k=0,1,\dots ,N/2-1$，已知频域奇偶分解结果 $E[k]$ 与 $O[k]$，可通过一次复数乘法与两次复数加减运算同时得到 $X[k]$ 与 $X[k+N/2]$。
由于 $E[k]$ 与 $O[k]$ 均以 $N/2$ 为周期，满足
[
E[k+N/2] = E[k],\quad O[k+N/2] = O[k]
]
因此 $X[k]$ 与 $X[k+N/2]$ 可共享同一组奇偶分解结果，只需计算一次 $E[k]$、$O[k]$ 即可完成一对频点的求解。

基于递归思想，对 $E[k](k=0,1,2,3)$ 继续拆分：

• $E[k]$ 对应偶数序号样本 $x[0],x[2],x[4],x[6]$，可视为长度 $N/2=4$ 的 DFT；
• 利用旋转因子关系
  $${\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}\cdot 2nk}={\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N/2}\cdot nk}$$
  将 $E[k]$ 进一步拆分为偶数项与奇数项，重复拆分操作直至每一部分仅包含单个样本。

关于蝶形图的详细解析可参考相关专项资料。

3.3 FFT 的复杂度

若对数据结构与算法、FFT 复杂度推导过程无学习需求，可跳过本节。

时间复杂度表征算法处理规模为 $N$ 的数据集时所需的运算次数，记为 $T(N)$。

3.3.1 递归关系推导

基于3.2节的拆分逻辑，FFT算法的时间复杂度满足递归关系：
$$T(N)=2T(N/2)+\mathcal{O}(N)$$
推导依据：将 $N$ 点 DFT 拆分为两个 $N/2$ 点 DFT（对应 $2T(N/2)$）；合并两个子问题结果时，需执行 $N$ 次复数加法和 $N/2$ 次复数乘法（$W_N^{k}O[k]$），这部分运算量为 $\mathcal{O}(N)$。

3.3.2 复杂度求解

复杂度分析的本质是“阶数”而非“精确数值”。

对递归公式展开（假设 $N=2^L$，$L={\log }_{2}N$）：
$$\begin{array}{rl}T(N)& =2T(N/2)+N\\ & =2\left(2T(N/4)+\frac{N}{2}\right)+N=4T(N/4)+2N\\ & =4\left(2T(N/8)+\frac{N}{4}\right)+2N=8T(N/8)+3N\\ & =\cdots \\ & =2^{L}T(1)+L\cdot N\end{array}$$

前提：单点 DFT（$N=1$ 时，仅需直接取样本值 $X[0]=x[0]$）无需任何复数乘法/加法运算，因此 $T(1)=0$。代入后可得：
$$T(N)=N{\log }_{2}N$$

由此可知，FFT 算法的时间复杂度为 $\mathcal{O}(N\log N)$。结合 FFT 蝶形运算的实际执行过程，可进一步明确具体运算次数：

• 复数乘法次数：$\frac{N}{2}{\log }_{2}N$（仅在合并奇偶子问题时需计算旋转因子与 $O[k]$ 的乘积，共 $N/2$ 次/层，总计 ${\log }_{2}N$ 层）；
• 复数加法次数：$N{\log }_{2}N$（每对 $X[k]$ 与 $X[k+N/2]$ 需 2 次加减运算，共 $N/2$ 对/层，总计 ${\log }_{2}N$ 层）。

4. 总结

1. 傅里叶变换的本质是将时域信号分解为不同频率复指数信号的叠加，DFT 是其时域与频域均离散化的工程实现，FFT 则是 DFT 的高效算法；
2. FFT 算法利用旋转因子的周期性与对称性，将 DFT 的 $\mathcal{O}(N^2)$ 复杂度降至 $\mathcal{O}(N\log N)$，手段为递归奇偶分解与蝶形运算复用；
3. FFT 算法通过分治与复用策略，大幅降低频域分析运算量，是数字信号处理中高效的频谱分析方法。

---

如何理解和掌握快速傅里叶变换的计算和概念？

杨宇翔 编辑于2014-05-10 18:29

离散傅里叶变换（DFT）的定义

离散傅里叶变换（Discrete Fourier Transform, DFT）的定义为：

$$X[k]=\sum _{n=0}^{N-1}x[n]{\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}kn},0\le k\le N-1$$

为表达简洁，令 $W_N={\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}}$，则上式可写为：

$$X[k]=\sum _{n=0}^{N-1}x[n]W_N^{kn}$$

FFT 的种类很多，本文以最简单的基 2 FFT（信号长度 $N$ 为 2 的整数幂）为例进行说明。

FFT 本质上是一种分治算法。该算法将长度为 $N$ 的信号分解为两个长度为 $\frac{N}{2}$ 的子信号进行处理，如此递归分解直至长度为 1。每一次分解均可减少计算量。理解 FFT 可分为以下三个步骤：

步骤 1：将信号 $x[n]$ 分解为两个子信号

• 偶数样本点信号：$x[2n]$
• 奇数样本点信号：$x[2n+1]$，其中 $n=0,1,\dots ,\frac{N}{2}-1$

则：

$$X[k]=\sum _{n=0}^{N-1}x[n]W_N^{kn}=\sum _{n=0}^{\frac{N}{2}-1}x[2n]W_N^{k(2n)}+\sum _{n=0}^{\frac{N}{2}-1}x[2n+1]W_N^{k(2n+1)}$$

步骤 2：将两个求和项理解为两个长度为 $\frac{N}{2}$ 的 DFT

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

其中 $W_N^{k(2n)}={\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N}\cdot 2n}={\mathrm{e}}^{-\mathrm{j}\frac{2\pi }{N/2}\cdot n}=W_{N/2}^{kn}$，因此：

$$X[k]=E[k]+W_N^{k}O[k],0\le k\le N-1$$

式中，
偶数样本点信号的 DFT 为：

$$E[k]=\sum _{n=0}^{\frac{N}{2}-1}x[2n]W_{N/2}^{kn}$$

奇数样本点信号的 DFT 为：

$$O[k]=\sum _{n=0}^{\frac{N}{2}-1}x[2n+1]W_{N/2}^{kn}$$

步骤 3：FFT 的具体计算过程（通过蝶形图可视化）

由于 $W_N^{kn}$ 对 $k$ 和 $n$ 均具有周期性：$W_N^{kn}=W_N^{k(n+N)}=W_N^{(k+N)n}$，故：

$$E[k]=E\left[k+\frac{N}{2}\right],O[k]=O\left[k+\frac{N}{2}\right]\text{（周期性）}$$

以 $N=8$ 为例，信号的 FFT 计算过程如下：

第 1 次分解：

第 2 次分解：

第 3 次分解：

FFT 与 DFT 的性能比较

根据 DFT 的定义：

• 对于任意 $k$，需进行 $N$ 次乘法运算，故 DFT 共需 $N^2$ 次乘法运算；
• 对于任意 $k$，需进行 $N-1$ 次加法运算，故 DFT 共需 $N(N-1)$ 次加法运算。

从 FFT 的蝶形图可知：

共需 $N({\log }_{2}N-1)$ 次乘法运算和 $N{\log }_{2}N$ 次加法运算。

---

如何理解和掌握快速傅里叶变换的计算和概念？

游凯超 编辑于 2023-06-18 09:17 北京

离散傅里叶变换（DFT）是信号处理中的重要方法，快速傅里叶变换（FFT）是其一种快速实现方式。本文从算法复杂度角度，阐述快速傅里叶变换"快"在何处。

从矩阵-向量乘法说起

对于大小为 $N$ 的矩阵 ${\mathbf{A}}_{N\times N}$，该矩阵作用于 $N$ 维向量 $\mathbf{x}$ 的矩阵-向量乘法 $\mathbf{A}\mathbf{x}$ 的复杂度显然为 $\mathcal{O}(N^2)$。

离散傅里叶变换 DFT

对长度为 $N$ 的向量 $\mathbf{x}$ 作离散傅里叶变换，本质上是用一个特殊的矩阵——离散傅里叶变换矩阵 ${\mathbf{D}}_{N\times N}$ 对 $\mathbf{x}$ 作矩阵-向量乘法 $\mathbf{X}=\mathbf{D}\mathbf{x}$。其中，离散傅里叶变换矩阵的元素具有规律性：$D_{ij}={\omega }^{ij}$，$0\le i,j\le N-1$（$\omega ={\mathrm{e}}^{-2\pi \mathrm{i}/N}$ 为复数方程 $z^N=1$ 的单位根，此处 $\mathrm{i}$ 为虚数单位，注意与行列下标 $i,j$ 区分）。

DFT 矩阵

快速傅里叶变换算法 FFT

由于离散傅里叶变换矩阵具有高度规律性，涉及该矩阵的乘法必然存在巧妙的计算方法以加速运算。快速傅里叶变换算法（FFT）正是这样一种加速 DFT 矩阵-向量乘法的算法。

FFT 算法的关键：奇偶分离（以 8 个数为例）

FFT 算法的原理十分简洁，仅需对矩阵乘法的列进行重新排列即可发现规律。以 8 元离散傅里叶变换为例，其原始形式如下：

8 元离散傅里叶变换

在矩阵乘法中，对右侧向量的行进行置换，同时对矩阵的列进行相应置换，结果保持不变。因此，可将偶数列（0, 2, 4, 6 列，程序员下标从 0 开始）移至左侧，对应地将偶数下标的 $x$（0, 2, 4, 6 行）移至上方，得到变换后的结果：

奇偶分离的 8 元离散傅里叶变换

观察上述变换后的分块矩阵，可见每一列的上下子矩阵之间均相差一个固定倍数。这并非偶然，正是 FFT 算法的关键所在。

奇偶分离的形式化描述

一般地（对于偶数 $N$），上述分块结果可表示为下图，其中 ${\mathbf{x}}_{\frac{N}{2}}$ 表示偶数行，${\mathbf{x}}_{\frac{N}{2}+1}$ 表示奇数行。

因 $\omega$ 为方程 $z^N=1$ 的单位根，故 ${\omega }^N=1$。为方便讨论，此处考虑 $N$ 为偶数的情形，于是 ${\omega }^{\frac{N}{2}}=-1$。由此可得四个分块矩阵的关系：

${\mathbf{D}}_2={\mathbf{D}}_1\cdot {\omega }^{Nj}={\mathbf{D}}_1$

${\mathbf{D}}_4={\mathbf{D}}_3\cdot {\omega }^{Nj+\frac{N}{2}(-N+1)}={\mathbf{D}}_3{\omega }^{\frac{N}{2}}=-{\mathbf{D}}_3$

借助以上关系，$N$ 个数的傅里叶变换结果可简化为：

因此，只需计算对偶数下标元素作用 ${\mathbf{D}}_1$、对奇数下标元素作用 ${\mathbf{D}}_3$ 所得结果 ${\mathbf{D}}_1{\mathbf{x}}_{\frac{N}{2}}$ 和 ${\mathbf{D}}_3{\mathbf{x}}_{\frac{N}{2}+1}$，即可得到 $\mathbf{D}\mathbf{x}$ 的完整结果。

离散傅里叶变换的递归结构

注意到 ${\mathbf{D}}_1={\omega }^{i\cdot (2j)}=({\omega }^2{)}^{ij}$ 即为 $\frac{N}{2}$ 个元素的 DFT 矩阵，故 ${\mathbf{D}}_1{\mathbf{x}}_{\frac{N}{2}}$ 是对下标为偶数的 $\frac{N}{2}$ 个元素进行离散傅里叶变换的结果。

注意到 ${\mathbf{D}}_3={\omega }^{i(2j-N+1)}=({\omega }^2{)}^{ij}\cdot {\omega }^i$，${\mathbf{D}}_3$ 与 ${\mathbf{D}}_1$ 仅相差一个行变换，因此 ${\mathbf{D}}_3{\mathbf{x}}_{\frac{N}{2}+1}$ 是对下标为奇数的 $\frac{N}{2}$ 个元素进行离散傅里叶变换后，再对第 $i$ 个元素乘以 ${\omega }^i$ 的结果。

由此得到离散傅里叶变换的递归结构：

$N$ 个元素的离散傅里叶变换结果可由奇数下标、偶数下标的 $\frac{N}{2}$ 个元素分别进行离散傅里叶变换的结果组合得到。

快速傅里叶变换的复杂度

记 $N$ 个元素的快速傅里叶变换复杂度为 $T(N)$，根据递归结构，可得递归公式：

$$T(N)=2T\left(\frac{N}{2}\right)+\mathcal{O}(N)$$

根据复杂度分析的主定理，立即可得快速傅里叶变换的时间复杂度为 $\mathcal{O}(N\log N)$。

离散傅里叶变换的应用

离散傅里叶变换矩阵常出现在一些特殊矩阵的分解中。

循环矩阵（Circulant Matrix）

循环矩阵由第一行经循环移位生成：

可以证明， C = D − 1 d i a g { D c } D \boldsymbol{C} = \boldsymbol{D}^{-1} \mathrm{diag}\{\boldsymbol{D}\boldsymbol{c}\} \boldsymbol{D} C=D−1diag{Dc}D，其中 $\mathbf{c}$ 为 $\mathbf{C}$ 的第一列。于是 C x = D − 1 d i a g { D c } D x \boldsymbol{C}\boldsymbol{x} = \boldsymbol{D}^{-1} \mathrm{diag}\{\boldsymbol{D}\boldsymbol{c}\} \boldsymbol{D}\boldsymbol{x} Cx=D−1diag{Dc}Dx，即 $\mathbf{D}\mathbf{c}$ 与 $\mathbf{D}\mathbf{x}$ 的逐元素乘法结果再做离散傅里叶逆变换 ${\mathbf{D}}^{-1}$。因此，涉及循环矩阵的矩阵-向量乘法亦可在 $\mathcal{O}(N\log N)$ 时间内完成。

Toeplitz 矩阵

Toeplitz 矩阵由第一行和第一列共同决定，其余元素沿对角线方向复制：

值得注意的是：大小为 $n\times n$ 的 Toeplitz 矩阵可嵌入大小为 $2n\times 2n$ 的循环矩阵中。下图的分块矩阵中，左上角和右下角为两个相同的 Toeplitz 矩阵。

因此，包含 Toeplitz 矩阵的矩阵-向量乘法可转化为循环矩阵的乘法，亦可在 $\mathcal{O}(N\log N)$ 时间内完成。

$$T\underline{v}=\left(\begin{array}{cc}{I}_n& 0_n\end{array}\right)A\left(\begin{array}{c}\underline{v}\\ 0_n\end{array}\right)$$

附加说明

• 为书写方便，本文仅讨论 $N$ 为偶数的情形。$N$ 为奇数的情形类似，读者可自行验证。
• 本文中 ${\mathbf{D}}_1={\omega }^{i\cdot (2j)}$ 的表示形式意指该矩阵的第 $i$ 行第 $j$ 列元素为 ${\omega }^{i\cdot (2j)}$，行列值从 0 开始计数。
• 本文参考了 MIT 课件 和 维基百科。部分文献中的离散傅里叶变换矩阵 $\mathbf{D}$ 包含归一化系数 $\frac{1}{\sqrt{N}}$ 或使用相反的辐角 $\omega ={\mathrm{e}}^{2\pi \mathrm{i}/N}$，这些区别不影响对 FFT 算法复杂度的分析。
• 复杂度分析的主定理参见 维基百科。
• 离散傅里叶变换的应用、循环矩阵和 Toeplitz 矩阵的分解参考了 MIT 课件。
• 循环矩阵在深度学习中的应用，可参见 ICCV 2015 论文《An Exploration of Parameter Redundancy in Deep Networks with Circulant Projections》。

---

如何理解和掌握快速傅里叶变换的计算和概念？

学理科的徐大喵 编辑于 2023-07-14 14:58 · 浙江

一、周期序列的离散傅里叶级数

首先给出周期性离散时间信号 $f_N(k)$ 的定义，下标 $N$ 表示其周期，即：
$$f_N(k)=f_N(k+lN),l=0,\pm 1,\pm 2,\dots$$

周期离散序列：序列 $f_N(k)$ 在整数 $k$ 上周期重复，周期为 $N$。即对于任意整数 $l$，序列值在 $k$ 与 $k+lN$ 处相等。这是离散傅里叶级数（DFS）和离散傅里叶变换（DFT）的基础假设。

根据前文所述，周期为 $T$ 的连续时间信号 $f(t)$ 可表示为一系列虚指数项 ${\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{T}t}$ 之和，称为傅里叶级数的复数形式：
$$\begin{array}{rl}f(t)& =\sum _{n=-\infty }^{+\infty }{c}_n{\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{T}t}\\ c_n& =\frac{1}{T}{\int }_{-T/2}^{T/2}f(t){\mathrm{e}}^{-\mathrm{j}n\frac{2\pi }{T}t}\mathrm{d}t\end{array}$$

类似地，这种傅里叶级数的分析方法亦可应用于离散周期序列。这是因为，离散周期序列 $f_N(k)$ 可由周期为 $T$ 的连续时间信号 $f(t)$ 经周期性采样获得（采样周期为 $T_s$）：

此时，得到的离散序列其周期可以表示为 $N=T/T_s$，因此周期序列 $f_N(k)$ 也可以表示为一系列虚指数项 ${\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{N}k}$ 和的形式：

$${\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{T}k{T}_s}={\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{N}k}$$

同时，根据复数运算，这组虚指数满足：

$${\mathrm{e}}^{\mathrm{j}(n+lN)\frac{2\pi }{N}k}={\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{N}k},l=0,\pm 1,\pm 2,\dots$$

因此，其中独立分布的虚指数仅有 $N$ 项，取 $n=0,1,\dots ,N-1$ 并写出 $f_N(k)$ 的展开式：

$$f_N(k)=\sum _{n=0}^{N-1}{c}_n{\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{N}k}$$

式中 $c_n$ 为待定系数。类比连续傅里叶级数求解待定系数的方法，上式两端同乘 ${\mathrm{e}}^{-\mathrm{j}m\frac{2\pi }{N}k}$ 并在一个周期内对 $k$ 求和（右端交换求和次序）：

$$\begin{array}{rl}\sum _{k=0}^{N-1}{f}_N(k){\mathrm{e}}^{-\mathrm{j}m\frac{2\pi }{N}k}& =\sum _{k=0}^{N-1}\left[{\mathrm{e}}^{-\mathrm{j}m\frac{2\pi }{N}k}\sum _{n=0}^{N-1}{c}_n{\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{N}k}\right]\\ & =\sum _{n=0}^{N-1}\left[c_n\sum _{k=0}^{N-1}{\mathrm{e}}^{\mathrm{j}(n-m)\frac{2\pi }{N}k}\right]\end{array}$$

根据等比数列求和公式，右端的求和项满足：

$$\sum _{k=0}^{N-1}{\mathrm{e}}^{\mathrm{j}(n-m)\frac{2\pi }{N}k}=\{\begin{array}{ll}N,& n=m\\ \frac{1-{\mathrm{e}}^{\mathrm{j}(n-m)\frac{2\pi }{N}\cdot N}}{1-{\mathrm{e}}^{\mathrm{j}(n-m)\frac{2\pi }{N}}}=0,& n,m=0,1,\dots ,N-1,n\ne m\end{array}$$

即当 $m=n$ 时，求得待定系数：

$$c_n=\frac{1}{N}\sum _{k=0}^{N-1}{f}_N(k){\mathrm{e}}^{-\mathrm{j}n\frac{2\pi }{N}k}$$

至此，离散周期序列可表示为一系列虚指数项之和的形式：

$$\begin{array}{rl}{f}_N(k)& =\sum _{n=0}^{N-1}{c}_n{\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{N}k}\\ c_n& =\frac{1}{N}\sum _{k=0}^{N-1}{f}_N(k){\mathrm{e}}^{-\mathrm{j}n\frac{2\pi }{N}k}\end{array}$$

给出一个实例，如下图所示为一离散周期序列 $N=4$，将其分解为离散傅里叶级数：

作出分解的结果示意图：

clear all
t = linspace(-2, 6, 1000);
f = 1/2 + 0.25*cos(t*pi/2) + 0.25*sin(t*pi/2) + 0.25*cos(3*t*pi/2) - 0.25*sin(3*t*pi/2);
plot(t, f, '--', 'LineWidth', 1.2)
grid on
hold on
k = [-1 0 1 2 3 4 5];
fk = [0 1 1 0 0 1 1];
stem(k, fk, '.', 'LineWidth', 1.2, 'MarkerSize', 20)
xlabel('k')
ylabel('f_N')

分解结果符合原周期序列。

二、快速傅里叶变换（FFT）

（a）离散傅里叶变换（DFT）

离散傅里叶变换主要用于计算机实现对信息的处理。由于计算机只能处理有限长度的离散数据，而周期序列的离散傅里叶级数恰好满足这些特征：只需对计算机输入在信号的一个周期 $T$ 内采样得到的 $N$ 个数据点（采样周期为 $T_s$），通过计算即可得到以虚指数项 ${\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{N}k}$ 表示的 $N$ 个谐波分量：

$$\begin{array}{rl}{f}_N(k)& =\sum _{n=0}^{N-1}{c}_n{\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{N}k}\\ c_n& =\frac{1}{N}\sum _{k=0}^{N-1}{f}_N(k){\mathrm{e}}^{-\mathrm{j}n\frac{2\pi }{N}k}\end{array}$$

根据各虚指数项的系数 $c_n$，定义周期序列的离散傅里叶变换为：

$$F_N(n)=\sum _{k=0}^{N-1}{f}_N(k){\mathrm{e}}^{-\mathrm{j}n\frac{2\pi }{N}k},n=0,1,\dots ,N-1$$
为还原该周期序列，其逆变换为：

$$f_N(k)=\frac{1}{N}\sum _{n=0}^{N-1}{F}_N(n)\cdot e^{jn\frac{2\pi }{N}k}$$

（b）MATLAB 应用实例

给出一连续时间信号：

$$f(t)=\frac{1}{2}+2\sin (10\omega t)+\sin (20\omega t)$$

其中，角频率 $\omega =100\pi$。取采样频率为基波频率的 100 倍，对其进行离散化采样得到周期序列 $f_N(k)$，使用离散傅里叶变换进行分析：

clear all
N = 100;        % 100 个采样点
w = 2*pi*50;    % 基波 50 Hz
T = 1/50;       % 基波周期
Ts = T/N;       % 采样周期
kN = 0:(N-1);
% 原连续时间信号
% ft = 1/2 + 2*sin(10*w*t) + sin(20*w*t);
% 采样获得有限长度离散数据
fN = 1/2 + 2*sin(10*w*kN*Ts) + sin(20*w*kN*Ts);
% 定义离散傅里叶变换
Fn = zeros(1, N);
for n = 0:N-1
    for k = 0:N-1
        Fn(n+1) = Fn(n+1) + fN(k+1)*exp(-1j*n*2*pi*k/N);
    end
end
A = abs(Fn);
plot(kN, A, 'LineWidth', 1.5)
grid on
xlabel('n')
ylabel('F_n')

作出离散傅里叶变换所得的幅度谱：

观察该幅度谱，可以发现其具有以下特点：

1） 在各次谐波分量中，除直流分量外，幅度谱关于 $N/2$ 对称。观察这一系列虚指数：

即其中关于 $N/2$ 对称的虚指数是一对共轭复数，它们幅值相等、相位相反。因此根据复数运算法则，关于 $N/2$ 对称的 $F(n)$ 也是一对共轭复数，于是可知幅度谱关于 $N/2$ 对称。

在此基础上可推得：关于 $N/2$ 对称的项 $F(n){\mathrm{e}}^{\mathrm{j}n\frac{2\pi }{N}k}$ 也是一对共轭复数，这意味着它们的实部相等。在离散傅里叶级数中，当等式左端 $f_N(k)$ 为实数序列时，等式右端复数求和项中的虚部会相消，仅保留实部（可参考第二部分实例）。因此关于 $N/2$ 对称的虚指数项可看作同一次谐波，即离散傅里叶变换最高可分解出 $(N/2-1)$ 次谐波（计算中向上取整）。

在本节实例中，$N=100$ 则最高可分解出 49 次谐波；若 $N=101$ 则最高可分解出 50 次谐波。

2） 原连续时间信号中的直流分量 $A_0$、基波分量和各次谐波分量的幅值 $A(n)$ 与离散傅里叶变换各项幅值之间的关系满足：

其中，系数 $1/N$ 由离散傅里叶变换的定义显然可知；系数 $2$ 与上一条中将关于 $N/2$ 对称的虚指数项看作同一次谐波相对应。

因此由离散傅里叶变换得到原连续时间信号各次谐波的幅值：

clear all
N = 100;        % 100 个采样点
w = 2*pi*50;    % 基波 50 Hz
T = 1/50;       % 基波周期
Ts = T/N;       % 采样周期
f = 50;
kN = 0:(N-1);
% 原连续时间信号
% ft = 1/2 + 2*sin(10*w*t) + sin(20*w*t);
% 采样获得有限长度离散数据
fN = 1/2 + 2*sin(10*w*kN*Ts) + cos(20*w*kN*Ts);
% 定义离散傅里叶变换
Fn = zeros(1, N);
for n = 0:N-1
    for k = 0:N-1
        Fn(n+1) = Fn(n+1) + fN(k+1)*exp(-1j*n*2*pi*k/N);
    end
end
A = abs(Fn);
if mod(N, 2) == 0
    fp = (0:N/2-1)*f;
    Ap(1) = A(1)/N;
    Ap(2:N/2-1+1) = 2*A(2:N/2-1+1)/N;         % MATLAB 数组标号起始为 1
else
    fp = (0:(N+1)/2-1)*f;
    Ap(1) = A(1)/N;
    Ap(2:(N+1)/2-1+1) = 2*A(2:(N+1)/2-1+1)/N; % MATLAB 数组标号起始为 1
end
plot(fp, Ap, 'LineWidth', 1.5)
grid on
xlabel('f/Hz')
ylabel('A_n')

（c）采样定理

连续时间信号的基波周期为 $T$，基波频率为 ${\omega }_T=2\pi /T$；

离散化采样周期为 $T_s$，采样频率为 ${\omega }_s=2\pi /T_s$；

其中，离散序列的周期 $N=T/T_s$，即 $N={\omega }_s/{\omega }_T$。

根据前述离散傅里叶变换幅度谱的特点（1），在此条件下能够分解出 $0\sim (N/2-1)$ 次谐波分量。定义连续时间信号中的最高次谐波频率为 ${\omega }_m$，那么要保证傅里叶变换能够分解出原信号的所有谐波分量，必须满足：

由于 $N$ 只能取整数，不等式亦可写为：

此即符合采样定理。

（d）快速傅里叶变换（FFT）

快速傅里叶变换（FFT）利用了虚指数项的对称性等特点，对离散傅里叶变换（DFT）的计算实现了简化，从而提高计算机的求解速度，其计算结果与离散傅里叶变换（DFT）完全一致。本节介绍一种快速傅里叶变换算法，适用于序列周期为 2 的整数次幂（$N=2^l$）的情况——

首先将离散傅里叶变换中的奇数项和偶数项分离（假设 $N=2^l$）：

上式中，将原周期序列 $f_N(k)$ 中的奇数项和偶数项分别看作以 $N/2$ 为周期的新序列 $f_{N/2}^{\prime }$ 和 $f_{N/2}^{\prime \prime }$，那么当 $n$ 的取值范围为 $0\sim N/2-1$ 时，上式中的两个求和项显然满足离散傅里叶变换的定义式，即：

此时，就将 $N$ 个数据点的离散傅里叶变换（$N$ 点 DFT）转化为了两个 $N/2$ 个数据点的离散傅里叶变换（$N/2$ 点 DFT）。但考虑到 $n$ 的取值范围，显然只计算出了"$N$ 点 DFT"的前半部分，还需找到后半部分的值。根据虚指数的对称性：

“$N$ 点 DFT” 后半部分的值为：

因此，"$N$ 点 DFT"的计算结果可用奇偶分离后两个"$N/2$ 点 DFT"的计算结果来表示。以 $N=8$ 为例，画出示意图（图中 $W_N^n={\mathrm{e}}^{\mathrm{j}n2\pi /N}$）：

对于"$N/2$ 点 DFT"还可对该序列继续进行奇偶分离，将其简化为"$N/4$ 点 DFT"：

继续进行奇偶分离，直至分解为"1 点 DFT"：

根据定义，原"$N$ 点 DFT"的计算量为：$N^2$ 次复数乘法运算和 $N(N-1)$ 次加法运算：

而参考前述示意图，当 $N=2^l$ 时，奇偶分离后的快速傅里叶变换（$N$ 点 FFT）包含 $l$ 级，每级进行 $N/2$ 次复数乘法和 $N$ 次加法运算，显然计算得到了简化。

奇偶分离的过程天然适用于递归算法，给出递归算法实现快速傅里叶变换（FFT）的 MATLAB 程序：

function Fn = myDFTx(x)
% x 为输入序列
N = length(x);      % 序列长度
L = 2^nextpow2(N);
if L == N           % 序列长度为 2 的幂次，进入 FFT 算法
    if N == 1
        Fn = x;
    else
        n = 0:N/2-1;
        wn = exp(-1j*2*pi/N);
        Wn = wn.^(n');
        x_even = x(1:2:end); % 提取偶数项
        x_odd = x(2:2:end);  % 提取奇数项
        Xe = myDFT2(x_even);
        Xo = myDFT2(x_odd);
        Fn = [Xe + Wn.*Xo; Xe - Wn.*Xo];
    end
else                % 序列长度非 2 的幂次，进入 DFT 算法
    Fn = zeros(1, N);
    for n = 0:N-1
        for k = 0:N-1
            Fn(n+1) = Fn(n+1) + x(k+1)*exp(-1j*n*2*pi*k/N);
        end
    end
end

注：本节所述快速傅里叶算法假设了数据点数量满足 $N=2^l$，因此在数据点数量不满足条件时，程序中仍采用定义法求解离散傅里叶变换。为加快求解速度，当采样点数量不满足 $N=2^l$ 时，可对序列补零后使用快速傅里叶算法。补零后，快速傅里叶变换（FFT）计算出的结果存在偏差，更加优化的操作还请读者继续查阅资料。

三、总结

在 Simulink 仿真中，我们常使用 FFT 分析工具，截取一个基波周期的波形，在数据处理后能够得到直流分量以及各次谐波分量相对于基波分量的大小：

亦可选择将波形数据导入 MATLAB 工作区，在数组中截取一个基波周期的 $N$ 个数据点，使用自带的 fft 函数进行分析。如前所述：

1. fft 计算结果中 $n$ 表示谐波的次数，最多可分解出 $N/2-1$ 次及以下谐波；
2. 计算结果的幅值与实际直流分量以及各次谐波幅值之间的关系满足如下关系式：

---

参考文献

• 快速理解 FFT 算法（完整无废话）. 知乎, 2021. https://zhuanlan.zhihu.com/p/407885496
• 如何理解和掌握快速傅里叶变换的计算和概念？. 知乎, 2014. https://www.zhihu.com/question/20456490/answer/25440654
  • 快速傅里叶变换（FFT）基本原理与应用实例 - 知乎
    https://zhuanlan.zhihu.com/p/624717583
  • 周期信号的单边频谱与双边频谱 - 知乎
    https://zhuanlan.zhihu.com/p/564444855
  • 线性连续系统的离散化实例（冲激响应不变法） - 知乎
    https://zhuanlan.zhihu.com/p/578006816