简介这份资源面向具备一定C语言与信号处理基础的学习者聚焦快速傅里叶变换FFT与逆变换IFFT的工程实现帮助读者理解离散傅里叶变换的数学原理并落地为可编译运行的代码。内容围绕Cooley-Tukey基2算法展开涵盖序列分解、蝶形运算、位反序处理以及IFFT中除以N、取共轭等关键步骤并涉及复数运算、内存管理与浮点误差等实践细节。资源包共47个文件约2.5MB以cpp源码、h头文件为核心辅以obj、pdb、ilk等编译中间产物与exe可执行文件另有sln、vcxproj等工程配置及日志文件构成一套完整的Visual Studio工程。目前已有2757人学习下载适合信号处理、图像分析与数字滤波方向的学习者参考可借助现成工程快速调试、对照理论验证算法正确性并在此基础上进行性能与精度优化。1. 从一段卡了三天的采样数据说起为什么我劝你手写一遍 FFT去年帮一个做振动信号采集的朋友看代码他拿到的是一块 24 位 ADC 连续采的 4096 点数据想提取几个关键频点的幅值。他一开始用了个现成的库结果在嵌入式板子上跑一次要十几毫秒实时性直接崩了。后来我让他把 FFT C 的实现思路用 C 语言重写一遍定点化之后单次运算压到 1.2 毫秒以内问题才算解决。这件事让我意识到很多人对 FFT 的认知停留在“调个库就行”但真到了资源受限、需要抠每一个蝶形运算的时候手写一遍的价值就出来了。这份资源就是围绕“利用 C 语言实现 FFT 运算”展开的它不是一个简单的代码片段而是一套从复数结构体定义、位反转重排、蝶形运算到结果验证的完整实现。适合两类人一类是刚学完数字信号处理、想把公式变成可运行代码的学生或转行者另一类是在 MCU、DSP 或边缘设备上做信号处理需要自己掌控内存和运算量的工程师。如果你只是想在 PC 上做个频谱图用 Python 的 numpy.fft 就够了但如果你想理解每一层循环在干什么、为什么基 2 比基 4 好写、旋转因子表怎么省内存那这份东西值得你花一个下午跟着敲一遍。2. 拆开 FFT 的黑匣子从 DFT 到蝶形运算的降维逻辑2.1 为什么不能直接算 DFTN² 到 N log N 的账先算一笔账。一个 N 点的离散傅里叶变换定义是 X[k] Σ x[n] * W_N^(nk)其中 W_N e^(-j2π/N)。如果直接按定义写双重循环复杂度是 O(N²)。N1024 时复数乘法和加法各要一百多万次N4096 时直接飙到一千六百万次以上。在 PC 上可能感觉不明显但在主频几十兆的 MCU 上这就是秒级和毫秒级的区别。FFT 的核心思路是利用旋转因子的周期性和对称性把长序列拆成短序列反复利用已经算过的结果。最经典的基 2 按时间抽取DIT算法把 N 点序列按奇偶下标分成两组每组 N/2 点然后递归下去。这样每一级有 N/2 个蝶形运算一共 log₂N 级总运算量降到 (N/2) log₂N 次复数乘法。N1024 时复数乘法从一百多万次降到 5120 次差了 200 倍。这就是为什么 FFT 是信号处理里绕不开的东西。常见做法是先把输入序列补零到 2 的整数次幂因为基 2 算法要求 N 2^m。如果你的采样点数不是 2 的幂要么补零要么用混合基算法但后者实现复杂度高很多。我一般建议在采集端就控制好点数比如 1024、2048、4096省得后面折腾。2.2 位反转重排为什么输入要“打乱”才能算基 2 DIT 算法有一个反直觉的地方它要求输入序列按位反转的顺序排列。举个例子N8 时原始下标 0 到 7 的二进制是 000、001、010、011、100、101、110、111位反转之后变成 000、100、010、110、001、101、011、111对应十进制就是 0、4、2、6、1、5、3、7。也就是说x[1] 要和 x[4] 换位置x[3] 要和 x[6] 换位置。为什么要这么做因为 DIT 算法在递归拆分时先按奇偶分再按奇偶的奇偶分最后得到的顺序正好是位反转的。如果你不做这一步直接按自然顺序输入算出来的频谱会是乱的。这个坑我见过不止一个人踩代码跑通了但频谱峰值对不上查了半天才发现是没做重排。实现位反转有两种常见方式。一种是查表法预先算好每个下标对应的反转值运行时直接查表交换速度快但占内存。另一种是原地计算法用位操作逐位反转省内存但每次要多几个指令周期。在 RAM 紧张的 MCU 上我一般用第二种下面这段代码就是原地反转的写法// 原地位反转重排n 必须是 2 的整数次幂 void bit_reverse(complex_t *x, int n) { int j 0; for (int i 0; i n - 1; i) { if (i j) { // 交换 x[i] 和 x[j] complex_t temp x[i]; x[i] x[j]; x[j] temp; } // 从高位开始找第一个 0把它变成 1后面的 1 变成 0 int k n 1; while (k j) { j - k; k 1; } j k; } }这段代码的逻辑是用 j 跟踪当前应该交换到的位置每次 i 递增时j 按位反转的规律更新。k n 1是从最高位开始试探while (k j)循环把 j 中对应位为 1 的部分减掉最后j k把第一个 0 位变成 1。整个过程不需要额外的反转表空间复杂度 O(1)。参数 n 必须是 2 的幂否则 k 的移位会出错。如果你不确定输入长度先在调用前做个判断。2.3 旋转因子表提前算好还是边算边用旋转因子 W_N^k cos(2πk/N) - j*sin(2πk/N) 在每一级蝶形运算里都要用。如果每次都用 sin/cos 函数现算在嵌入式环境里是灾难——三角函数库本身可能就占几 KB Flash而且调用开销大。常见做法是预先算一张表存 N/2 个复数运行时直接查。但这里有个取舍存全表需要 N/2 * 2 * sizeof(float) 字节。N4096 时单精度浮点要 16KB对很多小容量 MCU 来说太多了。我一般用两种省内存的策略。第一种是只存第一象限的 N/4 个值利用对称性推导其他象限内存减半。第二种是边算边用递推利用 W_N^(k1) W_N^k * W_N^1 的关系每次乘一个固定因子但这样会有累积误差点数多了之后精度会掉。下面是一个旋转因子表的生成代码用 double 算完再转 float 存减少累积误差// 预计算旋转因子表只存前 N/2 个 void init_twiddle(complex_t *twiddle, int n) { for (int i 0; i n / 2; i) { double angle -2.0 * M_PI * i / n; twiddle[i].real (float)cos(angle); twiddle[i].imag (float)sin(angle); } }参数说明twiddle 数组长度至少 N/2n 是 FFT 点数。用 double 计算再截断到 float比全程 float 的精度好一截代价只是初始化时多几个周期。如果你的平台没有 double 硬件支持那就只能用 float但建议把 N 控制在 2048 以内否则误差可能影响频谱幅值的准确性。3. 手把手写一个可运行的基 2 FFT结构体、蝶形与三层循环3.1 复数结构体与内存布局别小看对齐这件事C 语言没有内置复数类型C99 虽然加了complex.h但在很多嵌入式工具链里支持不好所以自己定义一个结构体是最稳妥的。定义本身很简单typedef struct { float real; float imag; } complex_t;但内存布局有讲究。如果结构体里两个 float 是连续存放的大多数编译器会按 4 字节对齐整个结构体占 8 字节没有填充。但如果你在结构体里混了 double 和 float或者加了其他字段就可能出现填充字节导致数组占用的内存比预期大。我一般会在定义之后加一句编译期断言确认大小// 编译期检查如果大小不是 8 字节就报错 _Static_assert(sizeof(complex_t) 8, complex_t size mismatch);这个断言在 C11 里可用C99 可以用 typedef 数组的技巧模拟。别嫌麻烦我见过因为结构体对齐问题导致 DMA 传输数据错位的案例查起来非常痛苦。另外如果你的平台支持 SIMD 或者有复数运算指令可以考虑用编译器内置的向量类型但那就不是纯 C 了。这份资源的目标是纯 C 实现所以结构体方案最通用。3.2 蝶形运算的三种写法从直观到优化蝶形运算是 FFT 的基本单元输入两个复数 a 和 b输出 a bW 和 a - bW。最直观的写法是直接展开// 最直观的蝶形运算可读性好但临时变量多 void butterfly_naive(complex_t *a, complex_t *b, complex_t w) { complex_t t; t.real b-real * w.real - b-imag * w.imag; t.imag b-real * w.imag b-imag * w.real; complex_t a_old *a; a-real a_old.real t.real; a-imag a_old.imag t.imag; b-real a_old.real - t.real; b-imag a_old.imag - t.imag; }这段代码逻辑清晰但用了三个临时复数编译器优化后可能还是多几条指令。第二种写法是减少临时变量把中间结果直接算到目标位置// 减少临时变量的蝶形运算 void butterfly_compact(complex_t *a, complex_t *b, complex_t w) { float tr b-real * w.real - b-imag * w.imag; float ti b-real * w.imag b-imag * w.real; float ar a-real; float ai a-imag; a-real ar tr; a-imag ai ti; b-real ar - tr; b-imag ai - ti; }这种写法只用了四个 float 临时变量在寄存器充足的平台上能全部放进寄存器速度更快。第三种是定点化写法把 float 换成 int32_t 或 int16_t配合 Q 格式。定点化的好处是速度快、功耗低但需要仔细处理溢出和精度。我一般建议先在浮点上验证算法正确再根据平台需求决定是否定点化。参数说明w 是当前级的旋转因子a 和 b 是参与蝶形运算的两个复数指针。注意 b 在运算后会被覆盖为差值所以如果后续还要用原始 b 的值需要提前备份。3.3 三层循环的完整实现与调用示例把位反转、旋转因子和蝶形运算串起来就是完整的 FFT 函数。核心是一个三层循环最外层遍历 log₂N 个级中间层遍历每一级的 N/2 个蝶形最内层做具体的蝶形运算。下面是一个可编译运行的完整实现#include stdio.h #include math.h #include stdint.h #define PI 3.14159265358979323846 typedef struct { float real; float imag; } complex_t; // 位反转重排 void bit_reverse(complex_t *x, int n) { int j 0; for (int i 0; i n - 1; i) { if (i j) { complex_t temp x[i]; x[i] x[j]; x[j] temp; } int k n 1; while (k j) { j - k; k 1; } j k; } } // 基 2 DIT FFTn 必须是 2 的幂 void fft(complex_t *x, int n) { bit_reverse(x, n); // 预计算旋转因子表 complex_t twiddle[n / 2]; for (int i 0; i n / 2; i) { double angle -2.0 * PI * i / n; twiddle[i].real (float)cos(angle); twiddle[i].imag (float)sin(angle); } // 三层循环 for (int stage 1; stage (int)log2(n); stage) { int half 1 (stage - 1); // 当前级蝶形的半跨度 int step 1 stage; // 当前级的分组跨度 for (int group 0; group n; group step) { for (int k 0; k half; k) { int idx_a group k; int idx_b group k half; // 旋转因子下标k * (n / step) int tw_idx k * (n / step); complex_t w twiddle[tw_idx]; float tr x[idx_b].real * w.real - x[idx_b].imag * w.imag; float ti x[idx_b].real * w.imag x[idx_b].imag * w.real; float ar x[idx_a].real; float ai x[idx_a].imag; x[idx_a].real ar tr; x[idx_a].imag ai ti; x[idx_b].real ar - tr; x[idx_b].imag ai - ti; } } } } // 测试生成一个单频正弦波看频谱峰值位置 int main() { int n 1024; complex_t data[n]; float fs 1024.0f; // 采样率 float f0 128.0f; // 信号频率 for (int i 0; i n; i) { data[i].real cosf(2.0f * PI * f0 * i / fs); data[i].imag 0.0f; } fft(data, n); // 找最大幅值对应的下标 float max_mag 0.0f; int max_idx 0; for (int i 0; i n / 2; i) { float mag sqrtf(data[i].real * data[i].real data[i].imag * data[i].imag); if (mag max_mag) { max_mag mag; max_idx i; } } printf(峰值下标: %d, 对应频率: %.2f Hz, 幅值: %.2f\n, max_idx, max_idx * fs / n, max_mag); return 0; }这段代码可以直接用 gcc 编译gcc fft.c -o fft -lm。运行后应该输出峰值下标 128对应频率 128 Hz幅值约 512。如果输出不对先检查 n 是不是 2 的幂再看旋转因子的下标计算有没有越界。参数说明stage从 1 到 log₂Nhalf是当前级蝶形运算中两个输入点的距离step是分组的跨度。tw_idx k * (n / step)这个公式是关键它保证了每一级用的旋转因子是正确的。我见过有人在这里直接写tw_idx k结果只有第一级对后面全错。4. 避坑与排查那些让我熬夜的 FFT 翻车现场4.1 频谱峰值对不上先查位反转再查旋转因子现象代码编译运行都没报错但频谱峰值出现在错误的位置或者幅值明显偏小。原因最常见的是位反转没做或者旋转因子的下标算错了。解决先用一个已知频率的正弦波做输入比如 N1024、fs1024、f0128理论上峰值应该在 128。如果峰值在别的地方先注释掉位反转看结果是不是更乱如果更乱说明位反转逻辑有问题。再检查tw_idx的计算确保每一级的旋转因子下标是k * (n / step)而不是k。4.2 浮点精度累积点数越大越明显现象N4096 时频谱底噪比 N1024 时高很多弱信号被淹没。原因旋转因子用 float 存储蝶形运算反复乘加误差累积。解决旋转因子表用 double 计算再转 float或者干脆全程用 double。如果平台支持可以在关键累加步骤用 Kahan 求和。另一个办法是减少级数用基 4 算法但实现复杂度会上升。我一般建议 N 不超过 4096 时用 float 够用超过就上 double。4.3 内存越界旋转因子表长度不够现象程序在 PC 上跑得好好的移植到嵌入式平台就 HardFault。原因旋转因子表只开了 N/2 的长度但某些级的tw_idx可能达到 N/2 - 1如果数组边界没对齐或者编译器优化导致越界访问就会触发硬件异常。解决把 twiddle 数组长度开到 N/2 1或者在访问前加断言assert(tw_idx n / 2)。另外如果 twiddle 是局部数组N 很大时可能爆栈建议改成静态数组或全局数组。4.4 输入不是 2 的幂补零还是换算法现象采样点数 1000直接调用 fft 函数后结果完全不对。原因基 2 算法要求 N 2^m1000 不是 2 的幂位反转和级数计算都会出错。解决最简单的办法是补零到 1024但补零会改变频谱的幅值需要在结果里做补偿。另一种是改用混合基算法支持 N 2^a * 3^b * 5^c但代码量会大很多。我一般建议在采集端就控制点数实在不行就补零补零后的幅值乘以 N_original / N_padded 做修正。4.5 定点化溢出Q 格式选择不当现象定点化之后某些频点的幅值变成负数或者饱和到最大值。原因Q 格式的整数位不够蝶形运算的中间结果溢出。解决先估算输入信号的最大幅值再确定 Q 格式的小数位数。比如输入幅值不超过 1.0用 Q15 时整数位只有 1 位蝶形运算的加法可能溢出。可以改用 Q14 或者 Q13牺牲一点精度换动态范围。另外每一级运算后可以做一次右移防止累积溢出但会引入截断误差。5. 进阶技巧用查表法把旋转因子内存砍掉一半5.1 对称性利用只存 N/4 个值旋转因子有一个很好的对称性W_N^(k N/4) -j * W_N^kW_N^(N/2 - k) conj(W_N^k)。这意味着你只需要存前 N/4 个值其他都可以通过取负、交换实虚部、取共轭得到。具体来说对于下标 k 在 [0, N/4) 范围内直接查表对于 [N/4, N/2)用 -j * table[k - N/4]对于 [N/2, 3N/4)用 -table[k - N/2]对于 [3N/4, N)用 j * table[k - 3N/4]。这样内存直接减半N4096 时从 16KB 降到 8KB。实现的时候我一般写一个内联函数来取旋转因子// 从 N/4 长度的表中取第 k 个旋转因子 static inline complex_t get_twiddle(const complex_t *table, int k, int n) { int quarter n / 4; if (k quarter) { return table[k]; } else if (k 2 * quarter) { complex_t t table[k - quarter]; // 乘以 -j实部变虚部虚部变负实部 complex_t r { t.imag, -t.real }; return r; } else if (k 3 * quarter) { complex_t t table[k - 2 * quarter]; complex_t r { -t.real, -t.imag }; return r; } else { complex_t t table[k - 3 * quarter]; complex_t r { -t.imag, t.real }; return r; } }参数说明table 长度是 N/4k 是旋转因子下标范围 [0, N/2)。n 是 FFT 点数。这个函数在每一级的内层循环里调用如果编译器能内联开销很小。注意 quarter 的计算要在循环外做好别每次都算。5.2 验证方法用已知信号做端到端对拍写完优化版之后怎么确认结果和原始版本一致我一般用两种方法。第一种是生成一个多频正弦叠加信号分别用原始版本和优化版本跑把输出逐点相减看最大误差是否在 1e-5 以内。第二种是跟 numpy.fft 对拍把 C 程序的输出存成文本用 Python 读进来做对比。下面是一个简单的对拍脚本import numpy as np # 读取 C 程序输出的频谱 c_result np.loadtxt(c_fft_output.txt, dtypenp.complex64) # 生成同样的输入信号 n 1024 fs 1024.0 f0 128.0 t np.arange(n) signal np.cos(2 * np.pi * f0 * t / fs) # numpy 的 FFT np_result np.fft.fft(signal) # 比较前 N/2 个点 diff np.abs(c_result[:n//2] - np_result[:n//2]) print(f最大误差: {diff.max():.6f}) print(f平均误差: {diff.mean():.6f})如果最大误差在 1e-4 量级说明浮点精度没问题如果误差很大先检查两边用的输入信号是否完全一致再检查 C 程序的输出有没有做归一化。numpy 的 fft 不做归一化C 程序也不做所以可以直接比。5.3 一个习惯每次改完都跑一遍单频测试从那以后我每次改 FFT 代码不管是换旋转因子表还是改蝶形运算都强制走一遍单频测试生成一个 128 Hz 的正弦波跑完看峰值是不是在 128幅值是不是接近 N/2。这个习惯帮我省了很多调试时间因为单频测试对位反转、旋转因子、蝶形运算的错误都很敏感一旦哪里写错峰值立刻跑偏。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取