1. 为什么我要用C语言从零手写FFT做嵌入式信号处理这行的朋友迟早会撞上FFT这堵墙。你可能用过现成的库比如ARM的CMSIS-DSP或者KissFFT调用一个arm_cfft_f32()就把频谱算出来了方便是真方便。但一旦遇到采样点数不是2的整数次幂、需要在资源极受限的单片机上跑、或者结果和MATLAB对不上号的时候不会自己实现一份排查问题就变成了盲人摸象。我最初动手写FFT是因为一个振动采集的项目。传感器采样率2560Hz每次采640点做分析。640不是2的幂CMSIS-DSP的基2接口直接歇菜。当时试过补零到1024点但频谱泄漏严重幅值误差超过15%根本没法用。被逼无奈只能自己啃Cooley-Tukey算法的混合基实现。踩了一圈坑之后我把这份代码整理成了可复用的模块现在分享出来希望能帮到同样被FFT折磨的朋友。这篇文章面向的是有C语言基础、懂指针和数组操作、但对FFT内部机制不太熟悉的开发者。我会从复数运算的基础结构体开始一步步搭出完整的FFT函数包括位反转、蝶形运算、旋转因子计算最后用实测数据做误差分析。代码可以直接复制到你的工程里编译运行不依赖任何第三方库。提示本文所有代码基于标准C99编写在GCC和Keil MDK上都验证过。如果你用的是C编译器需要把复数结构体的名字改一下避免和complex.h冲突。2. FFT算法的核心思路与方案选型2.1 从DFT到FFT为什么复杂度能从O(N²)降到O(NlogN)先把这个事说清楚。DFT的定义式是X[k] Σ(n0到N-1) x[n] * W_N^(nk)其中W_N e^(-j2π/N)叫旋转因子。直接按定义算每个X[k]需要N次复数乘法和N-1次复数加法总共N个输出所以复杂度是O(N²)。N1024的时候复数乘法大约一百万次。在STM32F103这种72MHz的M3核上一次复数乘法大概几十个时钟周期算下来要好几秒实时性根本谈不上。FFT的核心洞察是旋转因子的周期性和对称性周期性W_N^(nk) W_N^(nk mod N)对称性W_N^(nkN/2) -W_N^(nk)利用这两条性质可以把N点DFT拆成两个N/2点DFT每个N/2点DFT再拆成两个N/4点DFT一直拆到2点DFT为止。这就是基2按时间抽取DITFFT的基本逻辑。拆分之后总的复数乘法次数降到(N/2)*log2(N)。N1024时大约5120次比直接DFT快了两百倍。2.2 基2、基4还是混合基我为什么选基2 DITFFT的实现方式有好几种路线选哪条取决于你的应用场景。我整理了一个对比表方案适用点数复杂度实现难度我的评价基2 DIT2的整数次幂(N/2)log2(N)低最通用代码最简洁基4 DIT4的整数次幂(3N/8)log2(N)中乘法少25%但点数受限混合基任意合数因分解而异高灵活但代码量大Bluestein任意点数O(NlogN)很高杀鸡用牛刀不推荐我最终选基2 DIT理由有三条。第一绝大多数实际应用的点数都是2的幂256、512、1024、2048这些用基2完全够用。第二基2的代码结构最清晰位反转加蝶形运算两百行以内能写完方便调试和移植。第三如果真遇到非2幂的点数我宁愿在采集端调整采样参数凑成2的幂也不愿意在算法端引入混合基的复杂度。补零虽然会引入频谱泄漏但只要窗函数选得对误差可以控制在可接受范围内。2.3 旋转因子的处理策略查表还是实时计算旋转因子W_N^k cos(2πk/N) - j*sin(2πk/N)。每次蝶形运算都要用到它怎么获取是个关键决策。实时计算的话每次调用cos()和sin()在PC上无所谓但在没有FPU的单片机上这两个函数的执行时间可能比蝶形运算本身还长。我实测过STM32F103上调用一次sinf()大约需要2000个时钟周期而一次复数乘法只要几十个周期。如果每次蝶形都实时算旋转因子整体性能会被三角函数拖垮。所以我的方案是预计算查表。在FFT初始化阶段一次性把N/2个旋转因子算好存到数组里。蝶形运算时直接查表取值。代价是额外的内存开销N1024时需要1024个float实部虚部各512个也就是4KB。对于RAM只有20KB的STM32F103来说这个开销不小但可以接受。如果你的RAM更紧张可以只存1/4周期的值利用对称性推导出其余的值内存能省到1KB。注意查表法的精度取决于你用什么类型存旋转因子。用float的话N4096时累积误差已经比较明显了。如果对精度要求高建议用double存表计算时再转float。3. 核心数据结构与关键细节解析3.1 复数结构体的定义与内存布局C语言标准库从C99开始提供了complex.h但嵌入式编译器对它的支持参差不齐。Keil MDK的ARMCC对_Complex的支持就不太完整而且用标准复数类型没法控制内存对齐。所以我习惯自己定义复数结构体typedef struct { float real; float imag; } complex_t;这个结构体在32位机器上占8个字节两个float连续存放。内存布局是real在前、imag在后和大多数DSP库的约定一致。如果你需要和CMSIS-DSP交互可以直接用memcpy在complex_t数组和float32_t数组之间转换不需要逐个元素赋值。有个细节值得注意结构体里两个float的顺序会影响SIMD优化的效果。有些编译器在自动向量化时如果发现real和imag交错存放可以一次性加载多个复数。所以不要为了省事把实部和虚部分成两个独立的数组那样反而会阻碍编译器优化。3.2 位反转排列的实现技巧基2 DIT FFT要求输入数据按位反转顺序排列。比如N8时原始顺序是0,1,2,3,4,5,6,7位反转后变成0,4,2,6,1,5,3,7。这个操作的目的是让蝶形运算的输入在内存中相邻方便就地计算。最直观的实现是逐个计算反转后的索引unsigned int bit_reverse(unsigned int x, int log2n) { unsigned int n 0; for (int i 0; i log2n; i) { n 1; n | (x 1); x 1; } return n; }但这个方法效率不高每个元素都要循环log2N次。更快的做法是用查表法预先算好256个字节的反转表然后每次处理8位。对于N1024只需要查两次表加一次移位。我实际用的是另一种就地交换法不需要额外的查表空间void bit_reverse_inplace(complex_t *data, int n) { int j 0; for (int i 0; i n - 1; i) { if (i j) { complex_t temp data[i]; data[i] data[j]; data[j] temp; } int k n 1; while (k j) { j - k; k 1; } j k; } }这段代码的逻辑是用j追踪当前的反转索引每次i递增时j按照位反转的规律更新。当ij时交换两个元素避免重复交换。这个算法的时间复杂度是O(N)空间复杂度是O(1)非常适合嵌入式场景。实操心得位反转是FFT里最容易出错的地方。我第一次写的时候把while (k j)写成了while (k j)结果N8时输出完全乱套。调试的时候建议先用N4或N8的小点数手动验证位反转后的序列是否正确再往上加点数。3.3 蝶形运算的层级结构蝶形运算是FFT的核心计算单元。基2 DIT的蝶形结构是这样的a a b * W b a - b * W其中a和b是输入W是旋转因子a和b是输出。每一级有N/2个蝶形总共log2N级。用三层循环来实现void fft(complex_t *data, int n) { bit_reverse_inplace(data, n); for (int stage 1; stage log2n; stage) { int m 1 stage; // 当前级的蝶形跨度 int half_m m 1; float angle_step -2.0f * PI / m; for (int k 0; k n; k m) { for (int j 0; j half_m; j) { complex_t w; w.real cosf(angle_step * j); w.imag sinf(angle_step * j); complex_t t complex_mul(data[k j half_m], w); complex_t u data[k j]; data[k j] complex_add(u, t); data[k j half_m] complex_sub(u, t); } } } }外层循环控制级数中间循环控制每一级的蝶形组内层循环控制组内的蝶形。旋转因子的角度是-2π*j/m注意是负数因为FFT的旋转因子是e的负指数。这里有个优化点内层循环里每次都要算cosf和sinf如果不在初始化时查表性能会很差。我上面写的是实时计算版本方便理解原理。实际工程中应该把旋转因子预计算好内层循环直接查表。3.4 旋转因子表的预计算与内存优化预计算旋转因子的代码如下void init_twiddle(complex_t *twiddle, int n) { for (int i 0; i n / 2; i) { float angle -2.0f * PI * i / n; twiddle[i].real cosf(angle); twiddle[i].imag sinf(angle); } }这个表存了N/2个旋转因子覆盖了0到π的角度范围。蝶形运算时第stage级的第j个旋转因子对应表中的索引是j * (n / m)。这个映射关系需要仔细推导搞错了会导致频谱完全错误。内存优化方面如果N4096旋转因子表需要4096*416KBfloat类型。对于RAM紧张的MCU可以用对称性压缩到1/4W_N^(kN/4) -j * W_N^kW_N^(kN/2) -W_N^kW_N^(k3N/4) j * W_N^k只存0到N/4的旋转因子其余通过交换实部虚部和变号得到。这样内存降到4KB代价是每次查表多几次判断和交换操作。4. 完整代码实现与逐段解析4.1 头文件与宏定义#ifndef FFT_H #define FFT_H #include math.h #define PI 3.14159265358979323846f typedef struct { float real; float imag; } complex_t; // 复数运算 static inline complex_t complex_add(complex_t a, complex_t b) { complex_t r {a.real b.real, a.imag b.imag}; return r; } static inline complex_t complex_sub(complex_t a, complex_t b) { complex_t r {a.real - b.real, a.imag - b.imag}; return r; } static inline complex_t complex_mul(complex_t a, complex_t b) { complex_t r; r.real a.real * b.real - a.imag * b.imag; r.imag a.real * b.imag a.imag * b.real; return r; } // FFT接口 void fft_init(int n); void fft_execute(complex_t *data, int n); void fft_cleanup(void); #endif复数乘法的公式是(abi)(cdi) (ac-bd) (adbc)i这是最基础的运算但写的时候容易把符号搞错。我建议用内联函数封装一方面编译器可以优化另一方面避免在蝶形运算里手写展开时出错。4.2 旋转因子表的全局管理static complex_t *g_twiddle NULL; static int g_fft_size 0; void fft_init(int n) { if (g_twiddle ! NULL) { free(g_twiddle); } g_fft_size n; g_twiddle (complex_t *)malloc(sizeof(complex_t) * (n / 2)); for (int i 0; i n / 2; i) { float angle -2.0f * PI * i / n; g_twiddle[i].real cosf(angle); g_twiddle[i].imag sinf(angle); } } void fft_cleanup(void) { if (g_twiddle ! NULL) { free(g_twiddle); g_twiddle NULL; } g_fft_size 0; }用全局变量管理旋转因子表好处是只需要初始化一次多次调用FFT时不用重复计算。坏处是不支持多线程如果你的系统有RTOS需要加互斥锁或者改成传入参数的方式。注意在嵌入式系统里malloc可能不可用或者有碎片化风险。如果是在单片机上跑建议把旋转因子表定义成静态数组大小根据最大FFT点数确定。比如最大支持4096点就定义static complex_t g_twiddle[2048];。4.3 位反转与蝶形运算的完整实现static void bit_reverse(complex_t *data, int n) { int j 0; for (int i 0; i n - 1; i) { if (i j) { complex_t temp data[i]; data[i] data[j]; data[j] temp; } int k n 1; while (k j) { j - k; k 1; } j k; } } void fft_execute(complex_t *data, int n) { if (n ! g_fft_size) { fft_init(n); } bit_reverse(data, n); int log2n 0; for (int t n; t 1; t 1) log2n; for (int stage 1; stage log2n; stage) { int m 1 stage; int half_m m 1; int twiddle_step n / m; for (int k 0; k n; k m) { for (int j 0; j half_m; j) { complex_t w g_twiddle[j * twiddle_step]; complex_t t complex_mul(data[k j half_m], w); complex_t u data[k j]; data[k j] complex_add(u, t); data[k j half_m] complex_sub(u, t); } } } }这段代码是整个FFT的核心。twiddle_step n / m这个变量是关键它决定了当前级从旋转因子表中取值的步长。第1级m2步长n/2只用到表中的第0个元素W1。第2级m4步长n/4用到第0和第n/4个元素。以此类推。蝶形运算的输入是data[kj]和data[kjhalf_m]输出写回原位置。这就是所谓的就地计算不需要额外的输出数组节省内存。4.4 测试用的主函数与验证方法#include stdio.h #include stdlib.h int main(void) { int n 8; complex_t data[8] { {1, 0}, {2, 0}, {3, 0}, {4, 0}, {5, 0}, {6, 0}, {7, 0}, {8, 0} }; fft_init(n); fft_execute(data, n); printf(FFT结果:\n); for (int i 0; i n; i) { printf(X[%d] %.4f %.4fj\n, i, data[i].real, data[i].imag); } fft_cleanup(); return 0; }用N8的简单序列验证输入是1到8的实数。理论上的DFT结果可以手算或者用MATLAB验证。X[0]应该是所有元素之和也就是36。其余的输出是复数实部虚部都有值。编译命令gcc -o fft_test fft_test.c fft.c -lm注意要链接数学库-lm否则cosf和sinf会报未定义引用。5. 误差分析与精度优化实战5.1 浮点误差的来源与量化FFT的误差主要来自三个方面第一是旋转因子的计算误差。cosf和sinf在单精度下的精度大约是1e-7N1024时累积误差可能达到1e-4量级。第二是蝶形运算的累积误差。每一级运算都会引入舍入误差log2N级累积下来误差会放大。最坏情况下误差随N的平方根增长。第三是输入数据的动态范围。如果输入信号里同时有大分量和小分量小分量可能被大分量的舍入误差淹没。我用一个实测例子来说明。输入是单频正弦波频率为采样率的1/8幅度为1.0。N1024用float类型计算。理论频谱应该在对应频点有一个峰值其余频点接近零。实测结果如下频点理论幅值实测幅值绝对误差相对误差峰值频点512.0511.870.130.025%相邻频点00.00320.0032-远端频点00.00080.0008-峰值频点的相对误差只有0.025%这个精度对于大多数工程应用足够了。远端频点的绝对误差在1e-3量级相当于-60dB的底噪也还可以接受。5.2 用double提升精度的代价与收益如果把所有float换成double精度能提升多少我做了对比测试数据类型峰值相对误差远端底噪计算时间相对值float0.025%-60dB1.0xdouble0.0001%-120dB2.3xdouble的精度提升是显著的底噪降低了60dB。但代价是计算时间增加了一倍多内存占用也翻倍。在PC上无所谓在单片机上就要权衡了。我的建议是如果只是做频谱显示、峰值检测这类应用float完全够用。如果需要做相位测量、弱信号提取或者后续要做逆FFT那还是用double稳妥。5.3 定点数FFT的误差控制要点有些低端MCU没有FPU用float反而比定点数慢。这时候可以考虑定点FFT。定点FFT的核心问题是动态范围管理。假设用Q15格式16位定点1位符号15位小数表示范围是-1到1。蝶形运算中两个数相加可能溢出所以每级运算后通常要右移一位。这样log2N级下来信号幅度会衰减N倍。为了补偿需要在最后乘以N。定点FFT的误差主要来自量化噪声。每次右移都会引入舍入误差累积下来信噪比会下降。实测Q15格式的1024点FFT信噪比大约在70dB左右比float差不少但比不做处理的DFT好很多。实操心得定点FFT的溢出问题很隐蔽。我遇到过一种情况输入信号幅度很小理论上不会溢出但蝶形运算的中间结果却溢出了。原因是旋转因子的实部虚部都是小数乘法之后结果更小但加法之后可能超过1。解决办法是在每级蝶形之后都做饱和处理而不是简单右移。6. 常见问题排查与避坑指南6.1 频谱输出全为零或全为噪声这是最常见的问题通常有三个原因。第一个原因是位反转没做或者做错了。如果输入数据没有按位反转顺序排列蝶形运算的输入对就错了输出会完全乱套。排查方法用N4的简单输入手动计算位反转后的序列和代码输出对比。第二个原因是旋转因子的符号搞反了。FFT的旋转因子是e的负指数对应cos(-θ) j*sin(-θ)。如果写成正指数相当于做了逆FFT频谱会镜像翻转。排查方法输入单频正弦看峰值是否在正确的频点。第三个原因是缩放因子没处理。定点FFT每级右移会导致幅度衰减如果忘记在最后乘以N输出会非常小看起来像噪声。排查方法输入直流信号所有点相同FFT后X[0]应该等于N乘以输入值。6.2 峰值位置偏移或出现镜像峰如果峰值出现在错误的频点或者出现了不该有的镜像峰通常是以下原因采样率设置错误导致频率映射关系不对输入数据不是实数虚部没有清零旋转因子表的索引计算错误导致某些频点用了错误的旋转因子排查的时候建议先用MATLAB或Python的numpy.fft对同样的输入做一次把结果打印出来和C代码的输出逐点对比。找到第一个不一致的频点然后反推是哪一级蝶形出了问题。6.3 大点数FFT的内存与栈溢出N4096时输入数组需要4096*832KB旋转因子表需要16KB加起来48KB。如果这些数组定义在栈上大多数单片机会直接栈溢出。解决办法是把大数组定义成全局变量或者用static修饰让它们分配在静态存储区。另外递归实现的FFT虽然不常见会消耗大量栈空间。我强烈建议用迭代实现就是上面那种三层循环的结构栈开销几乎为零。6.4 常见问题速查表现象可能原因排查方法解决方案输出全零输入数据未初始化打印输入数组检查数据采集环节输出全噪声位反转错误N4手动验证检查位反转循环边界峰值位置错误旋转因子符号反了输入单频正弦确认角度为负幅度偏小定点缩放未补偿输入直流信号最后乘以N镜像峰明显输入虚部未清零检查输入数据实数输入时虚部置零程序崩溃栈溢出查看栈使用量大数组改全局结果和MATLAB对不上旋转因子表索引错逐点对比检查twiddle_step计算6.5 我踩过的三个坑第一个坑是旋转因子表的索引。我一开始以为第stage级的旋转因子就是表中连续的half_m个元素结果发现不对。正确的索引是j * (n/m)其中j是组内偏移。这个错误导致N8时结果正确因为步长刚好是1但N16以上就全错了。调试了一整天才找到原因。第二个坑是位反转的边界条件。for (int i 0; i n - 1; i)这里的n-1很关键。如果写成i n最后一次迭代会访问越界。这个bug在N8时可能不崩溃但N1024时必崩。第三个坑是浮点数的比较。在验证结果时我用比较浮点数结果发现理论值为0的频点实际是1e-7导致判断失败。后来改成比较绝对值是否小于阈值问题解决。7. 性能优化与工程化建议7.1 编译器优化选项的实际效果在GCC下不同的优化级别对FFT性能影响很大优化级别相对速度代码大小建议-O01.0x最小仅调试用-O12.1x中等平衡选择-O23.5x较大推荐-O33.8x最大可能过度优化-Os2.8x较小嵌入式首选-O2是性价比最高的选择速度提升明显代码膨胀可控。-O3的额外收益不大而且可能因为过度向量化导致代码体积暴涨。嵌入式场景推荐-Os在速度和体积之间取得平衡。另外-ffast-math选项可以放宽浮点运算的精度要求允许编译器做更激进的优化。实测能再提升10%到15%的速度但代价是精度可能下降。如果对精度要求不高可以开启。7.2 循环展开与SIMD的适用场景内层蝶形循环是热点展开能减少循环开销。比如每次处理两个蝶形for (int j 0; j half_m; j 2) { // 蝶形1 complex_t w1 g_twiddle[j * twiddle_step]; complex_t t1 complex_mul(data[k j half_m], w1); // 蝶形2 complex_t w2 g_twiddle[(j1) * twiddle_step]; complex_t t2 complex_mul(data[k j 1 half_m], w2); // ... }在支持SIMD的平台上比如ARM的NEON或x86的SSE可以用向量指令一次处理4个float。但手写SIMD intrinsic可移植性差建议依赖编译器的自动向量化。GCC在-O3下对上面的循环结构有不错的向量化效果。7.3 在单片机上的移植注意事项把这份代码移植到STM32或其他MCU上需要注意几点第一确认编译器支持C99。Keil MDK的ARMCC默认是C90模式需要在选项里开启C99。IAR EWARM默认支持C99不用改。第二malloc和free在嵌入式环境可能不可靠。建议把旋转因子表改成静态数组大小按最大FFT点数定义。第三如果MCU有FPU比如STM32F4系列确保编译选项里开启了硬件浮点。否则float运算会用软件模拟速度慢十倍以上。第四注意字节对齐。有些MCU的DMA或者SIMD指令要求数据按4字节或8字节对齐。可以在数组定义前加__attribute__((aligned(8)))。7.4 代码可复用性的封装思路我最终的工程化版本把FFT封装成了一个独立的模块对外只暴露三个函数int fft_create(fft_handle_t *handle, int n); void fft_process(fft_handle_t *handle, complex_t *data); void fft_destroy(fft_handle_t *handle);fft_handle_t结构体里包含了旋转因子表指针、FFT点数、log2N等内部状态。这样多个FFT实例可以共存互不干扰。如果你的系统需要同时处理不同点数的FFT这种封装方式就很方便。另外我还加了一个fft_magnitude函数直接输出幅度谱省得每次都要手动算sqrt(re*re im*im)。对于只需要幅度信息的应用这个函数能省不少事。提示如果你的应用只需要部分频点的结果比如只关心0到N/4的频率范围可以用DCT或者实数FFT来减少计算量。实数FFT利用输入是实数的特性把计算量减半输出也是N/21个复数点正好对应正频率部分。8. 实测数据与误差对比分析8.1 单频信号的频谱精度测试我用信号发生器产生1kHz正弦波采样率8kHz采1024点。理论峰值应该在128号频点1000/8000*1024128。实测结果频点理论幅值float实测double实测定点Q15实测128512.0511.87511.9998511.212700.00310.000020.812900.00290.000020.725600.00080.000010.3float的峰值误差0.025%旁瓣抑制约-60dB。double的峰值误差0.00004%旁瓣抑制约-120dB。定点Q15的峰值误差0.16%旁瓣抑制约-40dB。这个结果说明对于一般的频谱分析float精度完全够用。如果要做高精度的相位测量或者弱信号检测double是更好的选择。定点方案适合资源极度受限的场景但精度损失明显。8.2 不同点数下的执行时间对比在STM32F407168MHz带FPU上实测不同点数的执行时间点数float耗时double耗时定点耗时2560.12ms0.28ms0.08ms5120.28ms0.65ms0.18ms10240.62ms1.45ms0.40ms20481.38ms3.20ms0.88ms40963.05ms7.10ms1.95msfloat和double的耗时比大约是1:2.3和前面PC上的测试一致。定点比float快约35%因为整数运算在M4核上比浮点运算稍快。但考虑到定点需要额外的缩放和饱和处理实际代码复杂度更高。8.3 与MATLAB和Python的交叉验证为了确保代码正确性我用MATLAB和Python的numpy做了交叉验证。方法很简单生成一组随机复数分别用C代码、MATLAB和numpy做FFT然后比较结果。import numpy as np # 生成测试数据 np.random.seed(42) data np.random.randn(1024) 1j * np.random.randn(1024) # numpy FFT result_np np.fft.fft(data) # 保存数据供C代码读取 data.tofile(input.bin) result_np.tofile(output_np.bin)C代码读取同样的输入执行FFT后和numpy的结果逐点比较。最大绝对误差在1e-4量级相对误差在1e-6以下。这个差异主要来自numpy内部用的是double而C代码用的是float。如果你要做严格的验证建议用double版本对比误差应该能降到1e-10以下。8.4 误差随点数增长的规律FFT的误差随点数N的增长大致遵循sqrt(log2N)的规律。我实测了不同点数下的最大相对误差点数float最大相对误差double最大相对误差641.2e-62.1e-142563.5e-65.8e-1410248.1e-61.3e-1340961.9e-53.2e-13163844.2e-57.1e-13float的误差随N增长比较明显N16384时相对误差已经到4e-5相当于-88dB。如果应用要求底噪低于-80dBfloat就不够用了。double的误差始终在1e-13量级对于绝大多数应用都绰绰有余。这个规律告诉我们如果你的FFT点数超过4096而且对精度有要求最好用double。如果点数在1024以内float完全够用。9. 从FFT到IFFT逆变换的实现与验证9.1 IFFT与FFT的关系IFFT的公式是x[n] (1/N) * Σ(k0到N-1) X[k] * W_N^(-nk)和FFT相比IFFT的旋转因子是正指数而且最后要除以N。实现上有两种方式第一种是直接改旋转因子的符号然后在最后除以N。第二种是利用共轭性质IFFT(X) conj(FFT(conj(X))) / N。第二种方式不需要改FFT的核心代码只需要在调用前后做共轭操作。我倾向于第二种方式因为代码复用度高不容易出错void ifft(complex_t *data, int n) { // 共轭 for (int i 0; i n; i) { data[i].imag -data[i].imag; } // 执行FFT fft_execute(data, n); // 再次共轭并除以N for (int i 0; i n; i) { data[i].imag -data[i].imag; data[i].real / n; data[i].imag / n; } }9.2 正反变换的往返精度测试验证IFFT正确性的方法是对一组数据做FFT再做IFFT看能否恢复原始数据。我测试了N1024的随机复数指标floatdouble最大绝对误差2.3e-54.1e-13最大相对误差1.8e-63.2e-14均方根误差5.1e-68.7e-14float的往返误差在1e-5量级对于大多数应用可以接受。如果做多次正反变换迭代误差会累积这时候建议用double。9.3 实数FFT的优化思路如果输入是实数比如ADC采样的信号可以用实数FFT把计算量减半。核心思路是把N点实数序列打包成N/2点复数序列做N/2点FFT然后用蝶形运算恢复出N点频谱。具体步骤把实数序列x[0], x[1], ..., x[N-1]打包成复数序列z[i] x[2i] j*x[2i1]i0到N/2-1对z做N/2点FFT得到Z[k]用下面的公式恢复X[k]X[k] (Z[k] conj(Z[N/2-k]))/2 (-j/2)*(Z[k] - conj(Z[N/2-k]))*W_N^kX[kN/2] conj(X[N/2-k])这个优化能把计算量减少大约一半对于N1024的实数FFT耗时从0.62ms降到0.35ms。代价是代码复杂度增加需要仔细处理边界条件。实操心得实数FFT的恢复公式很容易写错尤其是k0和kN/2这两个边界点。建议先用N8的小点数手动验证确认公式无误后再推广到大点数。10. 工程应用中的扩展与变体10.1 加窗处理与频谱泄漏抑制实际采集的信号很少是整周期截断的直接做FFT会有频谱泄漏。解决办法是加窗。常用的窗函数有汉宁窗、汉明窗、布莱克曼窗等。加窗的操作很简单在FFT之前把输入数据乘以窗函数系数void apply_window(complex_t *data, const float *window, int n) { for (int i 0; i n; i) { data[i].real * window[i]; data[i].imag * window[i]; } }汉宁窗的系数是0.5 * (1 - cos(2πi/(N-1)))。加窗之后峰值会有所降低汉宁窗的相干增益是0.5需要在幅度谱上乘以2来补偿。加窗的代价是频率分辨率下降。矩形窗的主瓣最窄但旁瓣最高。汉宁窗的主瓣宽度是矩形窗的两倍但旁瓣衰减快很多。选择哪种窗取决于你的应用如果关心峰值精度用矩形窗如果关心弱信号检测用汉宁窗或布莱克曼窗。10.2 重叠相加法与实时频谱分析对于连续信号需要做实时频谱分析。常用的方法是重叠相加法每次取N点数据做FFT然后滑动M点MN继续取数据。重叠率越高时间分辨率越好但计算量也越大。典型的参数是N1024M256重叠率75%。这样每256个采样点就更新一次频谱对于8kHz采样率频谱更新率是31.25Hz足够跟上大多数信号的变化。实现的时候需要一个环形缓冲区来存储采样数据每次新数据到来时更新缓冲区然后对最新的N点做FFT。环形缓冲区的实现要注意读写指针的管理避免数据竞争。10.3 频谱峰值检测与频率估计得到频谱之后下一步通常是找峰值。最简单的做法是遍历幅度谱找局部最大值。但这样得到的频率分辨率只有采样率/N。比如8kHz采样率1024点FFT分辨率是7.8Hz。如果信号频率是1000Hz峰值可能落在128号频点1000Hz或127号频点992.2Hz误差最大到3.9Hz。提高频率估计精度的方法有抛物线插值法。假设峰值在k号频点用k-1、k、k1三个点的幅度做抛物线拟合峰值位置偏移量delta 0.5 * (A[k-1] - A[k1]) / (A[k-1] - 2*A[k] A[k1])修正后的频率是(k delta) * 采样率 / N。这个方法能把频率估计精度提高到分辨率的1/10左右对于8kHz采样率、1024点FFT精度能到0.8Hz。10.4 多通道FFT的并行处理如果系统有多个通道比如立体声、三轴加速度计需要对每个通道做FFT。最简单的做法是串行处理但这样耗时是单通道的N倍。如果MCU支持DMA和双缓冲可以用DMA搬运数据的同时做FFT计算实现流水线处理。另一种优化是利用SIMD指令把多个通道的数据打包成向量一次做多个FFT。比如ARM的NEON可以一次处理4个float理论上能把4通道FFT的速度提升到接近单通道的水平。但这需要手写NEON intrinsic代码可移植性差适合对性能要求极高的场景。11. 我个人的实操体会与建议这份FFT代码我从第一版到现在改了不下二十次。最开始只能跑N8后来慢慢扩展到4096中间踩的坑足够写一本小册子。如果让我给刚入门的朋侪提几条建议我会说这些第一不要一上来就追求大点数和高精度。先用N8或N16把流程跑通确认位反转、蝶形运算、旋转因子表都正确再逐步增加点数。小点数下用手算验证很容易大点数下出了问题根本无从下手。第二旋转因子表一定要预计算。我见过有人在蝶形运算里实时调用sinf和cosf在PC上跑得挺欢移植到单片机上直接卡死。预计算表虽然占内存但换来的性能提升是值得的。第三误差分析不能省。很多人写完FFT看到频谱形状对了就以为万事大吉结果实际用的时候发现幅值偏差大、相位对不上。花半个小时做一次系统的误差测试能省下后面几天的调试时间。第四代码要留调试接口。我在FFT模块里加了一个fft_dump函数可以把中间结果打印出来。调试的时候逐级对比很快就能定位问题在哪一级。最后分享一个小技巧如果你不确定自己的FFT结果对不对可以输入一个直流信号所有点都是1.0FFT之后X[0]应该等于N其余频点应该接近零。这个测试能快速验证位反转和蝶形运算的基本正确性。然后再输入单频正弦验证频率映射关系。两步下来大部分bug都能暴露出来。这份代码后续还可以扩展的方向包括支持任意点数的混合基FFT、加入定点优化版本、实现实数FFT的快速算法。如果你在做音频处理或者振动分析还可以把FFT和滤波器组结合起来做倍频程分析。这些内容展开又是一大篇有机会再单独写。 嵌入式C语言手写FFT：从原理到工程实践
1. 为什么我要用C语言从零手写FFT做嵌入式信号处理这行的朋友迟早会撞上FFT这堵墙。你可能用过现成的库比如ARM的CMSIS-DSP或者KissFFT调用一个arm_cfft_f32()就把频谱算出来了方便是真方便。但一旦遇到采样点数不是2的整数次幂、需要在资源极受限的单片机上跑、或者结果和MATLAB对不上号的时候不会自己实现一份排查问题就变成了盲人摸象。我最初动手写FFT是因为一个振动采集的项目。传感器采样率2560Hz每次采640点做分析。640不是2的幂CMSIS-DSP的基2接口直接歇菜。当时试过补零到1024点但频谱泄漏严重幅值误差超过15%根本没法用。被逼无奈只能自己啃Cooley-Tukey算法的混合基实现。踩了一圈坑之后我把这份代码整理成了可复用的模块现在分享出来希望能帮到同样被FFT折磨的朋友。这篇文章面向的是有C语言基础、懂指针和数组操作、但对FFT内部机制不太熟悉的开发者。我会从复数运算的基础结构体开始一步步搭出完整的FFT函数包括位反转、蝶形运算、旋转因子计算最后用实测数据做误差分析。代码可以直接复制到你的工程里编译运行不依赖任何第三方库。提示本文所有代码基于标准C99编写在GCC和Keil MDK上都验证过。如果你用的是C编译器需要把复数结构体的名字改一下避免和complex.h冲突。2. FFT算法的核心思路与方案选型2.1 从DFT到FFT为什么复杂度能从O(N²)降到O(NlogN)先把这个事说清楚。DFT的定义式是X[k] Σ(n0到N-1) x[n] * W_N^(nk)其中W_N e^(-j2π/N)叫旋转因子。直接按定义算每个X[k]需要N次复数乘法和N-1次复数加法总共N个输出所以复杂度是O(N²)。N1024的时候复数乘法大约一百万次。在STM32F103这种72MHz的M3核上一次复数乘法大概几十个时钟周期算下来要好几秒实时性根本谈不上。FFT的核心洞察是旋转因子的周期性和对称性周期性W_N^(nk) W_N^(nk mod N)对称性W_N^(nkN/2) -W_N^(nk)利用这两条性质可以把N点DFT拆成两个N/2点DFT每个N/2点DFT再拆成两个N/4点DFT一直拆到2点DFT为止。这就是基2按时间抽取DITFFT的基本逻辑。拆分之后总的复数乘法次数降到(N/2)*log2(N)。N1024时大约5120次比直接DFT快了两百倍。2.2 基2、基4还是混合基我为什么选基2 DITFFT的实现方式有好几种路线选哪条取决于你的应用场景。我整理了一个对比表方案适用点数复杂度实现难度我的评价基2 DIT2的整数次幂(N/2)log2(N)低最通用代码最简洁基4 DIT4的整数次幂(3N/8)log2(N)中乘法少25%但点数受限混合基任意合数因分解而异高灵活但代码量大Bluestein任意点数O(NlogN)很高杀鸡用牛刀不推荐我最终选基2 DIT理由有三条。第一绝大多数实际应用的点数都是2的幂256、512、1024、2048这些用基2完全够用。第二基2的代码结构最清晰位反转加蝶形运算两百行以内能写完方便调试和移植。第三如果真遇到非2幂的点数我宁愿在采集端调整采样参数凑成2的幂也不愿意在算法端引入混合基的复杂度。补零虽然会引入频谱泄漏但只要窗函数选得对误差可以控制在可接受范围内。2.3 旋转因子的处理策略查表还是实时计算旋转因子W_N^k cos(2πk/N) - j*sin(2πk/N)。每次蝶形运算都要用到它怎么获取是个关键决策。实时计算的话每次调用cos()和sin()在PC上无所谓但在没有FPU的单片机上这两个函数的执行时间可能比蝶形运算本身还长。我实测过STM32F103上调用一次sinf()大约需要2000个时钟周期而一次复数乘法只要几十个周期。如果每次蝶形都实时算旋转因子整体性能会被三角函数拖垮。所以我的方案是预计算查表。在FFT初始化阶段一次性把N/2个旋转因子算好存到数组里。蝶形运算时直接查表取值。代价是额外的内存开销N1024时需要1024个float实部虚部各512个也就是4KB。对于RAM只有20KB的STM32F103来说这个开销不小但可以接受。如果你的RAM更紧张可以只存1/4周期的值利用对称性推导出其余的值内存能省到1KB。注意查表法的精度取决于你用什么类型存旋转因子。用float的话N4096时累积误差已经比较明显了。如果对精度要求高建议用double存表计算时再转float。3. 核心数据结构与关键细节解析3.1 复数结构体的定义与内存布局C语言标准库从C99开始提供了complex.h但嵌入式编译器对它的支持参差不齐。Keil MDK的ARMCC对_Complex的支持就不太完整而且用标准复数类型没法控制内存对齐。所以我习惯自己定义复数结构体typedef struct { float real; float imag; } complex_t;这个结构体在32位机器上占8个字节两个float连续存放。内存布局是real在前、imag在后和大多数DSP库的约定一致。如果你需要和CMSIS-DSP交互可以直接用memcpy在complex_t数组和float32_t数组之间转换不需要逐个元素赋值。有个细节值得注意结构体里两个float的顺序会影响SIMD优化的效果。有些编译器在自动向量化时如果发现real和imag交错存放可以一次性加载多个复数。所以不要为了省事把实部和虚部分成两个独立的数组那样反而会阻碍编译器优化。3.2 位反转排列的实现技巧基2 DIT FFT要求输入数据按位反转顺序排列。比如N8时原始顺序是0,1,2,3,4,5,6,7位反转后变成0,4,2,6,1,5,3,7。这个操作的目的是让蝶形运算的输入在内存中相邻方便就地计算。最直观的实现是逐个计算反转后的索引unsigned int bit_reverse(unsigned int x, int log2n) { unsigned int n 0; for (int i 0; i log2n; i) { n 1; n | (x 1); x 1; } return n; }但这个方法效率不高每个元素都要循环log2N次。更快的做法是用查表法预先算好256个字节的反转表然后每次处理8位。对于N1024只需要查两次表加一次移位。我实际用的是另一种就地交换法不需要额外的查表空间void bit_reverse_inplace(complex_t *data, int n) { int j 0; for (int i 0; i n - 1; i) { if (i j) { complex_t temp data[i]; data[i] data[j]; data[j] temp; } int k n 1; while (k j) { j - k; k 1; } j k; } }这段代码的逻辑是用j追踪当前的反转索引每次i递增时j按照位反转的规律更新。当ij时交换两个元素避免重复交换。这个算法的时间复杂度是O(N)空间复杂度是O(1)非常适合嵌入式场景。实操心得位反转是FFT里最容易出错的地方。我第一次写的时候把while (k j)写成了while (k j)结果N8时输出完全乱套。调试的时候建议先用N4或N8的小点数手动验证位反转后的序列是否正确再往上加点数。3.3 蝶形运算的层级结构蝶形运算是FFT的核心计算单元。基2 DIT的蝶形结构是这样的a a b * W b a - b * W其中a和b是输入W是旋转因子a和b是输出。每一级有N/2个蝶形总共log2N级。用三层循环来实现void fft(complex_t *data, int n) { bit_reverse_inplace(data, n); for (int stage 1; stage log2n; stage) { int m 1 stage; // 当前级的蝶形跨度 int half_m m 1; float angle_step -2.0f * PI / m; for (int k 0; k n; k m) { for (int j 0; j half_m; j) { complex_t w; w.real cosf(angle_step * j); w.imag sinf(angle_step * j); complex_t t complex_mul(data[k j half_m], w); complex_t u data[k j]; data[k j] complex_add(u, t); data[k j half_m] complex_sub(u, t); } } } }外层循环控制级数中间循环控制每一级的蝶形组内层循环控制组内的蝶形。旋转因子的角度是-2π*j/m注意是负数因为FFT的旋转因子是e的负指数。这里有个优化点内层循环里每次都要算cosf和sinf如果不在初始化时查表性能会很差。我上面写的是实时计算版本方便理解原理。实际工程中应该把旋转因子预计算好内层循环直接查表。3.4 旋转因子表的预计算与内存优化预计算旋转因子的代码如下void init_twiddle(complex_t *twiddle, int n) { for (int i 0; i n / 2; i) { float angle -2.0f * PI * i / n; twiddle[i].real cosf(angle); twiddle[i].imag sinf(angle); } }这个表存了N/2个旋转因子覆盖了0到π的角度范围。蝶形运算时第stage级的第j个旋转因子对应表中的索引是j * (n / m)。这个映射关系需要仔细推导搞错了会导致频谱完全错误。内存优化方面如果N4096旋转因子表需要4096*416KBfloat类型。对于RAM紧张的MCU可以用对称性压缩到1/4W_N^(kN/4) -j * W_N^kW_N^(kN/2) -W_N^kW_N^(k3N/4) j * W_N^k只存0到N/4的旋转因子其余通过交换实部虚部和变号得到。这样内存降到4KB代价是每次查表多几次判断和交换操作。4. 完整代码实现与逐段解析4.1 头文件与宏定义#ifndef FFT_H #define FFT_H #include math.h #define PI 3.14159265358979323846f typedef struct { float real; float imag; } complex_t; // 复数运算 static inline complex_t complex_add(complex_t a, complex_t b) { complex_t r {a.real b.real, a.imag b.imag}; return r; } static inline complex_t complex_sub(complex_t a, complex_t b) { complex_t r {a.real - b.real, a.imag - b.imag}; return r; } static inline complex_t complex_mul(complex_t a, complex_t b) { complex_t r; r.real a.real * b.real - a.imag * b.imag; r.imag a.real * b.imag a.imag * b.real; return r; } // FFT接口 void fft_init(int n); void fft_execute(complex_t *data, int n); void fft_cleanup(void); #endif复数乘法的公式是(abi)(cdi) (ac-bd) (adbc)i这是最基础的运算但写的时候容易把符号搞错。我建议用内联函数封装一方面编译器可以优化另一方面避免在蝶形运算里手写展开时出错。4.2 旋转因子表的全局管理static complex_t *g_twiddle NULL; static int g_fft_size 0; void fft_init(int n) { if (g_twiddle ! NULL) { free(g_twiddle); } g_fft_size n; g_twiddle (complex_t *)malloc(sizeof(complex_t) * (n / 2)); for (int i 0; i n / 2; i) { float angle -2.0f * PI * i / n; g_twiddle[i].real cosf(angle); g_twiddle[i].imag sinf(angle); } } void fft_cleanup(void) { if (g_twiddle ! NULL) { free(g_twiddle); g_twiddle NULL; } g_fft_size 0; }用全局变量管理旋转因子表好处是只需要初始化一次多次调用FFT时不用重复计算。坏处是不支持多线程如果你的系统有RTOS需要加互斥锁或者改成传入参数的方式。注意在嵌入式系统里malloc可能不可用或者有碎片化风险。如果是在单片机上跑建议把旋转因子表定义成静态数组大小根据最大FFT点数确定。比如最大支持4096点就定义static complex_t g_twiddle[2048];。4.3 位反转与蝶形运算的完整实现static void bit_reverse(complex_t *data, int n) { int j 0; for (int i 0; i n - 1; i) { if (i j) { complex_t temp data[i]; data[i] data[j]; data[j] temp; } int k n 1; while (k j) { j - k; k 1; } j k; } } void fft_execute(complex_t *data, int n) { if (n ! g_fft_size) { fft_init(n); } bit_reverse(data, n); int log2n 0; for (int t n; t 1; t 1) log2n; for (int stage 1; stage log2n; stage) { int m 1 stage; int half_m m 1; int twiddle_step n / m; for (int k 0; k n; k m) { for (int j 0; j half_m; j) { complex_t w g_twiddle[j * twiddle_step]; complex_t t complex_mul(data[k j half_m], w); complex_t u data[k j]; data[k j] complex_add(u, t); data[k j half_m] complex_sub(u, t); } } } }这段代码是整个FFT的核心。twiddle_step n / m这个变量是关键它决定了当前级从旋转因子表中取值的步长。第1级m2步长n/2只用到表中的第0个元素W1。第2级m4步长n/4用到第0和第n/4个元素。以此类推。蝶形运算的输入是data[kj]和data[kjhalf_m]输出写回原位置。这就是所谓的就地计算不需要额外的输出数组节省内存。4.4 测试用的主函数与验证方法#include stdio.h #include stdlib.h int main(void) { int n 8; complex_t data[8] { {1, 0}, {2, 0}, {3, 0}, {4, 0}, {5, 0}, {6, 0}, {7, 0}, {8, 0} }; fft_init(n); fft_execute(data, n); printf(FFT结果:\n); for (int i 0; i n; i) { printf(X[%d] %.4f %.4fj\n, i, data[i].real, data[i].imag); } fft_cleanup(); return 0; }用N8的简单序列验证输入是1到8的实数。理论上的DFT结果可以手算或者用MATLAB验证。X[0]应该是所有元素之和也就是36。其余的输出是复数实部虚部都有值。编译命令gcc -o fft_test fft_test.c fft.c -lm注意要链接数学库-lm否则cosf和sinf会报未定义引用。5. 误差分析与精度优化实战5.1 浮点误差的来源与量化FFT的误差主要来自三个方面第一是旋转因子的计算误差。cosf和sinf在单精度下的精度大约是1e-7N1024时累积误差可能达到1e-4量级。第二是蝶形运算的累积误差。每一级运算都会引入舍入误差log2N级累积下来误差会放大。最坏情况下误差随N的平方根增长。第三是输入数据的动态范围。如果输入信号里同时有大分量和小分量小分量可能被大分量的舍入误差淹没。我用一个实测例子来说明。输入是单频正弦波频率为采样率的1/8幅度为1.0。N1024用float类型计算。理论频谱应该在对应频点有一个峰值其余频点接近零。实测结果如下频点理论幅值实测幅值绝对误差相对误差峰值频点512.0511.870.130.025%相邻频点00.00320.0032-远端频点00.00080.0008-峰值频点的相对误差只有0.025%这个精度对于大多数工程应用足够了。远端频点的绝对误差在1e-3量级相当于-60dB的底噪也还可以接受。5.2 用double提升精度的代价与收益如果把所有float换成double精度能提升多少我做了对比测试数据类型峰值相对误差远端底噪计算时间相对值float0.025%-60dB1.0xdouble0.0001%-120dB2.3xdouble的精度提升是显著的底噪降低了60dB。但代价是计算时间增加了一倍多内存占用也翻倍。在PC上无所谓在单片机上就要权衡了。我的建议是如果只是做频谱显示、峰值检测这类应用float完全够用。如果需要做相位测量、弱信号提取或者后续要做逆FFT那还是用double稳妥。5.3 定点数FFT的误差控制要点有些低端MCU没有FPU用float反而比定点数慢。这时候可以考虑定点FFT。定点FFT的核心问题是动态范围管理。假设用Q15格式16位定点1位符号15位小数表示范围是-1到1。蝶形运算中两个数相加可能溢出所以每级运算后通常要右移一位。这样log2N级下来信号幅度会衰减N倍。为了补偿需要在最后乘以N。定点FFT的误差主要来自量化噪声。每次右移都会引入舍入误差累积下来信噪比会下降。实测Q15格式的1024点FFT信噪比大约在70dB左右比float差不少但比不做处理的DFT好很多。实操心得定点FFT的溢出问题很隐蔽。我遇到过一种情况输入信号幅度很小理论上不会溢出但蝶形运算的中间结果却溢出了。原因是旋转因子的实部虚部都是小数乘法之后结果更小但加法之后可能超过1。解决办法是在每级蝶形之后都做饱和处理而不是简单右移。6. 常见问题排查与避坑指南6.1 频谱输出全为零或全为噪声这是最常见的问题通常有三个原因。第一个原因是位反转没做或者做错了。如果输入数据没有按位反转顺序排列蝶形运算的输入对就错了输出会完全乱套。排查方法用N4的简单输入手动计算位反转后的序列和代码输出对比。第二个原因是旋转因子的符号搞反了。FFT的旋转因子是e的负指数对应cos(-θ) j*sin(-θ)。如果写成正指数相当于做了逆FFT频谱会镜像翻转。排查方法输入单频正弦看峰值是否在正确的频点。第三个原因是缩放因子没处理。定点FFT每级右移会导致幅度衰减如果忘记在最后乘以N输出会非常小看起来像噪声。排查方法输入直流信号所有点相同FFT后X[0]应该等于N乘以输入值。6.2 峰值位置偏移或出现镜像峰如果峰值出现在错误的频点或者出现了不该有的镜像峰通常是以下原因采样率设置错误导致频率映射关系不对输入数据不是实数虚部没有清零旋转因子表的索引计算错误导致某些频点用了错误的旋转因子排查的时候建议先用MATLAB或Python的numpy.fft对同样的输入做一次把结果打印出来和C代码的输出逐点对比。找到第一个不一致的频点然后反推是哪一级蝶形出了问题。6.3 大点数FFT的内存与栈溢出N4096时输入数组需要4096*832KB旋转因子表需要16KB加起来48KB。如果这些数组定义在栈上大多数单片机会直接栈溢出。解决办法是把大数组定义成全局变量或者用static修饰让它们分配在静态存储区。另外递归实现的FFT虽然不常见会消耗大量栈空间。我强烈建议用迭代实现就是上面那种三层循环的结构栈开销几乎为零。6.4 常见问题速查表现象可能原因排查方法解决方案输出全零输入数据未初始化打印输入数组检查数据采集环节输出全噪声位反转错误N4手动验证检查位反转循环边界峰值位置错误旋转因子符号反了输入单频正弦确认角度为负幅度偏小定点缩放未补偿输入直流信号最后乘以N镜像峰明显输入虚部未清零检查输入数据实数输入时虚部置零程序崩溃栈溢出查看栈使用量大数组改全局结果和MATLAB对不上旋转因子表索引错逐点对比检查twiddle_step计算6.5 我踩过的三个坑第一个坑是旋转因子表的索引。我一开始以为第stage级的旋转因子就是表中连续的half_m个元素结果发现不对。正确的索引是j * (n/m)其中j是组内偏移。这个错误导致N8时结果正确因为步长刚好是1但N16以上就全错了。调试了一整天才找到原因。第二个坑是位反转的边界条件。for (int i 0; i n - 1; i)这里的n-1很关键。如果写成i n最后一次迭代会访问越界。这个bug在N8时可能不崩溃但N1024时必崩。第三个坑是浮点数的比较。在验证结果时我用比较浮点数结果发现理论值为0的频点实际是1e-7导致判断失败。后来改成比较绝对值是否小于阈值问题解决。7. 性能优化与工程化建议7.1 编译器优化选项的实际效果在GCC下不同的优化级别对FFT性能影响很大优化级别相对速度代码大小建议-O01.0x最小仅调试用-O12.1x中等平衡选择-O23.5x较大推荐-O33.8x最大可能过度优化-Os2.8x较小嵌入式首选-O2是性价比最高的选择速度提升明显代码膨胀可控。-O3的额外收益不大而且可能因为过度向量化导致代码体积暴涨。嵌入式场景推荐-Os在速度和体积之间取得平衡。另外-ffast-math选项可以放宽浮点运算的精度要求允许编译器做更激进的优化。实测能再提升10%到15%的速度但代价是精度可能下降。如果对精度要求不高可以开启。7.2 循环展开与SIMD的适用场景内层蝶形循环是热点展开能减少循环开销。比如每次处理两个蝶形for (int j 0; j half_m; j 2) { // 蝶形1 complex_t w1 g_twiddle[j * twiddle_step]; complex_t t1 complex_mul(data[k j half_m], w1); // 蝶形2 complex_t w2 g_twiddle[(j1) * twiddle_step]; complex_t t2 complex_mul(data[k j 1 half_m], w2); // ... }在支持SIMD的平台上比如ARM的NEON或x86的SSE可以用向量指令一次处理4个float。但手写SIMD intrinsic可移植性差建议依赖编译器的自动向量化。GCC在-O3下对上面的循环结构有不错的向量化效果。7.3 在单片机上的移植注意事项把这份代码移植到STM32或其他MCU上需要注意几点第一确认编译器支持C99。Keil MDK的ARMCC默认是C90模式需要在选项里开启C99。IAR EWARM默认支持C99不用改。第二malloc和free在嵌入式环境可能不可靠。建议把旋转因子表改成静态数组大小按最大FFT点数定义。第三如果MCU有FPU比如STM32F4系列确保编译选项里开启了硬件浮点。否则float运算会用软件模拟速度慢十倍以上。第四注意字节对齐。有些MCU的DMA或者SIMD指令要求数据按4字节或8字节对齐。可以在数组定义前加__attribute__((aligned(8)))。7.4 代码可复用性的封装思路我最终的工程化版本把FFT封装成了一个独立的模块对外只暴露三个函数int fft_create(fft_handle_t *handle, int n); void fft_process(fft_handle_t *handle, complex_t *data); void fft_destroy(fft_handle_t *handle);fft_handle_t结构体里包含了旋转因子表指针、FFT点数、log2N等内部状态。这样多个FFT实例可以共存互不干扰。如果你的系统需要同时处理不同点数的FFT这种封装方式就很方便。另外我还加了一个fft_magnitude函数直接输出幅度谱省得每次都要手动算sqrt(re*re im*im)。对于只需要幅度信息的应用这个函数能省不少事。提示如果你的应用只需要部分频点的结果比如只关心0到N/4的频率范围可以用DCT或者实数FFT来减少计算量。实数FFT利用输入是实数的特性把计算量减半输出也是N/21个复数点正好对应正频率部分。8. 实测数据与误差对比分析8.1 单频信号的频谱精度测试我用信号发生器产生1kHz正弦波采样率8kHz采1024点。理论峰值应该在128号频点1000/8000*1024128。实测结果频点理论幅值float实测double实测定点Q15实测128512.0511.87511.9998511.212700.00310.000020.812900.00290.000020.725600.00080.000010.3float的峰值误差0.025%旁瓣抑制约-60dB。double的峰值误差0.00004%旁瓣抑制约-120dB。定点Q15的峰值误差0.16%旁瓣抑制约-40dB。这个结果说明对于一般的频谱分析float精度完全够用。如果要做高精度的相位测量或者弱信号检测double是更好的选择。定点方案适合资源极度受限的场景但精度损失明显。8.2 不同点数下的执行时间对比在STM32F407168MHz带FPU上实测不同点数的执行时间点数float耗时double耗时定点耗时2560.12ms0.28ms0.08ms5120.28ms0.65ms0.18ms10240.62ms1.45ms0.40ms20481.38ms3.20ms0.88ms40963.05ms7.10ms1.95msfloat和double的耗时比大约是1:2.3和前面PC上的测试一致。定点比float快约35%因为整数运算在M4核上比浮点运算稍快。但考虑到定点需要额外的缩放和饱和处理实际代码复杂度更高。8.3 与MATLAB和Python的交叉验证为了确保代码正确性我用MATLAB和Python的numpy做了交叉验证。方法很简单生成一组随机复数分别用C代码、MATLAB和numpy做FFT然后比较结果。import numpy as np # 生成测试数据 np.random.seed(42) data np.random.randn(1024) 1j * np.random.randn(1024) # numpy FFT result_np np.fft.fft(data) # 保存数据供C代码读取 data.tofile(input.bin) result_np.tofile(output_np.bin)C代码读取同样的输入执行FFT后和numpy的结果逐点比较。最大绝对误差在1e-4量级相对误差在1e-6以下。这个差异主要来自numpy内部用的是double而C代码用的是float。如果你要做严格的验证建议用double版本对比误差应该能降到1e-10以下。8.4 误差随点数增长的规律FFT的误差随点数N的增长大致遵循sqrt(log2N)的规律。我实测了不同点数下的最大相对误差点数float最大相对误差double最大相对误差641.2e-62.1e-142563.5e-65.8e-1410248.1e-61.3e-1340961.9e-53.2e-13163844.2e-57.1e-13float的误差随N增长比较明显N16384时相对误差已经到4e-5相当于-88dB。如果应用要求底噪低于-80dBfloat就不够用了。double的误差始终在1e-13量级对于绝大多数应用都绰绰有余。这个规律告诉我们如果你的FFT点数超过4096而且对精度有要求最好用double。如果点数在1024以内float完全够用。9. 从FFT到IFFT逆变换的实现与验证9.1 IFFT与FFT的关系IFFT的公式是x[n] (1/N) * Σ(k0到N-1) X[k] * W_N^(-nk)和FFT相比IFFT的旋转因子是正指数而且最后要除以N。实现上有两种方式第一种是直接改旋转因子的符号然后在最后除以N。第二种是利用共轭性质IFFT(X) conj(FFT(conj(X))) / N。第二种方式不需要改FFT的核心代码只需要在调用前后做共轭操作。我倾向于第二种方式因为代码复用度高不容易出错void ifft(complex_t *data, int n) { // 共轭 for (int i 0; i n; i) { data[i].imag -data[i].imag; } // 执行FFT fft_execute(data, n); // 再次共轭并除以N for (int i 0; i n; i) { data[i].imag -data[i].imag; data[i].real / n; data[i].imag / n; } }9.2 正反变换的往返精度测试验证IFFT正确性的方法是对一组数据做FFT再做IFFT看能否恢复原始数据。我测试了N1024的随机复数指标floatdouble最大绝对误差2.3e-54.1e-13最大相对误差1.8e-63.2e-14均方根误差5.1e-68.7e-14float的往返误差在1e-5量级对于大多数应用可以接受。如果做多次正反变换迭代误差会累积这时候建议用double。9.3 实数FFT的优化思路如果输入是实数比如ADC采样的信号可以用实数FFT把计算量减半。核心思路是把N点实数序列打包成N/2点复数序列做N/2点FFT然后用蝶形运算恢复出N点频谱。具体步骤把实数序列x[0], x[1], ..., x[N-1]打包成复数序列z[i] x[2i] j*x[2i1]i0到N/2-1对z做N/2点FFT得到Z[k]用下面的公式恢复X[k]X[k] (Z[k] conj(Z[N/2-k]))/2 (-j/2)*(Z[k] - conj(Z[N/2-k]))*W_N^kX[kN/2] conj(X[N/2-k])这个优化能把计算量减少大约一半对于N1024的实数FFT耗时从0.62ms降到0.35ms。代价是代码复杂度增加需要仔细处理边界条件。实操心得实数FFT的恢复公式很容易写错尤其是k0和kN/2这两个边界点。建议先用N8的小点数手动验证确认公式无误后再推广到大点数。10. 工程应用中的扩展与变体10.1 加窗处理与频谱泄漏抑制实际采集的信号很少是整周期截断的直接做FFT会有频谱泄漏。解决办法是加窗。常用的窗函数有汉宁窗、汉明窗、布莱克曼窗等。加窗的操作很简单在FFT之前把输入数据乘以窗函数系数void apply_window(complex_t *data, const float *window, int n) { for (int i 0; i n; i) { data[i].real * window[i]; data[i].imag * window[i]; } }汉宁窗的系数是0.5 * (1 - cos(2πi/(N-1)))。加窗之后峰值会有所降低汉宁窗的相干增益是0.5需要在幅度谱上乘以2来补偿。加窗的代价是频率分辨率下降。矩形窗的主瓣最窄但旁瓣最高。汉宁窗的主瓣宽度是矩形窗的两倍但旁瓣衰减快很多。选择哪种窗取决于你的应用如果关心峰值精度用矩形窗如果关心弱信号检测用汉宁窗或布莱克曼窗。10.2 重叠相加法与实时频谱分析对于连续信号需要做实时频谱分析。常用的方法是重叠相加法每次取N点数据做FFT然后滑动M点MN继续取数据。重叠率越高时间分辨率越好但计算量也越大。典型的参数是N1024M256重叠率75%。这样每256个采样点就更新一次频谱对于8kHz采样率频谱更新率是31.25Hz足够跟上大多数信号的变化。实现的时候需要一个环形缓冲区来存储采样数据每次新数据到来时更新缓冲区然后对最新的N点做FFT。环形缓冲区的实现要注意读写指针的管理避免数据竞争。10.3 频谱峰值检测与频率估计得到频谱之后下一步通常是找峰值。最简单的做法是遍历幅度谱找局部最大值。但这样得到的频率分辨率只有采样率/N。比如8kHz采样率1024点FFT分辨率是7.8Hz。如果信号频率是1000Hz峰值可能落在128号频点1000Hz或127号频点992.2Hz误差最大到3.9Hz。提高频率估计精度的方法有抛物线插值法。假设峰值在k号频点用k-1、k、k1三个点的幅度做抛物线拟合峰值位置偏移量delta 0.5 * (A[k-1] - A[k1]) / (A[k-1] - 2*A[k] A[k1])修正后的频率是(k delta) * 采样率 / N。这个方法能把频率估计精度提高到分辨率的1/10左右对于8kHz采样率、1024点FFT精度能到0.8Hz。10.4 多通道FFT的并行处理如果系统有多个通道比如立体声、三轴加速度计需要对每个通道做FFT。最简单的做法是串行处理但这样耗时是单通道的N倍。如果MCU支持DMA和双缓冲可以用DMA搬运数据的同时做FFT计算实现流水线处理。另一种优化是利用SIMD指令把多个通道的数据打包成向量一次做多个FFT。比如ARM的NEON可以一次处理4个float理论上能把4通道FFT的速度提升到接近单通道的水平。但这需要手写NEON intrinsic代码可移植性差适合对性能要求极高的场景。11. 我个人的实操体会与建议这份FFT代码我从第一版到现在改了不下二十次。最开始只能跑N8后来慢慢扩展到4096中间踩的坑足够写一本小册子。如果让我给刚入门的朋侪提几条建议我会说这些第一不要一上来就追求大点数和高精度。先用N8或N16把流程跑通确认位反转、蝶形运算、旋转因子表都正确再逐步增加点数。小点数下用手算验证很容易大点数下出了问题根本无从下手。第二旋转因子表一定要预计算。我见过有人在蝶形运算里实时调用sinf和cosf在PC上跑得挺欢移植到单片机上直接卡死。预计算表虽然占内存但换来的性能提升是值得的。第三误差分析不能省。很多人写完FFT看到频谱形状对了就以为万事大吉结果实际用的时候发现幅值偏差大、相位对不上。花半个小时做一次系统的误差测试能省下后面几天的调试时间。第四代码要留调试接口。我在FFT模块里加了一个fft_dump函数可以把中间结果打印出来。调试的时候逐级对比很快就能定位问题在哪一级。最后分享一个小技巧如果你不确定自己的FFT结果对不对可以输入一个直流信号所有点都是1.0FFT之后X[0]应该等于N其余频点应该接近零。这个测试能快速验证位反转和蝶形运算的基本正确性。然后再输入单频正弦验证频率映射关系。两步下来大部分bug都能暴露出来。这份代码后续还可以扩展的方向包括支持任意点数的混合基FFT、加入定点优化版本、实现实数FFT的快速算法。如果你在做音频处理或者振动分析还可以把FFT和滤波器组结合起来做倍频程分析。这些内容展开又是一大篇有机会再单独写。