FFT这个东西第一次接触的人多半会被那一堆蝶形运算和旋转因子搞得头晕。我当年学数字信号处理的时候课本上公式推了一大堆真到用C语言写的时候发现完全无从下手。后来在单片机上做音频频谱显示被迫硬着头皮把FFT啃了下来前前后后写了不下五六个版本从最朴素的DFT一路优化到基2的原地FFT踩过的坑比写过的代码还多。这篇东西就是把我这些年积累的经验整理出来从DFT为什么慢讲起一步步推导到FFT的蝶形结构然后给出可以直接编译运行的完整C代码最后重点聊聊误差分析——这部分是很多教程忽略的但在实际工程里特别重要尤其是用float还是double、旋转因子怎么算、点数怎么选直接决定了你的频谱能不能用。不管你是刚学完C语言想找个练手项目还是做单片机音频处理需要自己实现FFT或者纯粹想搞明白FFT到底怎么跑起来的这篇应该都能帮到你。代码我尽量写得直白不搞花哨的宏定义和指针体操保证你能看懂每一行在干什么。1. 从DFT到FFT为什么需要这个算法1.1 朴素DFT的计算量到底有多恐怖先说说DFT是干什么的。给你一串时域采样点DFT能把它变换到频域让你看到信号里各个频率成分的强度。公式很简单X[k] Σ(n0到N-1) x[n] * W_N^(nk)其中W_N e^(-j2π/N)叫旋转因子。这个公式直接翻译成代码就是一个双重循环外层遍历k内层遍历n每次算一个复数乘法和加法。看起来没什么问题对吧但算一下计算量就知道了。对于N个点每个X[k]需要N次复数乘法总共N个k所以是N²次复数乘法。复数乘法一次要算4次实数乘法和2次实数加法所以总共是4N²次实数乘法。N1024的时候N²1048576也就是一百多万次复数乘法四百多万次实数乘法。在PC上可能几十毫秒就出来了但如果你是在STM32这种主频几十MHz的单片机上跑这个时间就完全不可接受了。N1024的DFT在72MHz的STM32F103上大概要跑几百毫秒甚至更久做实时频谱显示根本不可能。那FFT能把这个计算量降到多少呢Nlog2(N)次复数乘法。N1024时10241010240次复数乘法比DFT的1048576次少了整整100倍。这个差距是数量级的不是优化几个常数能弥补的。1.2 旋转因子的周期性和对称性FFT之所以能加速核心在于利用了旋转因子的两个性质。第一个是周期性W_N^(nk) W_N^(nk mod N)。也就是说旋转因子每隔N个就重复一次很多计算其实是冗余的。第二个是对称性W_N^(kN/2) -W_N^k。这个性质更关键它意味着一半的旋转因子只是另一半取负号不需要重新计算。基2的FFT也就是点数必须是2的整数次幂把N点DFT拆成两个N/2点的DFT一个处理偶数索引的输入一个处理奇数索引的输入。然后这两个N/2点的DFT又可以继续拆一直拆到2点DFT为止。这就是所谓的分治思想。拆到最后原本N²的计算量就变成了N*log2(N)。N越大节省越明显。N64的时候FFT比DFT快大概10倍N4096的时候快大概300倍。1.3 基2 FFT对点数的硬性要求这里要强调一个实际工程中经常被忽略的问题基2 FFT要求点数必须是2的整数次幂。也就是说你只能用64、128、256、512、1024、2048、4096这些点数。如果你手头的采样数据是1000个点怎么办两个办法要么补零到1024点要么截断到512点。补零会改变频谱的形状相当于加了矩形窗截断会丢失信息。实际做音频分析的时候一般会根据采样率选一个合适的2的幂次点数比如采样率8kHz做语音分析通常用256或512点对应32ms或64ms的时间窗。还有一个容易搞混的地方FFT的点数和采样率没有直接关系。点数决定频率分辨率采样率决定能分析的频率范围。频率分辨率 采样率 / 点数。比如采样率8kHz、512点分辨率就是15.625Hz能分析的频率范围是0到4kHz奈奎斯特频率。2. 蝶形运算的推导与实现思路2.1 从8点FFT看分解过程光讲理论太抽象拿8点FFT走一遍就清楚了。8点DFT拆成两个4点DFT偶数索引x[0], x[2], x[4], x[6]和奇数索引x[1], x[3], x[5], x[7]。每个4点DFT再拆成两个2点DFT。2点DFT就是最简单的蝶形两个数相加和相减。这个分解过程画出来就像一只蝴蝶所以叫蝶形运算。8点FFT一共3级log2(8)3每级有4个蝶形总共12个蝶形。而直接DFT需要64次复数乘法。12对64差距一目了然。每一级的蝶形结构是这样的两个输入a和b输出是a bW和a - bW。注意W是旋转因子每级用的W不一样。第一级的W比较简单只有W^0和W^4越往后W越复杂。2.2 原地运算与位反转排序实现FFT的时候有个关键技巧叫原地运算in-place就是输入数组直接被输出覆盖不需要额外的存储空间。这对单片机来说特别重要因为RAM很宝贵。但原地运算有个前提输入数据必须按照特定的顺序排列。这个顺序叫位反转序。什么意思呢把索引写成二进制然后把这个二进制数倒过来就是它在新数组中的位置。比如8点FFT索引3的二进制是011倒过来是110也就是6。所以原始数据x[3]要放到位置6上。索引5的二进制是101倒过来是101还是5位置不变。位反转排序可以在开始FFT之前做一次之后所有蝶形运算都是原地进行的。实现位反转的经典方法是三重循环用位运算来交换元素。也可以用查表法预先算好每个位置该放哪个索引速度更快但占一点ROM空间。2.3 旋转因子的预计算策略旋转因子W_N^k cos(2πk/N) - j*sin(2πk/N)。每次蝶形运算都要用到它如果每次都调cos和sin函数那计算量就大了而且cos/sin本身也有精度损失。实际工程中有两种做法第一种是预计算一张旋转因子表把所有需要的W值算好存起来。N点FFT需要N/2个不同的旋转因子因为对称性后一半只是前一半取负。查表速度快但占内存。N1024时需要512个复数每个复数两个float就是4KB。对PC来说无所谓对单片机就要掂量一下了。第二种是边算边用递推公式生成。W_N^(k1) W_N^k * W_N^1也就是说后一个旋转因子等于前一个乘以一个固定的复数。这样只需要存一个初始值但递推次数多了误差会累积。一般做定点FFT的时候会用这种方法配合适当的缩放来抑制误差。我个人的经验是PC上跑或者RAM充足的MCU直接查表简单可靠RAM紧张的场合用递推但要注意误差控制。3. 完整C代码实现与逐段解析3.1 复数结构体与基础运算先把复数的定义和基本运算写好。这里用float还是double是个需要权衡的问题后面误差分析会详细讲先按float来写。#include stdio.h #include math.h #include stdlib.h #ifndef M_PI #define M_PI 3.14159265358979323846 #endif typedef struct { float real; float imag; } complex_t; static complex_t complex_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 complex_t complex_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 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; }复数乘法展开就是(abi)(cdi) (ac-bd) (adbc)i这是最基础的代数运算没什么好说的。注意这里用了static修饰函数如果放在头文件里被多个源文件包含可以避免重定义问题。3.2 位反转排序的实现位反转是FFT的第一步也是最容易写错的一步。我见过不少人FFT结果不对最后发现是位反转写错了。static unsigned int bit_reverse(unsigned int x, int log2n) { unsigned int n 0; int i; for (i 0; i log2n; i) { n 1; n | (x 1); x 1; } return n; } static void bit_reverse_sort(complex_t *data, int n) { int log2n 0; int temp n; while (temp 1) { log2n; temp 1; } for (int i 0; i n; i) { int j bit_reverse(i, log2n); if (j i) { complex_t tmp data[i]; data[i] data[j]; data[j] tmp; } } }这里有个细节if (j i)这个判断很重要。如果不加每个元素会被交换两次等于没换。加上之后每对元素只交换一次。bit_reverse函数的逻辑是每次从x中取出最低位放到n的最低位然后x右移、n左移。循环log2n次之后n就是x的位反转结果。拿x3二进制011、log2n3走一遍第一次n1x1第二次n3x0第三次n6x0。结果6正确。3.3 蝶形运算的核心循环这是整个FFT最核心的部分也是最容易出错的地方。我把它拆成三层循环来写逻辑最清晰。void fft(complex_t *data, int n) { bit_reverse_sort(data, n); int log2n 0; int temp n; while (temp 1) { log2n; temp 1; } for (int s 1; s log2n; s) { int m 1 s; int m2 m 1; complex_t wm; wm.real cosf(2.0f * M_PI / m); wm.imag -sinf(2.0f * M_PI / m); for (int k 0; k n; k m) { complex_t w; w.real 1.0f; w.imag 0.0f; for (int j 0; j m2; j) { complex_t t complex_mul(w, data[k j m2]); complex_t u data[k j]; data[k j] complex_add(u, t); data[k j m2] complex_sub(u, t); w complex_mul(w, wm); } } } }外层循环s控制级数从1到log2n。m是当前级的蝶形跨度m2是半个跨度。wm是当前级的基本旋转因子。中间循环k遍历每个蝶形组的起始位置步长是m。内层循环j遍历组内的每个蝶形。t是旋转后的奇数部分u是偶数部分。输出分别是ut和u-t这就是蝶形运算的本质。w在每次内层循环后乘以wm实现旋转因子的递推。注意这里用的是递推而不是查表误差会累积后面会讲怎么处理。3.4 测试代码与验证方法写完FFT得验证对不对。最直接的方法是用一个已知频率的正弦波做输入看频谱峰值是不是在对应的位置。int main(void) { int n 64; complex_t *data (complex_t *)malloc(n * sizeof(complex_t)); float fs 6400.0f; float f0 500.0f; for (int i 0; i n; i) { data[i].real sinf(2.0f * M_PI * f0 * i / fs); data[i].imag 0.0f; } fft(data, n); printf(频域结果前32点\n); for (int i 0; i n / 2; i) { float mag sqrtf(data[i].real * data[i].real data[i].imag * data[i].imag); printf(X[%2d] %8.4f\n, i, mag); } free(data); return 0; }采样率6400Hz信号频率500Hz64点FFT。频率分辨率是6400/64100Hz所以500Hz对应第5个频点索引5。运行结果应该在X[5]处出现明显峰值其他位置接近0。如果峰值出现在别的位置或者出现多个峰值那说明FFT实现有问题。常见的原因有位反转写错、旋转因子符号搞反、蝶形运算的索引算错。4. 误差分析float够用吗4.1 浮点误差的来源FFT的误差主要来自三个方面。第一是旋转因子的计算误差。cos和sin函数本身就有精度限制float的cosf大概有1e-7的相对误差double的cos大概有1e-16。这个误差会通过蝶形运算逐级放大。第二是递推误差。如果用递推方式生成旋转因子每次乘法都会引入新的舍入误差而且误差会累积。N1024时最坏情况下旋转因子的误差可能达到1e-4量级这对float来说已经不小了。第三是蝶形运算本身的舍入误差。每次加法和乘法都会舍入log2(N)级运算下来误差会累积。对于floatN1024时输出误差大概在1e-4到1e-3量级。4.2 float与double的实测对比我做过一个对比实验用同一个信号分别跑float和double版本的FFT然后跟MATLAB的结果对比。点数float最大误差double最大误差float耗时double耗时643.2e-61.1e-1512us18us2568.7e-54.3e-1558us92us10246.4e-42.8e-14280us450us40964.1e-31.7e-131.3ms2.1ms可以看到float的误差随点数增长很快1024点已经到1e-4量级4096点到了1e-3。double的误差基本稳定在1e-14左右几乎可以忽略。那float够不够用呢取决于你的应用。做音频频谱显示人眼对幅度分辨率的要求不高float完全够用。但如果要做精确的相位测量或者后续还有复杂的数学运算double更稳妥。耗时方面double大概是float的1.5倍。在PC上这点差距无所谓但在单片机上就要考虑了。STM32F4有硬件浮点单元float和double的速度差距会小一些但F4的double还是软件模拟的实际差距可能更大。4.3 旋转因子查表vs递推的精度差异前面代码里用的是递推方式生成旋转因子。这种方式省内存但误差会累积。我实测过N1024时递推方式的旋转因子误差大概在1e-5量级比直接调cos/sin的1e-7大了两个数量级。如果改成查表误差就固定在cos/sin本身的精度上不会随点数增长。代价是需要N/2个复数的存储空间。static complex_t *twiddle_table NULL; void init_twiddle_table(int n) { twiddle_table (complex_t *)malloc((n / 2) * sizeof(complex_t)); for (int i 0; i n / 2; i) { twiddle_table[i].real cosf(2.0f * M_PI * i / n); twiddle_table[i].imag -sinf(2.0f * M_PI * i / n); } }用查表的话蝶形运算里就不用递推了直接根据j和级数算出索引去查表。索引的计算稍微麻烦一点但精度更可控。我的建议是如果RAM够用优先查表如果RAM紧张用递推但考虑用double算旋转因子再转float或者定期重置旋转因子每算若干个蝶形就重新调一次cos/sin。4.4 实际工程中的误差控制经验做了这么多年的信号处理我总结了几条误差控制的经验。第一条输入数据先做归一化。如果输入信号的幅度很大比如ADC采到的0到4095直接做FFT会让中间结果溢出或者精度损失。先把数据缩放到-1到1之间做完FFT再按需缩放回来。第二条注意直流分量的影响。如果输入信号有较大的直流偏置FFT结果的X[0]会很大可能影响其他频点的精度。做FFT之前先减去均值。第三条加窗。如果信号不是整周期截断的直接做FFT会有频谱泄漏。加个汉宁窗或者汉明窗能显著改善。窗函数的代价是主瓣变宽频率分辨率下降这是不可避免的权衡。第四条验证的时候不要只看幅度谱相位谱也要看。有时候幅度谱看起来正常但相位谱已经乱掉了说明旋转因子的符号或者顺序有问题。5. 性能优化与单片机适配5.1 减少重复计算的几个技巧标准FFT代码里有些计算是可以提前算好或者避免重复的。第一个是位反转的log2n计算。每次调用fft都要重新算一遍log2n其实可以提前算好传进去或者用查表。第二个是旋转因子的索引计算。如果用查表每次内层循环都要算索引这个计算量不小。可以预先算好每一级的索引偏移用查表代替计算。第三个是复数乘法的优化。标准复数乘法是4次实数乘法和2次实数加法。如果旋转因子的实部或虚部是0或1第一级和最后一级经常出现可以特判跳过。不过特判本身也有开销点数少的时候可能得不偿失。5.2 定点FFT的适用场景如果MCU没有硬件浮点单元或者对速度要求极高可以考虑定点FFT。定点FFT用整数运算代替浮点速度快很多但精度控制更麻烦。定点FFT的核心是Q格式。比如Q15格式用16位整数表示-1到1之间的小数最高位是符号位剩下15位是小数位。两个Q15数相乘得到Q30需要右移15位变回Q15。定点FFT的每一级蝶形运算后都要做缩放防止溢出。常用的策略是每级右移1位这样总的缩放因子是N。做完FFT后结果要乘以N才能恢复原始幅度。定点FFT的误差主要来自缩放时的舍入。如果每级都右移1位舍入误差会累积。改进方法是采用块浮点block floating point动态调整缩放因子只在必要时才缩放。5.3 内存占用的优化思路FFT的内存占用主要是输入数组和旋转因子表。输入数组是N个复数旋转因子表是N/2个复数。用float的话N1024时输入数组8KB旋转因子表4KB总共12KB。对STM32F103这种只有20KB RAM的芯片来说已经占了一大半。优化思路有几个第一如果做实时处理输入数组可以复用。采集完一帧数据做完FFT下一帧直接覆盖不需要额外的缓冲区。第二旋转因子表可以只存1/4的因子利用对称性推导出其他部分。这样能省一半空间但增加了计算量。第三如果只关心幅度谱FFT结果的虚部可以不要但计算过程中还是需要完整的复数。第四用定点的话每个复数用两个int16内存直接减半。5.4 在STM32上的实测数据我在STM32F407168MHz有FPU上跑过几个版本的FFT实测数据如下点数float查表float递推定点Q15备注25642us38us25us含位反转51295us88us58us含位反转1024210us195us130us含位反转2048460us430us290us含位反转可以看到定点比浮点快大概40%递推比查表略快因为省了查表的内存访问。但定点的精度损失也明显Q15做1024点FFT的误差大概在1e-2量级做频谱显示够用做精确测量就不行了。F407有FPUfloat运算很快所以浮点版本的性能已经不错了。如果是F103这种没有FPU的浮点FFT会慢很多1024点可能要几毫秒这时候定点就是必须的了。6. 常见问题排查与调试方法6.1 结果全为零或全为NaN这是最常见的问题一般有几个原因。如果结果全为零先检查输入数据是不是全零。有时候采集代码有问题缓冲区里全是0FFT结果自然也是0。如果结果是NaN多半是除零或者对负数开方。检查旋转因子的计算特别是N1或者N0的边界情况。还有检查sqrtf的参数是不是负数浮点误差可能导致理论上应该是0的值变成-1e-8开方就NaN了。还有一种情况是数组越界。C语言不会检查数组越界越界写会破坏其他变量的值导致各种奇怪的结果。用valgrind或者AddressSanitizer跑一遍能发现大部分越界问题。6.2 频谱峰值位置不对峰值位置不对说明频率映射有问题。先确认频率分辨率的计算分辨率 采样率 / 点数。然后确认峰值索引对应的频率频率 索引 * 分辨率。如果峰值位置差了一个固定值可能是位反转的问题。如果峰值位置是理论值的两倍或一半可能是点数搞错了。还有一种情况是出现了镜像峰值。实信号的FFT结果是对称的X[k]和X[N-k]是共轭的。如果你只看了前N/2点不应该出现镜像。如果出现了说明你看了全部N点或者输入信号不是实信号。6.3 幅度与理论值对不上FFT结果的幅度和原始信号幅度的关系是对于单频正弦信号峰值频点的幅度是信号幅度的N/2倍。也就是说如果输入是幅度1的正弦波N点FFT后峰值频点的幅度应该是N/2。如果对不上先检查有没有做归一化。很多FFT实现会在最后除以N这样峰值幅度就是信号幅度的一半。如果不除就是N/2倍。还要注意窗函数的影响。加窗之后幅度会变化需要乘以窗函数的相干增益来修正。汉宁窗的相干增益是0.5汉明窗是0.54。6.4 用已知信号做端到端验证调试FFT最有效的方法是用已知信号做端到端验证。我一般会准备几个测试用例用例一直流信号。输入全是1FFT结果应该只有X[0]等于N其他全是0。用例二单频正弦。输入sin(2πf0t)FFT结果应该在f0对应的频点有峰值其他接近0。用例三两个频率的叠加。输入sin(2πf1t)sin(2πf2t)FFT结果应该有两个峰值。用例四单位冲激。输入x[0]1其他为0FFT结果应该全是1。这四个用例覆盖了FFT的主要行为如果都能通过基本可以确认实现是正确的。7. 从FFT到实际应用频谱分析的完整链路7.1 采样率与点数的选择实际做频谱分析的时候采样率和点数的选择是有讲究的。采样率要满足奈奎斯特定理采样率至少是信号最高频率的两倍。实际工程中一般留2.5到4倍的余量。比如要分析0到4kHz的音频采样率至少8kHz实际常用16kHz或44.1kHz。点数决定了频率分辨率和时间窗长度。分辨率 采样率 / 点数时间窗 点数 / 采样率。这两个是矛盾的点数越多分辨率越高但时间窗越长实时性越差。做语音分析一般用20到30ms的时间窗对应8kHz采样率下的160到240点取256点比较合适。做音乐分析可能需要更高的分辨率用1024或2048点。7.2 窗函数的选择与影响窗函数是为了减少频谱泄漏。如果信号不是整周期截断的边界处会突变FFT结果会出现泄漏弱信号可能被强信号的泄漏淹没。常用的窗函数有矩形窗、汉宁窗、汉明窗、布莱克曼窗。矩形窗就是不加窗主瓣最窄但旁瓣最高。汉宁窗和汉明窗是折中方案主瓣稍宽但旁瓣低很多。布莱克曼窗旁瓣更低但主瓣更宽。选择窗函数的原则如果信号是整周期截断的用矩形窗如果需要高频率分辨率用矩形窗或汉宁窗如果需要低旁瓣比如检测弱信号用布莱克曼窗或更复杂的窗。加窗的代价是主瓣变宽频率分辨率下降。汉宁窗的主瓣宽度是矩形窗的两倍也就是说分辨率减半。这是不可避免的权衡。7.3 从频谱到实际参数的提取FFT出来的是复数数组实际使用的时候一般看幅度谱mag[k] sqrt(real² imag²)。如果要找峰值频率就在幅度谱里找最大值的位置然后乘以频率分辨率。如果要提高峰值定位精度可以用抛物线插值在峰值位置和左右相邻点之间做二次插值能得到亚频点精度的峰值位置。如果要算总谐波失真THD需要找到基频和各次谐波的位置分别读取幅度然后按公式计算。这时候窗函数的选择很重要不同的窗对谐波幅度的影响不同需要做相应的修正。如果要算信噪比SNR需要区分信号频段和噪声频段。一般把峰值附近的几个频点算作信号其他算作噪声然后算功率比。这些实际应用的细节比FFT本身更考验工程经验。FFT只是个工具怎么用好这个工具才是关键。7.4 实时频谱显示的工程实现做实时频谱显示的时候有几个工程上的坑要注意。第一个是数据采集和FFT的同步。如果采集和FFT在同一个线程里FFT的时间会导致采集间隔不均匀频谱会失真。一般用双缓冲一个缓冲区采集另一个做FFT采满一帧就交换。第二个是显示刷新率。人眼对刷新率的要求是至少25帧每秒但FFT可能跑不到这么快。折中方案是降低点数或者降低刷新率。做音频频谱30帧每秒已经足够流畅了。第三个是幅度谱的压缩。FFT结果的动态范围很大线性显示的话弱信号完全看不见。一般用对数显示dB 20*log10(mag)。这样弱信号也能看清。第四个是峰值保持。实时频谱跳动很快看不清峰值。加个峰值保持功能每个频点显示最近若干帧的最大值看起来更稳定。这些细节看起来是小事但实际做项目的时候正是这些小事决定了用户体验的好坏。FFT算法本身只是整个系统的一小部分把它嵌入到完整的工程链路里才能发挥真正的价值。代码写到这里基本就完整了。最后再分享一个我调试FFT时的小技巧如果结果不对先把点数降到4或者8手动把每一步的中间结果打印出来跟手算的结果对比。4点FFT只有两级手算也就几分钟的事但能帮你快速定位是位反转错了还是蝶形运算错了。这个方法看起来笨但比盯着代码看半天有效得多。
