4、FFT基本原理：离散傅里叶变换到FFT、蝶形运算、旋转因子、按时间抽取与按频率抽取
在CMSIS-DSP定点FFT的实现与优化中，理解FFT（快速傅里叶变换）的数学本质与运算结构是掌握缩放策略与溢出控制的前提。FFT并非一种新的变换，而是离散傅里叶变换（DFT）的一种高效计算算法。本章将从DFT的定义出发，剖析其计算瓶颈，进而引出FFT的核心思想——蝶形运算与旋转因子，并详细对比按时间抽取（DIT）与按频率抽取（DIF）两种经典分解路径。这些原理直接决定了定点实现中数据位宽的增长模式、缩放因子的插入位置以及溢出风险的分布特征。
上图展示了DFT直接计算与FFT分解的核心差异，以及旋转因子的数学特性。DFT的O(N²)复杂度在嵌入式实时信号处理中是不可接受的，而FFT通过分治策略将复杂度降至O(N·log₂N)。对于CMSIS-DSP定点实现而言，理解蝶形运算的数值行为比单纯记忆公式更为重要——因为定点数的有限字长效应会在每一级蝶形中累积，而缩放策略正是基于这种累积规律设计的。
4.1 离散傅里叶变换的计算瓶颈与FFT的诞生
4.1.1 核心原理：DFT的数学定义与计算量分析
对于长度为N的有限长序列x[n]，其N点DFT定义为：
X[k] = Σ_{n=0}^{N-1} x[n] · W_N^{nk}, k = 0, 1, ..., N-1其中W_N = e^{-j2π/N}为旋转因子。直接计算每个X[k]需要N次复数乘法和(N-1)次复数加法，计算全部N个频点则需要N²次复数乘法和N(N-1)次复数加法。当N=1024时，复数乘法次数超过100万次，这在Cortex-M系列处理器上即使以CMSIS-DSP优化库运行，也需要数毫秒量级，无法满足实时频谱分析需求。
FFT的核心洞察在于：旋转因子W_N^{nk}具有周期性和对称性，使得DFT运算中存在大量重复计算。通过将N点DFT分解为多个较短的DFT，可以显著减少乘法次数。以基-2 FFT为例，N点DFT可分解为两个N/2点DFT，每个N/2点DFT再分解为两个N/4点DFT，如此递归直至2点DFT。整个分解过程形成log₂N级蝶形运算，每级包含N/2个蝶形单元。
4.1.2 判断方法：何时必须使用FFT而非直接DFT
在实际工程中，选择FFT还是直接DFT取决于以下因素：
- 变换点数N：当N≥32时，FFT的优势已非常明显；N<16时直接DFT可能更简单。
- 实时性要求：若系统要求频谱更新率高于1kHz，N=256以上的DFT几乎不可能实时完成，必须使用FFT。
- 定点实现约束：CMSIS-DSP提供定点FFT函数（如arm_cfft_q15），其内部已针对Cortex-M内核的SIMD指令和饱和运算进行了优化，直接调用库函数比手写DFT更可靠。
4.1.3 案例拆解：从DFT到FFT的复杂度对比
假设一个电机振动监测系统需要以40kHz采样率采集1024点数据进行频谱分析。直接DFT需要约100万次复数乘法，在180MHz的Cortex-M4上每次复数乘法约需10个周期，总耗时约55ms。而基-2 FFT只需约5120次复数乘法（1024×log₂1024/2），耗时约0.28ms，性能提升近200倍。这一差距在定点实现中更为关键——DFT的累加次数多，定点溢出风险远高于FFT。
4.1.4 常见误区：忽略FFT的位宽增长特性
初学者常误以为FFT减少了运算量，也同时减少了定点溢出的风险。实际上，FFT的每级蝶形运算都会产生数据位宽增长。以Q15格式为例，输入数据为16位有符号数，经过一级蝶形后，加法结果可能达到17位，乘法结果可能达到31位。CMSIS-DSP的q15 FFT函数内部采用饱和运算和可选的缩放因子来管理这种增长，但用户必须理解缩放因子的作用，否则输出可能严重失真。
4.2 蝶形运算：FFT的基本计算单元
4.2.1 核心原理：蝶形运算的数学结构与数据流
基-2蝶形运算是FFT的最小计算单元，其输入为两个复数a和b，输出为两个复数A和B：
A = a + b·W_N^k
B = a - b·W_N^k其中W_N^k为旋转因子。从计算角度看，每个蝶形需要1次复数乘法（b×W_N^k）和2次复数加法/减法。在定点实现中，复数乘法涉及4次实数乘法和2次实数加减法，而复数加减法各涉及2次实数加减法。CMSIS-DSP的蝶形运算针对Cortex-M内核的饱和算术指令进行了优化，能够在单个周期内完成16位饱和加减法。
蝶形运算的关键特性是原位计算（in-place）：输出A和B可以直接覆盖输入a和b的存储位置，无需额外缓冲区。这一特性使得FFT的内存占用仅为输入数据的2倍（实部+虚部），非常适合内存受限的嵌入式系统。
4.2.2 操作流程：蝶形运算的定点实现步骤
在Q15定点格式下，蝶形运算的实现流程如下：
- 加载操作数：从内存读取a_re, a_im, b_re, b_im（均为q15格式）。
- 旋转因子查表：从预计算的旋转因子表中读取W_N^k的实部（cos值）和虚部（-sin值），同样为q15格式。
- 复数乘法：计算temp_re = b_re×cos - b_im×(-sin)，temp_im = b_re×(-sin) + b_im×cos。乘法结果需进行饱和左移（通常左移1位以保留符号位）。
- 复数加减：A_re = a_re + temp_re，A_im = a_im + temp_im；B_re = a_re - temp_re，B_im = a_im - temp_im。加减法使用饱和指令防止溢出。
- 存储结果：将A和B写回原地址。
4.2.3 案例拆解：8点FFT的蝶形运算全流程
以8点基-2 FFT为例，共需3级蝶形运算（log₂8=3），每级4个蝶形。第一级蝶形使用旋转因子W_8^0=1，第二级使用W_8^0和W_8^2，第三级使用W_8^0、W_8^1、W_8^2、W_8^3。下表展示了各级蝶形的旋转因子分布：
| 级数 | 蝶形索引 | 旋转因子 | 旋转因子值（Q15近似） |
|---|---|---|---|
| 第1级 | 0,1,2,3 | W_8^0 | 32767 + j0 |
| 第2级 | 0,1 | W_8^0 | 32767 + j0 |
| 第2级 | 2,3 | W_8^2 | 0 - j32767 |
| 第3级 | 0 | W_8^0 | 32767 + j0 |
| 第3级 | 1 | W_8^1 | 23170 - j23170 |
| 第3级 | 2 | W_8^2 | 0 - j32767 |
| 第3级 | 3 | W_8^3 | -23170 - j23170 |
在定点实现中，旋转因子表通常以const数组形式存储在Flash中，避免运行时计算三角函数。CMSIS-DSP提供了arm_cfft_q15函数，内部自动处理旋转因子查表和蝶形运算，用户只需提供输入数组和变换长度。
4.2.4 常见误区：忽视蝶形运算中的饱和处理
许多开发者在定点FFT实现中只关注算法逻辑，却忽略了饱和处理的重要性。当输入信号接近满量程（如Q15格式的±32767）时，蝶形运算的加法结果可能超过16位范围。CMSIS-DSP的饱和指令（如SSAT）会自动将结果钳位到最大/最小值，但这会引入非线性失真。正确的做法是在FFT前对输入数据进行适当的缩放（如右移2位），确保整个变换过程中数据不会频繁触及饱和边界。
4.3 旋转因子：FFT的相位旋转核心
4.3.1 核心原理：旋转因子的数学性质与查表优化
旋转因子W_N^k = e^{-j2πk/N} = cos(2πk/N) - j·sin(2πk/N)是FFT运算中唯一的乘法系数。其核心性质包括：
- 周期性：W_N^{k+N} = W_N^k，因此只需存储k=0到N/2-1的N/2个值。
- 对称性：W_N^{k+N/2} = -W_N^k，W_N^{N/2-k} = -j·W_N^k（共轭对称）。
- 可约性：W_{2N}^{2k} = W_N^k，用于混合基FFT。
在定点实现中，旋转因子以Q15格式存储。例如W_8^1 = cos(π/4) - j·sin(π/4) ≈ 0.7071 - j0.7071，Q15近似值为23170 - j23170。CMSIS-DSP的旋转因子表按特定顺序排列，以支持基-2、基-4和混合基FFT算法。
4.3.2 判断方法：旋转因子表的存储策略
选择旋转因子表的存储方式需考虑以下因素：
- Flash容量：N点FFT需要N/2个复数旋转因子，每个复数占4字节（2×q15），N=1024时需2KB Flash。对于资源紧张的MCU，可只存储1/4周期（N/4个值），其余通过对称性推导。
- 计算速度：查表比实时计算三角函数快约10-20倍。CMSIS-DSP默认使用全表存储，以换取最高执行速度。
- 精度需求：Q15格式的旋转因子精度约为0.00003，对于大多数信号处理应用足够。若需要更高精度，可使用Q31格式，但运算速度会下降。
4.3.3 案例拆解：旋转因子在CMSIS-DSP中的实际应用
在CMSIS-DSP的arm_cfft_q15函数中，旋转因子表由arm_cfft_q15_init函数初始化。用户需要提供旋转因子表数组和变换长度。例如，初始化1024点FFT：
arm_cfft_instance_q15 fft_inst;
arm_cfft_q15_init(&fft_inst, 1024, 0, 1);
// 参数: 实例指针, FFT长度, 正变换/反变换标志, 位反转标志该函数内部会生成或引用预定义的旋转因子表。对于q15格式，CMSIS-DSP使用16位旋转因子，并在蝶形运算中通过饱和乘法指令（如SMULBB）实现复数乘法。值得注意的是，CMSIS-DSP的旋转因子表是按位反转顺序排列的，这是为了配合DIT算法的输入位反转要求。
4.3.4 常见误区：旋转因子精度不足导致的频谱泄漏
当旋转因子以Q15格式存储时，其量化误差约为±1 LSB。对于N=1024的FFT，这种误差在各级蝶形中累积，可能导致输出信噪比下降约6dB。如果应用对频谱精度要求较高（如高精度谐波分析），建议使用Q31格式的旋转因子或浮点FFT。CMSIS-DSP同时提供arm_cfft_f32浮点版本，但执行速度约为定点版本的1/4。
4.4 按时间抽取（DIT）与按频率抽取（DIF）
4.4.1 核心原理：两种分解策略的数学等价性
DIT和DIF是FFT的两种基本分解策略，它们在数学上等价，但数据流方向不同：
- DIT（按时间抽取）：将输入序列x[n]按奇偶位置分解为两个N/2点子序列，分别进行DFT，再通过蝶形运算合并。输入需要位反转排列，输出为自然顺序。
- DIF（按频率抽取）：将输出序列X[k]按奇偶频率分解，输入为自然顺序，输出需要位反转排列。蝶形运算中先进行加减法，再进行旋转因子乘法。
两种方法的计算量完全相同（N/2·log₂N次复数乘法），但蝶形运算的结构不同。DIT的蝶形为“先乘后加”，DIF的蝶形为“先加后乘”。这一差异在定点实现中影响显著：DIT的乘法结果直接参与加法，可能产生较大的中间值；DIF的加法结果先被缩放，再进行乘法，数值范围更可控。
4.4.2 操作流程：DIT与DIF的定点实现对比
在CMSIS-DSP中，arm_cfft_q15函数默认使用DIT算法（基-4混合基）。其内部流程为：
- 位反转排列：将输入数据按位反转顺序重新排列。
- 基-4蝶形运算：每级处理4个输入，产生4个输出，使用3个旋转因子。
- 缩放控制：每级蝶形后执行一次右移2位的操作（相当于乘以0.25），防止溢出。
对于DIF算法，CMSIS-DSP未直接提供，但用户可通过修改旋转因子表顺序和蝶形结构自行实现。DIF的优势在于输入数据无需位反转，适合数据流式处理场景（如实时频谱分析仪）。
4.4.3 案例拆解：DIT与DIF在定点溢出控制中的差异
考虑一个输入信号幅度接近满量程的Q15数据序列。在DIT算法中，第一级蝶形的乘法结果（b×W_N^k）可能达到Q30格式，随后与a相加时需先右移15位，这可能导致精度损失。而在DIF算法中，第一级蝶形先执行a+b和a-b，结果范围扩大1倍（Q16格式），随后乘以旋转因子时需右移16位，精度损失更小。因此，对于输入信号动态范围较大的场景，DIF算法在定点实现中通常具有更好的数值稳定性。
| 特性 | DIT（按时间抽取） | DIF（按频率抽取） |
|---|---|---|
| 输入顺序 | 位反转 | 自然顺序 |
| 输出顺序 | 自然顺序 | 位反转 |
| 蝶形结构 | 先乘后加 | 先加后乘 |
| 定点溢出风险 | 乘法后加法可能溢出 | 加法后乘法范围更可控 |
| CMSIS-DSP支持 | arm_cfft_q15（默认） | 需自定义实现 |
| 适用场景 | 通用频谱分析 | 流式处理、高动态范围输入 |
4.4.4 常见误区：混淆位反转顺序与缩放策略
一个常见错误是认为DIT和DIF的位反转顺序可以互换。实际上，位反转顺序是算法结构决定的，不能随意更改。如果使用DIT算法但输入数据未进行位反转，输出将完全错误。CMSIS-DSP的arm_cfft_q15函数内部已处理位反转，用户只需确保输入数据按自然顺序排列即可。另一个误区是忽视缩放策略的差异——DIT和DIF的每级缩放因子可能不同，直接套用会导致溢出或精度下降。