C语言手写基2 FFT：从原理到嵌入式频谱分析实战

简介这是一份以C语言实现的快速傅立叶变换FFT源代码包面向数字信号处理学习者、嵌入式开发者以及需要频域分析的技术人员可帮助理解并落地DFT/FFT算法。压缩包共2个文件整体仅2KB其中cpp为算法主源码风格简洁明了txt为辅佐说明文件记录来源或背景资料。代码采用蝶形运算与分治策略将离散傅立叶变换DFT的计算复杂度从O(N^2)降至O(N log N)包含位反转、序列分解、前向与反向变换等关键步骤。离散傅立叶变换的典型公式X[k]Σ x[n]e^{-i2πkn/N}揭示了时域信号到频谱的本质转换而该源码正好可对照公式逐行理解也可直接提取用于图像处理、音频分析、通信调制等实际实验。目前已有251人学习浏览适合课程设计、算法研习与工程参考是入门快速傅立叶变换、提升C语言信号处理能力的紧凑实用资源。1. 从 DFT 到 FFT一段标准的数学捷径FFT快速傅里叶变换不是一种新的变换而是离散傅里叶变换DFT的高效计算方法。DFT 的复杂度是 O(N²)当 N1024 时大约需要一百万次复数乘法而 FFT 利用旋转因子的周期性和对称性把复杂度压到 O(N log N)同样的 1024 点变换只需要约一万次运算。这个差距在实时系统里是本质性的前者做一帧频谱分析可能耗时几十毫秒后者可以做到微秒级直接决定了嵌入式设备上能不能跑得起实时频谱显示。用 C 语言实现 FFT核心价值在于可控性——没有运行时依赖没有内存托管的意外停顿精确知道每一条数据通路和每一个蝴蝶运算的开销。本文覆盖从旋转因子计算到位反转寻址再到定点化与频谱校正的完整落地路径适合需要在 MCU、DSP 或 Linux 用户态程序里自己实现 FFT而非调用现成库的开发者。2. FFT 的数学结构与 C 语言实现前的三个关键决策2.1 按时间抽取DIT还是按频率抽取DIF基 2 FFT 有两个经典变体DITDecimation In Time和 DIFDecimation In Frequency。DIT 把输入序列按下标奇偶拆分先在时域做位反转再从前向后做蝴蝶运算输出自然顺序DIF 则相反输入自然顺序输出位反转蝴蝶的方向是从后往前。嵌入式工程里绝大多数场景选 DIT原因是输出顺序和输入顺序的约定更自然后续接频谱分析、窗函数时心智负担小。DIF 的优势是可以在原位做某些定点化优化但代码可读性差调试成本高除非是极端资源受限否则不推荐新手起步用 DIF。// 典型DIT基2 FFT的屎壳郎循环骨架完整实现见3.3节 for (len 2; len n; len 1) { step len 1; for (i 0; i n; i len) { for (j 0; j step; j) { // 蝴蝶运算对(ij)和(ijstep)两个点做复数旋转 } } }这段骨架不是可运行代码只是为了展示 DIT 的三层循环结构外层控制蝶形跨度len中层定位每一组蝶形的起始位置内层遍历组内的step个旋转因子。理解这个结构之后再往里填具体的复数运算是很自然的事。结构确定后真正的决策点有三个旋转因子索引怎么映射、数据是复数数组还是实虚分离、位反转用什么策略。这三个决策的质量决定移植性和数值精度也决定后续排错容易不容易。2.2 旋转因子表预计算还是实时计算旋转因子 W 定义为 e^(-j2πk/N)k0,1,…,N/2-1。常见做法有两种一是每次蝴蝶运算时调用cosf/sinf实时计算二是程序启动时把 N/2 个复常数算好存表。实时计算节省 RAM 但消耗 CPU且连续调用三角函数在低端 MCU 上可能产生不可忽略的延迟预计算则以 RAM 换速度1024 点 FFT 需要 512 个复数每个复数两个 float共 4KB对绝大多数嵌入式平台毫无压力。预计算表的生成有一个工程技巧不要为每一级单独算表而是用统一的 N/2 点主表每一级的旋转因子通过索引映射取用。索引规律是当蝶形跨度为2*steplen时第 m 组第 j 个旋转因子的数组下标为j * (N / len)。这个映射关系是 FFT 实现中最容易写错的地方写成公式放在代码里作为注释比任何时候都重要。// 预计算旋转因子表浮点版本n为FFT点数twiddle为输出表 void twiddle_init(int n, float complex *twiddle) { for (int k 0; k n / 2; k) { float angle -2.0f * PI * k / n; twiddle[k] cosf(angle) sinf(angle) * I; // C99复数语法 } }参数说明n必须是 2 的幂否则索引映射完全不成立twiddle数组长度为n/2只存半个周期是因为另外一半互为相反数可以通过取负得到不必浪费 RAM。角度为负是因为 FFT 正变换的定义里包含负指数若想做逆变换则取正号。2.3 数据排布复数交错数组与实虚分离数组C 语言没有原生复数类型的历史遗留C99 有complex但嵌入式编译器不一定完整支持工程上有两种排布方式。复数交错数组是把实部和虚部交错存放float data[2*n]第 i 个复数的实部在data[2*i]虚部在data[2*i1]。实虚分离则是两个独立数组float re[n]、float im[n]。交错数组的优点是数组下标运算统一传参时只需一个指针缺点是 SIMD 优化时数据混排不友好。实虚分离的优点是方便和 ADC 采集的纯实数序列对接——直接填实部数组虚部清零——同时某些 MCU 的 DSP 库天然采用分离格式。个人经验是除非你确定要用厂商提供的 DSP 库否则一律用交错数组。原因很简单调试方便打印一个复数只需要一条语句而分离数组要两句且位反转交换时两个数组要同时交换漏一个就是隐蔽的 bug。3. 手写基 2 FFT 的完整 C 代码与逐段拆解3.1 位反转寻址用循环移位代替查表DIT FFT 要求输入序列按下标二进制位反转后的顺序排列。N8 时原始下标 1001需要和第 4100个数据交换。教科书写法是逐位判断工程上有两种更优策略预先生成位反转表或运行时用循环右移累加生成下一个位反转下标。// 原位位反转n必须为2的幂n_log2为log2(n) void bit_reverse(float complex *buf, int n, int n_log2) { int rev 0; for (int i 1; i n; i) { // 从上一个反转值递推当前反转值 int bit n 1; while (rev bit) { rev ^ bit; bit 1; } rev ^ bit; if (i rev) { // 只交换一次避免重复换回原样 float complex tmp buf[i]; buf[i] buf[rev]; buf[rev] tmp; } } }逻辑说明rev始终是当前下标i-1的位反转值通过从最高位向最低位寻找第一个 0 位并置 1、之前的 1 全部清 0得到下一个位反转值。这个技巧避免查表换来的是每次迭代的几步位运算在 N 不大时比查表更省内存。if (i rev)条件保证每对元素只交换一次。手动验证 N8 的序列i1 时 rev4交换 buf[1] 和 buf[4]i3 时 rev6交换 buf[3] 和 buf[6]。完成后序列顺序为 0,4,2,6,1,5,3,7正好是原始下标的位反转。3.2 蝴蝶运算的计算流图与复数旋转实现基 2 蝴蝶的核心公式是一个复数加减乘组合设输入为 x1 和 x2旋转因子为 W则输出为 x1x1Wx2x2x1-Wx2。展开成实部虚部运算一共需要 4 次实数乘法和 6 次实数加减法。// 就地复数乘法*dst *src * ww为复常数 void cmul(float complex *dst, const float complex *src, float complex w) { float re crealf(*src) * crealf(w) - cimagf(*src) * cimagf(w); float im crealf(*src) * cimagf(w) cimagf(*src) * crealf(w); *dst re im * I; }参数说明crealf和cimagf是 C99 标准宏分别取复数的实部和虚部I是虚数单位。如果编译器不支持复数可以把float complex替换为struct { float re; float im; }乘法逻辑完全一样。这里的重点在于结果写回前必须用临时变量保存输入值不能直接覆盖否则后续运算读到的就是被污染的数据。蝴蝶运算的完整实现可以借助 C99 复数运算符直接写// 原位基2蝶形运算buf为数据指针j为组内偏移step为半跨度 buf[i j step] buf[i j] - twiddle[j * (N / len)] * buf[i j step]; buf[i j] buf[i j] twiddle[j * (N / len)] * buf[i j step];注意这里必须先把右边的buf[ijstep]读到临时变量否则第一行执行后第二行读到的已经是旋转后的值。C 语言的求值顺序和赋值副作用是 FFT 代码最经典的坑没有之一。3.3 完整可编译的 C 语言 FFT 主流程把前面两个模块组合起来再加上外层三层循环就是一个完整的基 2 FFT。以下代码在 GCC 和 Clang 下可直接编译运行N16 的测试输入可以手动验证。#include complex.h #include math.h #include stdio.h #define PI 3.14159265358979323846 void fft(float complex *buf, int n) { int n_log2 0; while ((1 n_log2) n) n_log2; // 位反转 预计算旋转因子 bit_reverse(buf, n, n_log2); float complex twiddle[n / 2]; for (int k 0; k n / 2; k) { twiddle[k] cexp(-2.0f * PI * I * k / n); } // 蝶形主循环 for (int len 2; len n; len 1) { int step len 1; for (int i 0; i n; i len) { for (int j 0; j step; j) { float complex w twiddle[j * (n / len)]; float complex t w * buf[i j step]; buf[i j step] buf[i j] - t; buf[i j] buf[i j] t; } } } }逻辑说明外层len从 2 开始倍增直到 nstep是每组蝶形的子集大小中层i以len为步长遍历所有组内层j遍历组内每个蝶形。旋转因子索引j * (n / len)的规律当 len2 时只有 W^0len4 时有 W^0 和 W^2len8 时有 W^0、W^2、W^4、W^6正好对应主表间隔n/len取样。这段代码没有做任何优化但结构清晰可以作为基准实现后续再做定点化或 SIMD 优化时以它的输出为对照。验证方法输入buf [1,0, 2,0, 3,0, 4,0]N4FFT 结果应为[100i, -22i, -20i, -2-2i]。如果结果不对优先检查位反转交换是否成对、旋转因子符号是否一致。4. 实数输入优化与浮点精度控制4.1 实数 FFT 的打包技巧一次变换算出两路实数序列频谱大多数实际输入是实数——ADC 采样的电压、音频 PCM 数据都是实数序列。直接喂给复数 FFT 是把虚部填 0一半的运算是白算的。标准优化手法是实序列打包把两个不同实数序列 x 和 y 分别放入复数数组的实部和虚部做一次复数 FFT再从结果中拆出两个实数序列各自的频谱。拆分的数学依据是共轭对称性。设 z FFT(x iy)则 x 的频谱 X[k] (Z[k] conj(Z[N-k])) / 2y 的频谱 Y[k] (Z[k] - conj(Z[N-k])) / (2i)。这个技巧把 FFT 调用次数减半代价是拆包时的几十次复数运算在 N 较大时收益明显。// 从复数FFT结果Z中拆分出实序列x和y的频谱 float complex X[N], Y[N]; for (int k 0; k N; k) { int k2 (N - k) (N - 1); // 模N取反索引防止k0时越界 X[k] (Z[k] conj(Z[k2])) / 2.0f; Y[k] (Z[k] - conj(Z[k2])) / 2.0f * (-I); }参数说明k2的计算用按位与代替取模必须在 N 为 2 的幂时才成立k0 时 k20正好对应直流分量的实部。拆分后 X[0] 是纯实数包含 x 的直流分量Y[0] 同理。这套技巧要求 FFT 点数为 2 的幂且两个输入序列长度相同否则索引对应关系完全混乱。4.2 浮点误差来源与逐级误差上限浮点 FFT 的误差来自两个地方旋转因子表的近似cosf 本身有约 1e-7 的相对误差和每级蝶形的舍入误差。逐级分析时每次蝶形产生大约 1 个ulp 的误差单精度约为 1.2e-7共 log2(N) 级所以总的相对误差上界约为 log2(N) 个 ulp。N1024 时理论最坏误差约为 1.2e-4实际统计上远小于这个数因为误差符号随机大部分相互抵消。控制误差的工程手段有三层。第一层旋转因子表用double计算再转存成float避免在 float 精度上逐次累加角度误差第二层蝶形内部的中间乘积用double累积只在写回float时截断相当于在关键路径上做了局部双精度第三层对于需要高动态范围的场景在变换前对输入做整体归一化把幅度缩放到 [0,1] 区间防止蝶形中间结果溢出 float 上限。// 变换前对输入做功率归一化避免大信号导致中间结果溢出 float max_val 0.0f; for (int i 0; i n; i) max_val fmaxf(max_val, cabsf(buf[i])); if (max_val 0.0f) { float scale 1.0f / max_val; for (int i 0; i n; i) buf[i] * scale; }这种预缩放会影响频谱幅值所以逆变换或幅值计算时要乘回max_val。幅度谱的正确修正是第 5 章的主题。5. 频谱幅值校正与加窗的工程细节5.1 幅值修正直流、幅值与均方根的关系FFT 输出的每个频点复数值与实际信号的幅值之间有固定的换算关系。对于幅值为 A 的正弦信号若采样点数 N 且不加窗则对应频点上的幅度谱值为 A*N/2正频率处。因此实际的信号幅值需要除以 N/2 才能恢复直流分量则只需除以 N。这个换算关系在测试时最容易出错直接用原始 FFT 输出画频谱图纵轴的数字大得离谱以为程序写错了其实是没做归一化。// 单边幅度谱计算freq_bin是FFT输出的频点下标 float amplitude 2.0f * cabsf(fft_out[freq_bin]) / N; if (freq_bin 0) amplitude cabsf(fft_out[0]) / N; // 直流分量特殊处理参数说明N是 FFT 点数乘 2 是因为把负频率的功率折回了正频率直流分量没有镜像不能乘 2。若信号频率恰好落在两个频点之间非整数周期采样幅值会被稀释到相邻多个频点上这是频谱泄漏不是计算错误。5.2 加窗选择汉宁窗还是平顶窗频谱泄漏是有限长截断的固有问题加窗是一种折衷主瓣变宽换取旁瓣降低。汉宁窗Hann是最常用的通用窗主瓣宽度为 4 个频点旁瓣衰减 31dB平顶窗Flat-top则把幅值精度放在第一位主瓣展宽到 8 个频点左右但幅值读数误差可以降到 0.01dB 量级。选择逻辑很直接做频谱幅值精确测量用平顶窗做谐波分析和噪声底测量用汉宁窗。// 汉宁窗应用win_buf为与输入等长的窗口系数数组 for (int i 0; i N; i) { float w 0.5f * (1.0f - cosf(2.0f * PI * i / (N - 1))); data[i] * w; win_buf[i] w; } // 加窗后做FFT幅值换算时需要除以窗口系数的平均值 float win_gain 0.0f; for (int i 0; i N; i) win_gain win_buf[i]; win_gain / N; // 汉宁窗均值约为0.5 // 实际幅值 cabsf(fft_out[k]) * 2 / (N * win_gain);注意加窗后的归一化系数不是 N/2而是 N 乘以窗口均值。汉宁窗的均值约 0.5所以等效归一化因子是 N/4如果忘记除以win_gain所有非直流频点的幅值会整体偏小约 6dB。这个坑在频谱分析仪应用的实测阶段几乎必现。5.3 频率分辨率与频谱细化FFT 的频率分辨率是 fs/N即采样率除以 FFT 点数。N1024、fs10240Hz 时的分辨率是 10Hz。提高分辨率有两个途径增加 N需要更长的采样时间或降低 fs可能造成混叠。频谱细化Zoom FFT的思路是把感兴趣频段搬移到零频附近再降低采样率实现上可以用数字下变频实现先用一个复正弦把目标频段混频到基带再经过低通滤波最后做点数较少的 FFT。这样既保持频率分辨率又不必做超长 FFT 消耗内存。这个方法只推荐在目标频段很窄且明确时使用否则预滤波和混频的代价可能超过去做一次完整 FFT。// 数字下变频把中心频率fc附近的频谱搬移到基带 for (int i 0; i N; i) { float phase -2.0f * PI * fc * i / fs; baseband[i] data[i] * (cosf(phase) sinf(phase) * I); } // 然后对baseband做低通滤波再降采样做FFT参数说明fc是目标频段的中心频率fs是采样率混频后的信号包含搬移到基带的原信号和搬移到 2*fc 处的镜像必须经低通滤波滤除镜像否则后续抽取会产生混叠。低通滤波器的截止频率设为目标频段带宽的一半滤波完再做 M 倍抽取FFT 点数可以减少为 N/M 而分辨率不变。6. 验证你的 FFT 实现GNU Octave 对照与边界用例最有效的验证方式不是只看波形而是用一个独立实现做逐点对照。GNU Octave 是免费且语法兼容 MATLAB 的工具完全够用。验证流程分三步先在 Octave 里生成测试信号并计算参考频谱再用 C 程序读入同一组数据计算 FFT最后逐点比较幅度谱和相位谱的误差。% GNU Octave 参考实现生成两个频率分量的测试信号 fs 10240; N 1024; t (0:N-1) / fs; x 1.5 * sin(2*pi*1000*t) 0.8 * sin(2*pi*3400*t); X fft(x); dlmwrite(test_input.csv, x, precision, %.8f); dlmwrite(test_output_ref.csv, [real(X) imag(X)], precision, %.8f); % 把数据和参考频谱导出为CSV供C程序对照使用说明dlmwrite输出的 CSV 可以直接被 C 程序用fscanf读入test_output_ref.csv第一列是实部第二列是虚部。注意 Octave 的 FFT 输出顺序和 C 实现一致都是自然频率顺序不需要额外调整。对照脚本的边界用例至少要覆盖三种场景纯直流输入所有频点除 0 号外都应为 0、单音正弦幅值应精确恢复、冲激输入频谱应为平坦的常数。其中冲激测试最能暴露位反转实现的错误——任何一个交换错误都会导致频谱出现周期性伪峰且伪峰位置和交换错误的位索引直接相关。# 用awk快速比较C程序输出和Octave参考值的最大误差 awk NRFNR{re[NR]$1; im[NR]$2; next} {dr$1-re[FNR]; di$2-im[FNR]; if(dr0)dr-dr; if(di0)di-di; if(drmax)maxdr; if(dimax)maxdi} END{printf max_err %.6e\n, max} test_output_ref.csv test_output_c.csv误差上界的判断标准是浮点单精度实现的最大误差在1e-4到1e-3量级属于正常如果达到1e-2以上基本可以断定算法或数据对齐有误不是舍入误差的范畴。对照通过后可以进一步验证复数打包技巧——用两路实数信号打包做一次复数 FFT再拆分与分别做实数 FFT 的结果比较误差应在同一量级。这个测试通过后你的 FFT 实现才算真正具备实用价值。本文还有配套的精品资源点击获取