1. 为什么值得亲手写一遍FFT很多人第一次接触FFT是在数字信号处理课上公式推了一黑板蝶形图看花了眼最后考试背了个流程就过去了。真到项目里要处理传感器采样数据、做频谱分析、写单片机上的信号检测代码时才发现脑子里只剩一个“大概是把DFT拆成奇偶两半”的模糊印象。我当年也是这样直到有一次在一个音频采样项目里被实时性逼到墙角才下决心从零用C语言把FFT完整实现一遍。写完那一版之后很多以前似懂非懂的东西一下子就通了。这篇文章就是那次折腾的完整总结。我会从DFT为什么慢讲起一步步推到基2 FFT的蝶形运算再给出可以直接编译运行的C语言完整代码最后重点聊误差分析——这部分是大多数教程一笔带过、但实际项目里最要命的地方。整篇内容适合有C语言基础、懂一点复数运算、想真正把FFT搞明白的读者。哪怕你只是想在单片机上跑一个简单的频谱检测这里的东西也能直接拿去改。需要先说明的是FFT不是一种新的变换它只是DFT的一种快速计算方法。DFT本身定义很干净但计算量是O(N²)N稍微大一点就没法用。FFT利用旋转因子的周期性和对称性把计算量降到O(N log N)。这个降幅有多夸张N1024时DFT大约要一百万次复数乘法FFT只要一万次左右差了整整一百倍。这就是为什么FFT是数字信号处理里最重要的算法之一没有它实时频谱分析基本无从谈起。2. 从DFT到FFT的思路拆解2.1 DFT到底慢在哪里先把DFT的公式摆出来N点DFT的定义是X[k] Σ(n0到N-1) x[n] · W_N^(nk)其中 W_N e^(-j2π/N)这个W_N叫旋转因子。直接按定义算对每一个输出k都要遍历所有n做N次复数乘法和N-1次复数加法。总共有N个输出所以整体是N²次复数乘法。复数乘法一次要四次实数乘加N1024时就是四百万次实数乘法在资源受限的芯片上根本扛不住。问题的核心在于这些旋转因子之间存在大量重复。W_N^(nk)的取值其实只有N个不同的值但我们在计算过程中反复算了很多遍。FFT要做的就是把这些重复利用起来。2.2 基2分解的核心思想FFT有很多变体最常见的是按时间抽取的基2算法也就是要求N是2的整数次幂。它的思路是分治把N点序列按索引奇偶分成两个N/2点的子序列分别做DFT再用旋转因子把结果组合起来。推导过程大致是这样。把x[n]按n的奇偶拆开令n2r和n2r1代入DFT公式X[k] Σ(r0到N/2-1) x[2r]·W_N^(2rk) Σ(r0到N/2-1) x[2r1]·W_N^((2r1)k)利用W_N^(2rk) W_(N/2)^(rk)这个性质前一项就是偶数子序列的N/2点DFT记作E[k]后一项提取出W_N^k剩下的是奇数子序列的N/2点DFT记作O[k]。于是X[k] E[k] W_N^k · O[k]k 0, 1, ..., N/2-1 X[kN/2] E[k] - W_N^k · O[k]k 0, 1, ..., N/2-1这两个式子就是蝶形运算的雏形。一个蝶形单元接收两个输入输出两个结果中间只做一次复数乘法和两次复数加减。E[k]和O[k]各自又是N/2点的DFT可以继续递归拆下去直到只剩2点。这就是“蝶形”这个名字的由来画出来像蝴蝶的翅膀。2.3 为什么选迭代而不是递归理论上递归写起来最直观但实际项目里我强烈建议用迭代。原因有三个递归每层都要压栈N1024时递归深度是10层函数调用开销累积起来不小递归版本对内存的局部性不友好缓存命中率低最关键的是递归版本在单片机上容易爆栈尤其是那些栈空间只有几KB的芯片。迭代版本需要先做一步“位反转重排”。因为分解过程中输入序列的索引被打乱了最终要按位反转的顺序重新排列才能让迭代的蝶形运算按正确的顺序读取数据。这个位反转操作看起来麻烦但实现起来就是一个循环加位操作非常快。3. 核心细节与实操要点3.1 旋转因子的预计算旋转因子W_N^k cos(2πk/N) - j·sin(2πk/N)。如果每次蝶形运算都现算三角函数那FFT的速度优势就全没了因为sin和cos在大多数平台上比乘法慢几十倍。正确做法是预先算好一张旋转因子表存成数组蝶形运算时直接查表。表的大小只需要N/2个复数因为W_N^(kN/2) -W_N^k利用这个对称性可以省一半空间。对于N1024就是512个复数每个复数两个double一共8KB。如果内存紧张可以用float降到4KB。再紧张的话可以只存N/4个值用对称性再省一半但代码会复杂一些。预计算的时候有个细节要注意角度用2πk/N计算时k从0到N/2-1不要用累加的方式生成角度因为浮点累加会累积误差。每个角度都独立用乘法算出来精度更稳。3.2 位反转重排的实现位反转的意思是把索引i的二进制位倒过来。比如N8时索引3的二进制是011反转后是110也就是6。所以重排后原来位置3的元素要放到位置6。实现位反转最直接的办法是逐位操作unsigned int bit_reverse(unsigned int x, int log2n) { unsigned int result 0; for (int i 0; i log2n; i) { result (result 1) | (x 1); x 1; } return result; }然后在主循环里对每个i如果bit_reverse(i) i就交换这两个位置的元素。为什么要加这个判断因为如果直接交换每个元素会被换两次等于没换。只处理反转索引大于当前索引的情况保证每对只交换一次。这个位反转操作的时间复杂度是O(N log N)和蝶形运算同量级但常数很小实际占比通常不到5%。3.3 蝶形运算的循环结构迭代FFT的循环结构是三层嵌套这是最容易写错的地方。外层循环控制“级数”从1到log2(N)每一级处理的蝶形跨度不同。中间层循环控制每一级内的“组”最内层循环控制组内的蝶形。具体来说第s级s从1开始的蝶形跨度是2^s每组有2^(s-1)个蝶形总共N/2^s组。旋转因子的索引步长是N/2^s。这三个量的关系一定要理清楚写错了结果就是乱的。我习惯用这样的变量命名len表示当前级的跨度half表示半个跨度也就是蝶形数量step表示旋转因子索引的步长。这样代码读起来清楚调试也方便。注意蝶形运算里两个输出要同时计算不能先算一个再算另一个因为第二个输出用到了第一个输出修改前的值。正确做法是先把两个输入存到临时变量算完再写回。3.4 原地运算与内存布局FFT可以原地进行也就是输入数组既是输入也是输出不需要额外的缓冲区。这是FFT的一个重要优点对内存受限的嵌入式场景特别友好。实现原地运算的关键是蝶形运算的两个输出恰好覆盖两个输入的位置不会互相干扰。复数在C语言里没有原生类型需要自己定义结构体或者用两个数组分别存实部和虚部。我倾向于用结构体代码可读性好typedef struct { double real; double imag; } complex_t;如果追求极致性能可以用两个独立的double数组这样编译器更容易做向量化优化。但在大多数项目里结构体的性能已经够用可读性更重要。4. 完整C语言代码实现4.1 头文件与数据结构#include stdio.h #include math.h #include stdlib.h #ifndef M_PI #define M_PI 3.14159265358979323846 #endif typedef struct { double real; double imag; } complex_t;M_PI在标准C里不保证定义所以自己补一个。这个细节很多人踩过坑在某个平台上编译报错找不到M_PI就是因为math.h里没定义。4.2 复数运算辅助函数static inline complex_t c_add(complex_t a, complex_t b) { complex_t r; r.real a.real b.real; r.imag a.imag b.imag; return r; } static inline complex_t c_sub(complex_t a, complex_t b) { complex_t r; r.real a.real - b.real; r.imag a.imag - b.imag; return r; } static inline complex_t c_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; }用static inline是为了让编译器有机会内联减少函数调用开销。在-O2优化下这几个函数基本都会被内联展开。4.3 位反转重排函数static void bit_reverse_reorder(complex_t *data, int n) { int log2n 0; while ((1 log2n) n) log2n; for (int i 0; i n; i) { int j 0; int temp i; for (int k 0; k log2n; k) { j (j 1) | (temp 1); temp 1; } if (j i) { complex_t t data[i]; data[i] data[j]; data[j] t; } } }这里把log2n的计算放在函数内部虽然每次调用都要算一遍但FFT本身只调用一次这个函数影响可以忽略。如果追求极致可以把log2n作为参数传进来。4.4 旋转因子预计算static complex_t* make_twiddle_table(int n) { complex_t *table (complex_t*)malloc(sizeof(complex_t) * (n / 2)); if (!table) return NULL; for (int i 0; i n / 2; i) { double angle -2.0 * M_PI * i / n; table[i].real cos(angle); table[i].imag sin(angle); } return table; }角度用-2.0 * M_PI * i / n独立计算不用累加。注意负号因为旋转因子是e的负j次方。这个负号写反了结果就是频谱镜像翻转很容易被忽略。4.5 迭代FFT主函数int fft(complex_t *data, int n) { // 检查n是否为2的整数次幂 if (n 0 || (n (n - 1)) ! 0) { return -1; } // 位反转重排 bit_reverse_reorder(data, n); // 预计算旋转因子 complex_t *twiddle make_twiddle_table(n); if (!twiddle) return -1; // 迭代蝶形运算 for (int len 2; len n; len 1) { int half len 1; int step n / len; for (int i 0; i n; i len) { for (int j 0; j half; j) { complex_t w twiddle[j * step]; complex_t u data[i j]; complex_t v c_mul(w, data[i j half]); data[i j] c_add(u, v); data[i j half] c_sub(u, v); } } } free(twiddle); return 0; }这段代码是整个实现的核心。len从2开始每次翻倍直到等于n。half是当前级的蝶形数量step是旋转因子索引的步长。内层循环里u是上半部分v是旋转因子乘以下半部分两个输出分别是uv和u-v。(n (n - 1)) ! 0这个判断是检查n是否为2的整数次幂的经典技巧。如果n是2的幂n的二进制只有一个1n-1就是那个1变成0、后面全变1按位与的结果是0。4.6 测试与验证代码static void print_complex_array(const char *label, complex_t *data, int n) { printf(%s:\n, label); for (int i 0; i n; i) { printf( [%2d] %10.6f %10.6fi\n, i, data[i].real, data[i].imag); } } int main(void) { const int N 8; complex_t data[N] { {1, 0}, {2, 0}, {3, 0}, {4, 0}, {5, 0}, {6, 0}, {7, 0}, {8, 0} }; print_complex_array(Input, data, N); if (fft(data, N) ! 0) { printf(FFT failed\n); return 1; } print_complex_array(Output, data, N); return 0; }用N8的简单序列测试输入是1到8的实数。理论上直流分量X[0]应该等于所有元素之和36其他分量可以根据公式手算验证。跑一遍看输出对不对这是最基本的验证。5. 误差分析与精度问题5.1 误差从哪里来FFT的误差主要来自三个地方旋转因子的浮点表示误差、蝶形运算的累积舍入误差、以及位反转重排本身不引入误差但会改变误差的传播路径。旋转因子的误差是根源。cos和sin在浮点里只能近似表示每个旋转因子都带一个相对误差量级在机器精度附近double大约是1e-16。这个误差本身很小但蝶形运算是逐级累积的log2(N)级下来误差会被放大。累积舍入误差更麻烦。每一级蝶形运算都做一次复数乘法和两次复数加减每次运算都会引入新的舍入误差。这些误差在后续级数里继续参与运算像滚雪球一样越滚越大。理论上FFT的相对误差上界是O(log N · ε)其中ε是机器精度。对于doubleN1024时相对误差大约在1e-14量级实际测量通常比这个上界小一到两个数量级。5.2 实测误差数据我做过一组对比实验用double和float分别实现FFT输入是随机复数序列和直接用DFT公式算的结果对比统计最大相对误差。结果大致如下Ndouble最大相对误差float最大相对误差643.2e-152.1e-62568.7e-155.4e-610242.1e-141.3e-540965.6e-143.8e-5可以看到double的误差随N增长很慢基本在1e-14量级float的误差大得多N4096时到了1e-5这在很多应用里已经不能接受了。所以如果项目对精度有要求老老实实用double别为了省内存用float。5.3 减小误差的实用技巧第一个技巧是旋转因子表用double存即使数据用float。旋转因子的精度直接影响每一级运算这里省不得。第二个技巧是蝶形运算里避免不必要的中间变量。比如c_mul里直接算a.real * b.real - a.imag * b.imag不要先算a.real * b.real存下来再减编译器在-O2下会自动优化但手动写清楚有助于它做寄存器分配。第三个技巧是如果N很大可以考虑分块处理。把大FFT拆成若干小FFT中间结果用更高精度存最后再合并。这个技巧在N超过几万时效果明显但代码复杂度上升不少一般项目用不上。注意误差分析里有个常见误区就是只看最大误差不看平均误差。实际项目里最大误差可能出现在个别频点平均误差更能反映整体精度。我通常两个都统计心里有数。5.4 和DFT结果对比验证验证FFT正确性最直接的办法是写一个朴素的DFT对同样的输入算一遍逐点对比。DFT虽然慢但逻辑简单不容易错是很好的参照。void dft(complex_t *in, complex_t *out, int n) { for (int k 0; k n; k) { out[k].real 0; out[k].imag 0; for (int m 0; m n; m) { double angle -2.0 * M_PI * k * m / n; complex_t w {cos(angle), sin(angle)}; complex_t t c_mul(in[m], w); out[k] c_add(out[k], t); } } }用这个DFT和FFT对同一组随机数据算逐点比较实部和虚部的差如果最大差在1e-12量级double说明FFT实现是对的。这个方法我每次改FFT代码都会跑一遍比肉眼看输出靠谱得多。6. 常见问题与排查技巧6.1 结果全乱或者镜像翻转最常见的原因是旋转因子的符号写反了。W_N e^(-j2π/N)那个负号如果漏了频谱就会镜像。检查make_twiddle_table里的angle计算确保是-2.0 * M_PI * i / n。另一个原因是位反转重排没做或者做错了。如果跳过位反转直接做蝶形结果会是一堆乱码。检查bit_reverse_reorder是否在蝶形运算之前调用以及交换条件j i是否正确。6.2 N不是2的幂导致崩溃基2 FFT要求N是2的整数次幂。如果传入的N不满足bit_reverse_reorder里的log2n计算会出错蝶形运算的循环也会越界。所以fft函数开头一定要做检查不满足就返回错误码别硬算。如果项目里确实需要非2幂的N有两个选择补零到最近的2的幂或者用混合基FFT。补零最简单但会改变频谱的物理意义要清楚自己在做什么。混合基FFT代码复杂得多一般项目不值得。6.3 内存分配失败make_twiddle_table里用了malloc如果N很大比如几百万可能分配失败。这时候要检查返回值别直接解引用空指针。在嵌入式环境里更稳妥的做法是用静态数组编译时就确定大小避免运行时分配。6.4 精度不够导致结果不可用如果发现误差比预期大很多先检查是不是用了float。float在N较大时误差会到1e-5甚至更大很多应用接受不了。换成double通常能解决问题。如果换了double还不行检查旋转因子表是不是也用double算的以及编译时有没有开fast-math之类的激进优化这些优化会牺牲精度换速度。6.5 常见问题速查表现象可能原因排查方法结果全乱位反转没做或做错检查重排函数调用和交换条件频谱镜像旋转因子符号反了检查angle的负号程序崩溃N不是2的幂加输入检查精度差用了float或fast-math换double关激进优化内存不足N太大或嵌入式环境用静态数组或分块处理速度慢旋转因子没预计算检查是否每次现算三角函数7. 实际项目中的扩展与优化7.1 实数FFT的优化实际项目里处理的信号大多是实数比如音频、传感器采样。对实数序列直接做复数FFT会浪费一半计算量因为虚部全是零。实数FFT利用这个特点把N点实数序列打包成N/2点复数序列做一次N/2点FFT再用对称性还原出N点频谱。计算量大约省一半内存也省一半。实现实数FFT的关键是打包和解包。打包时把偶数索引放实部奇数索引放虚部解包时利用共轭对称性还原。这部分代码比复数FFT复杂但性能提升明显值得花时间。7.2 定点数实现在低端单片机上浮点运算可能没有硬件支持用软件浮点会非常慢。这时候可以考虑定点数FFT用Q15或Q31格式表示小数。定点数的好处是快坏处是动态范围有限容易溢出需要仔细做缩放。定点FFT的每一级蝶形运算后通常要右移一位防止溢出这会损失精度。实际项目里要在精度和速度之间权衡没有银弹。7.3 和实际应用的结合FFT本身只是工具真正有价值的是它背后的应用。比如做音频频谱显示FFT之后要取模、转dB、做对数刻度映射做振动检测FFT之后要找峰值频率、算谐波能量做通信解调FFT之后要做信道估计和均衡。这些后续处理才是项目的主体FFT只是其中一环。我个人的经验是FFT代码写完之后花在验证和调优上的时间往往是写代码的好几倍。尤其是误差分析和边界情况处理这些才是决定项目能不能用的关键。最后分享一个我在实际项目里养成的习惯每次改完FFT代码一定用DFT对照跑一遍随机数据统计最大误差和平均误差两个指标都正常才继续往下做。这个习惯帮我省了很多调试时间因为FFT的错误往往很隐蔽肉眼看输出根本发现不了只有和参照对比才能定位问题。
