快速傅里叶变换（Fast Fourier Transform, FFT）是离散傅里叶变换（Discrete Fourier Transform, DFT）的一种高效计算方法。它并非一种全新的变换，而是DFT在计算效率上的一次革命性突破。理解FFT的原理，就是理解它如何巧妙地重组计算步骤，大幅减少运算量，从而使频域分析在各种工程和科学领域得以广泛应用。
FFT是什么？它的核心概念与DFT的关系
DFT的本质与FFT的诞生
要理解FFT，首先要从DFT说起。DFT是一种将离散时域信号转换为离散频域信号的数学工具。对于一个包含N个样本的时域序列$x[n]$（$n=0, 1, …, N-1$），其DFT定义为：
$X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j \frac{2\pi}{N} nk}$
其中，$k=0, 1, …, N-1$。
这个公式看起来简洁，但其计算量却相当庞大。直观来看，计算每一个频点$X[k]$都需要进行N次乘法和N-1次加法。总共有N个频点需要计算，因此总的计算复杂度大约是$O(N^2)$。当N变得非常大时（例如，几十万甚至上百万个样本），这种$O(N^2)$的计算量会使得DFT变得不切实际，耗时巨大。
FFT的核心：加速DFT
FFT的核心概念在于它是一种计算DFT的算法，通过“分治”策略将一个大规模的DFT问题分解成多个小规模的DFT问题，并巧妙地利用DFT的周期性和对称性来消除冗余计算，从而极大地提高了计算效率。它计算的仍然是与DFT相同的结果，但以指数级更快的速度完成。因此，可以这样理解：FFT是DFT的“快车道”。
FFT与DFT的关系
- 结果一致： FFT算法计算出的结果与直接按照DFT定义式计算出的结果是完全相同的（在数值精度允许的范围内）。
- 本质不变： 它们都完成了时域到频域的转换，都揭示了信号中包含的频率成分。
- 方法不同： DFT是“做什么”（变换），FFT是“怎么做”（算法）。FFT是DFT的一种高效实现方法。
为什么FFT会比DFT快？性能提升的根本原因是什么？
计算复杂度的巨大差异
这是FFT最引人注目的优势。对于N点序列的DFT，其计算复杂度是$O(N^2)$。而对于FFT，特别是基2（radix-2）FFT，其计算复杂度是$O(N \log_2 N)$。
为了直观感受这个差异，我们来看几个例子：
- 当N=1024时：
- DFT复杂度：$1024^2 \approx 100万$
- FFT复杂度：$1024 \times \log_2(1024) = 1024 \times 10 \approx 1万$
此时FFT比DFT快了大约100倍。
- 当N=65536（$2^{16}$）时：
- DFT复杂度：$65536^2 \approx 42亿$
- FFT复杂度：$65536 \times \log_2(65536) = 65536 \times 16 \approx 100万$
此时FFT比DFT快了大约4200倍。
随着N的增大，FFT的性能优势呈指数级增长，使得对大数据量进行频域分析成为可能。
性能提升的根本原因：分治策略与旋转因子周期性
FFT之所以能实现如此巨大的性能提升，主要归功于以下两个核心思想：
-
分治策略（Divide and Conquer）：
FFT算法将一个N点DFT分解为两个N/2点DFT，然后将这两个N/2点DFT的输出组合起来得到N点DFT的输出。这个分解过程可以递归进行，直到分解为2点或1点DFT。2点DFT的计算非常简单，1点DFT则直接等于输入值。通过这种分解，原本需要计算一个大矩阵乘法的问题，被拆解成了多个小矩阵乘法的组合。
-
旋转因子（Twiddle Factor）的周期性和对称性：
DFT公式中的$e^{-j \frac{2\pi}{N} nk}$被称为旋转因子，通常简写为$W_N^{nk}$。这个旋转因子具有非常重要的周期性和对称性，例如$W_N^{k+N} = W_N^k$ 和 $W_N^{k+N/2} = -W_N^k$。FFT算法巧妙地利用这些性质，使得在合并子DFT结果时，许多乘法计算可以被复用或避免，特别是通过“蝶形运算”结构，将重复的计算合并处理，极大地减少了乘法和加法的总次数。
这两点结合起来，构成了FFT算法的精髓，实现了从$O(N^2)$到$O(N \log N)$的复杂度飞跃。
FFT的基本工作原理：如何将大DFT分解成小DFT？
Cooley-Tukey算法：基2分解法
最常用的FFT算法是Cooley-Tukey算法，特别是其基2（radix-2）版本。这个算法要求输入的样本点数N必须是2的幂次方（即$N=2^M$，其中M为整数）。如果N不是2的幂次方，通常需要进行零填充（zero-padding）来凑成2的幂次方。
分解策略：偶数项与奇数项
Cooley-Tukey算法的核心思想是将N点DFT分解为两个N/2点DFT，分别处理原始序列的偶数索引项和奇数索引项。假设$X[k]$是N点DFT的结果，我们可以将求和项分为偶数n和奇数n：
$X[k] = \sum_{n \text{ even}} x[n] \cdot W_N^{nk} + \sum_{n \text{ odd}} x[n] \cdot W_N^{nk}$
设$n=2m$（偶数项）和$n=2m+1$（奇数项），其中$m=0, 1, …, N/2-1$。代入上式得到：
$X[k] = \sum_{m=0}^{N/2-1} x[2m] \cdot W_N^{(2m)k} + \sum_{m=0}^{N/2-1} x[2m+1] \cdot W_N^{(2m+1)k}$
利用旋转因子的性质$W_N^{2mk} = (W_N^2)^{mk} = W_{N/2}^{mk}$，我们可以重写上式：
$X[k] = \sum_{m=0}^{N/2-1} x[2m] \cdot W_{N/2}^{mk} + W_N^k \sum_{m=0}^{N/2-1} x[2m+1] \cdot W_{N/2}^{mk}$
令 $G[k] = \sum_{m=0}^{N/2-1} x[2m] \cdot W_{N/2}^{mk}$ 为偶数索引序列的N/2点DFT，
令 $H[k] = \sum_{m=0}^{N/2-1} x[2m+1] \cdot W_{N/2}^{mk}$ 为奇数索引序列的N/2点DFT。
那么，N点DFT的$X[k]$可以表示为：
$X[k] = G[k] + W_N^k H[k]$
这对于$k=0, 1, …, N/2-1$是成立的。但DFT的输出有N个点，即$k$的范围是$0, …, N-1$。对于$k \geq N/2$，我们利用$W_N^{k+N/2} = -W_N^k$ 以及$G[k+N/2]=G[k]$、$H[k+N/2]=H[k]$（因为$G[k]$和$H[k]$是N/2点DFT，它们以N/2为周期），可以得到：
$X[k+N/2] = G[k] – W_N^k H[k]$
这就是FFT的“蝶形运算”公式，它将两个N/2点DFT的结果$G[k]$和$H[k]$组合起来，通过一次乘法（与$W_N^k$相乘）和两次加减法，得到两个N点DFT的结果$X[k]$和$X[k+N/2]$。
蝶形运算（Butterfly Operation）
蝶形运算是FFT算法的核心计算单元。它得名于其计算图示形似蝴蝶。一个基本的蝶形运算接受两个输入（例如$A$和$B \cdot W_N^k$），然后产生两个输出（$A+B \cdot W_N^k$和$A-B \cdot W_N^k$）。
具体到Cooley-Tukey算法中：
- 输入是$G[k]$和$H[k]$。
- 将$H[k]$与旋转因子$W_N^k$相乘，得到$W_N^k H[k]$。
- 第一个输出是$X[k] = G[k] + W_N^k H[k]$。
- 第二个输出是$X[k+N/2] = G[k] – W_N^k H[k]$。
通过这种方式，原本需要独立计算的$X[k]$和$X[k+N/2]$可以利用相同的$G[k]$、$H[k]$和$W_N^k H[k]$，以极少的额外计算完成。这便是效率提升的关键。
位反转（Bit-Reversal Permutation）
在基2 FFT算法中，为了实现原地计算（in-place computation），输入序列通常需要进行“位反转”重排。这是因为在Cooley-Tukey分解的过程中，序列的偶数项和奇数项被不断地分离。当分解到最底层时，各个子序列的原始索引会变得混乱。位反转操作就是为了在算法开始前，将输入序列按照最终计算所需的数据流顺序进行预排列，使得在蝶形运算的每个阶段，所需的数据都能相邻或以固定步长访问。
例如，对于N=8（000到111），位反转如下：
- 000 -> 000 (0)
- 001 -> 100 (4)
- 010 -> 010 (2)
- 011 -> 110 (6)
- 100 -> 001 (1)
- 101 -> 101 (5)
- 110 -> 011 (3)
- 111 -> 111 (7)
所以，原始序列$[x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7]$会变成$[x_0, x_4, x_2, x_6, x_1, x_5, x_3, x_7]$作为FFT的输入。
迭代过程与阶段
一个N点（$N=2^M$）的FFT算法通常包含M个阶段（stages）。每个阶段都会执行N/2个蝶形运算。随着阶段的推进，蝶形运算的“跨度”会逐渐增大。例如，在第一个阶段，蝶形运算处理的是相距N/2的元素；在最后一个阶段，蝶形运算处理的是相邻元素。
FFT的计算发生在哪个阶段？N点FFT需要多少次乘法和加法？
FFT计算的发生阶段
FFT的计算通常发生在需要从时域信号提取频域特征的任何环节。它本身是一个算法，一个计算过程。在实际应用中，FFT是信号处理链条中的一个关键“模块”或“步骤”。
例如，在一个实时音频分析系统中：
- 采样： 模拟音频信号被模数转换器（ADC）转换为离散的数字样本序列（时域）。
- 分帧与加窗： 为了处理连续信号，通常会将其分割成短小的“帧”，并对每帧应用窗函数（如汉明窗），以减少频谱泄漏。
- FFT计算： 对加窗后的每一帧数据执行FFT算法，将时域样本转换为频域表示。
- 频域分析： 对FFT的输出结果进行进一步处理，如计算频谱幅度、功率谱密度、提取频率特征（如基频、谐波）、滤波等。
- 应用： 将频域分析结果用于均衡器调整、语音识别、音乐信息检索、振动分析等各种应用。
可见，FFT计算紧随在信号的采样和预处理之后，是时域到频域转换的核心步骤。
N点FFT的乘法和加法次数
对于一个基2的N点FFT（其中$N=2^M$），其计算量大致如下：
- 乘法次数： $ \frac{N}{2} \log_2 N $ 次复数乘法。
- 例如，对于$N=1024$，需要 $ \frac{1024}{2} \times 10 = 5120 $ 次复数乘法。
- 加法次数： $ N \log_2 N $ 次复数加法。
- 例如，对于$N=1024$，需要 $ 1024 \times 10 = 10240 $ 次复数加法。
需要注意的是，一个复数乘法通常需要4次实数乘法和2次实数加法。一个复数加法则需要2次实数加法。因此，总的实数运算量会更大。但与DFT的$O(N^2)$实数乘法和加法相比，FFT的$O(N \log N)$仍然是巨大的性能提升。
如何理解FFT的输出结果？（频率分辨率、Nyquist采样定理）
FFT输出的频域表示
FFT的输出是一个包含N个复数值的序列$X[0], X[1], …, X[N-1]$。每一个复数值$X[k]$都代表了原始时域信号中某个特定频率成分的幅度和相位信息。
- 幅度： $ |X[k]| $ 表示该频率成分的强度。
- 相位： $ \text{angle}(X[k]) $ 表示该频率成分相对于参考相位的偏移。
通常我们更关心幅度谱，即$|X[k]|$的分布。
频率分辨率（Frequency Resolution）
频率分辨率是FFT分析结果中一个非常重要的概念，它表示FFT能够区分的最小频率间隔。对于N点FFT，如果采样率为$F_s$，则频率分辨率$ \Delta f $为：
$ \Delta f = \frac{F_s}{N} $
这意味着FFT的输出频率点是 $0, \Delta f, 2\Delta f, …, (N-1)\Delta f$。频率分辨率越高（$\Delta f$越小），我们就能越精确地分辨出彼此接近的频率成分。要提高频率分辨率，可以通过增加采样点数N（例如，通过采集更长时间的数据或进行零填充）来实现。
Nyquist采样定理与频谱折叠（Aliasing）
奈奎斯特采样定理指出，为了无失真地从离散样本中重建原始模拟信号，采样频率$F_s$必须至少是信号中最高频率成分的两倍，即$F_s \ge 2 F_{max}$。这个最高频率$F_{max}$被称为奈奎斯特频率$F_{Nyquist} = F_s/2$。
FFT的输出中，前N/2个点（从$X[0]$到$X[N/2-1]$）通常代表正频率，从$0$到$F_s/2 – \Delta f$。而$X[N/2]$则代表$F_s/2$频率成分。对于实数输入信号，FFT输出的后半部分（从$X[N/2+1]$到$X[N-1]$）是前半部分的共轭对称，代表负频率成分，或者也可以解释为折叠到正频率范围的高频部分。
如果原始模拟信号中存在频率高于奈奎斯特频率$F_s/2$的成分，这些高频成分在采样后会“折叠”到$0$到$F_s/2$的频率范围内，导致低频部分的频谱失真，这种现象称为频谱混叠（Aliasing）。为了避免混叠，通常在进行ADC采样之前，会使用低通滤波器（称为抗混叠滤波器）滤除高于$F_s/2$的频率成分。
幅度和相位信息
-
DC分量 ($X[0]$)：
对应于零频率，即信号的直流（DC）分量或平均值。对于实数信号，它是一个实数。
-
幅度谱：
通常，我们绘制$|X[k]|$的图形，横轴是频率$k \cdot \Delta f$，纵轴是幅度。这直接显示了信号在不同频率上的能量分布。
需要注意的是，对于实数信号，由于频谱的对称性，我们通常只分析从$0$到$F_s/2$的频率范围（即$X[0]$到$X[N/2]$）。此时，除了$X[0]$和$X[N/2]$，其他频率成分的实际幅度需要乘以2，因为它们的能量被分摊到了正负频率两个对称分量中。
-
相位谱：
绘制$\text{angle}(X[k])$的图形，横轴是频率，纵轴是相位。相位谱揭示了不同频率成分在时间轴上的相对起始位置。在许多应用中，相位信息与幅度信息同样重要，例如在通信系统的解调或图像处理中。
如何选择合适的FFT点数N？
N的选择对分辨率和计算量的影响
选择合适的FFT点数N是进行有效频域分析的关键决策。
-
对频率分辨率的影响：
正如前面提到的，$ \Delta f = F_s/N $。更大的N意味着更小的$\Delta f$，即更高的频率分辨率。这允许我们更精确地区分信号中相邻的频率成分。
例如，如果有一个信号包含两个非常接近的频率（如100Hz和101Hz），要分辨它们，就需要$\Delta f \le 1$Hz。如果采样率是1000Hz，那么N至少需要1000（即$1000/1000=1$Hz）。
-
对计算量的影响：
FFT的计算复杂度是$O(N \log N)$。N越大，计算时间越长。在实时系统中，必须在分辨率要求和计算资源（处理器速度、内存）之间取得平衡。
-
对时间窗长的影响：
对于给定采样率$F_s$，N点FFT对应的时间窗长度（即采集数据的时间）为$T = N/F_s$。要提高频率分辨率，通常意味着需要采集更长时间的数据。
零填充（Zero-padding）的作用与副作用
当N不是2的幂次方时，或者当需要更精细的频谱“视图”但没有更多数据可供分析时，可以使用零填充。
作用：
- 凑成2的幂次方： 如果原始数据点数M不是2的幂次方，在数据末尾添加零直到N达到最接近且大于M的2的幂次方，以便使用高效的基2 FFT算法。
- 改善频谱的视觉效果（“插值”）： 零填充不会增加原始信号的实际频率分辨率，因为它没有引入新的信息。但它会在原始DFT的频率点之间“插值”出更多的频率点，使得频谱曲线看起来更平滑、更连续。这有助于更好地观察频谱峰值的位置和形状，尤其是在相邻频率之间。
副作用：
- 不增加实际信息： 零填充并不能真正提高系统分辨不同频率的能力（即真实频率分辨率$ \Delta f_{true} = F_s/M $，其中M是原始数据点数），它只是在现有信息的基础上增加了频谱的显示点数。
- 可能引入伪影： 如果原始信号没有经过适当的窗函数处理，零填充可能导致频谱中出现额外的旁瓣，这可能掩盖真实的信号特征。
因此，零填充是一个有用的工具，但必须理解其局限性：它提高了频谱的“显示分辨率”，而不是“物理分辨率”。真正的物理分辨率取决于原始数据的时间长度。
总结来说，FFT原理的核心在于其分治的算法思想和对旋转因子对称性的巧妙利用，这使得它能够以远超DFT的速度完成时域到频域的转换，从而成为现代数字信号处理的基石。在实际应用中，对FFT点数N的选择、奈奎斯特采样定理的遵循以及对频率分辨率的理解，都是正确解读和利用FFT结果的关键。
