C语言手写FFT：从DFT到O(N log N)的实战指南-北京尧图网络

简介这份资源面向信号处理、图像分析与数字滤波方向的学习者与开发者提供用C语言实现快速傅里叶变换FFT及逆变换IFFT的完整工程代码。内容围绕Cooley-Tukey基2算法展开涵盖序列分解、蝶形运算、位反序处理等关键环节并给出IFFT的共轭与归一化实现思路适合具备一定C/C基础、希望理解频域算法底层原理的读者参考。压缩包共47个文件约2.5MB以cpp源码、h头文件为核心辅以obj、pdb、ilk等编译中间产物及sln、vcxproj工程配置另有exe可执行文件与日志记录方便直接编译运行与调试。目前已有2757人学习下载。借助这套代码读者可对照复数运算与内存管理细节梳理FFT/IFFT的完整实现流程并在此基础上开展性能测试与精度验证为后续在信号处理、图像分析等场景中的算法移植与优化提供可复用的工程参考。1. 用 C 语言手写 FFT从 DFT 的 O(N²) 到 O(N log N) 到底省在哪如果你写过音频频谱显示、电机振动分析或者简易示波器上位机大概率绕不开一个坎采样 1024 个点做离散傅里叶变换DFT双重循环跑下来要一百多万次复数乘法在单片机上直接卡成幻灯片。这时候老工程师会告诉你上 FFT。FFT 不是一种新的变换它只是 DFT 的快速算法核心思路是把长序列按奇偶拆成两半利用旋转因子的周期性和对称性把 O(N²) 的乘法量压到 O(N log N)。1024 点 DFT 约 104 万次复数乘法基 2 FFT 只要约 5120 次差了 200 倍这就是为什么它值得你亲手用 C 写一遍。这篇笔记面向三类人正在学数字信号处理、想把公式变成能跑的代码的学生需要在 MCU 或 DSP 上做实时频谱、又不想背一整个库的嵌入式开发者以及想搞明白「库函数里那个 fft 到底干了什么」的工程师。我会从复数结构体和旋转因子表讲起给出可编译的基 2 时域抽取DIT实现再补上位倒序、窗函数、实数序列优化这些实战里躲不掉的东西。全程只依赖标准 C不引第三方库方便你直接搬进工程。2. 基 2 DIT-FFT 的推导与最小可运行实现2.1 为什么是「奇偶拆分」而不是别的拆法DFT 的定义是 X[k] Σₙ₌₀^{N-1} x[n]·W_N^{nk}其中 W_N e^{-j2π/N} 叫旋转因子。当 N 是 2 的整数次幂时可以把 x[n] 按 n 的奇偶分成两组偶数组 x[2r] 和奇数组 x[2r1]各 N/2 个点。代入后得到X[k] E[k] W_N^k · O[k] X[k N/2] E[k] − W_N^k · O[k]其中 E[k]、O[k] 分别是偶、奇子序列的 N/2 点 DFT。这两行是整个 FFT 的命门一个 N 点变换被拆成两个 N/2 点变换再花一次复数乘加合并。递归下去log₂N 层之后只剩 2 点变换蝴蝶运算就成型了。这就是时域抽取Decimation In Time名字的由来——在时域上不断抽取。选基 2 而不是基 4 或分裂基理由很实际基 2 代码最短、最容易验证N 必须是 2 的幂这个限制在工程里几乎不是问题采样点数取 256、512、1024、2048 都是常规操作。基 4 能再省约 25% 的乘法但代码复杂度和调试成本对多数人来说不划算除非你在做极致的定点优化。2.2 复数结构体与旋转因子表先把基础设施搭好。C 语言没有内置复数C99 的 complex.h 在不少嵌入式工具链里支持不全自己定义一个最省事。#include stdio.h #include math.h #ifndef M_PI #define M_PI 3.14159265358979323846 #endif typedef struct { double re; // 实部 double im; // 虚部 } complex_t; // 复数乘法: (abi)(cdi) (ac-bd) (adbc)i static inline complex_t c_mul(complex_t a, complex_t b) { complex_t r; r.re a.re * b.re - a.im * b.im; r.im a.re * b.im a.im * b.re; return r; } static inline complex_t c_add(complex_t a, complex_t b) { complex_t r { a.re b.re, a.im b.im }; return r; } static inline complex_t c_sub(complex_t a, complex_t b) { complex_t r { a.re - b.re, a.im - b.im }; return r; }旋转因子表预先算好别在蝴蝶循环里现算 cos/sin。三角函数调用一次几十到上百个周期放在内层循环里能把 FFT 的性能优势全吃掉。表长取 N/2 就够因为 W_N^{kN/2} −W_N^k对称性可以省一半存储。// 预计算旋转因子 W_N^k exp(-j*2*pi*k/N), k 0..N/2-1 static void init_twiddle(complex_t *w, int n) { for (int k 0; k n / 2; k) { double angle -2.0 * M_PI * k / n; w[k].re cos(angle); w[k].im sin(angle); } }参数说明n 是 FFT 点数必须是 2 的幂w 数组长度至少 n/2。用 double 是为了先把算法跑对定点版本后面再谈。如果你在 Cortex-M4 这类带 FPU 的芯片上float 版本速度也够看没有 FPU 的 M0/M3就得考虑 Q15 定点那是另一个话题。2.3 位倒序FFT 里最容易翻车的一步递归拆分的输入顺序是「偶偶偶…奇奇奇…」等价于把下标按二进制位翻转。比如 N8 时x[1]001要换到 x[4]100的位置。这一步写错输出频谱会整个乱掉而且不会报错属于典型的黑匣子故障。// 原地位倒序重排 static void bit_reverse(complex_t *x, int n) { int j 0; for (int i 1; i n; i) { int bit n 1; // 从高位开始找第一个可以进位的 0 for (; j bit; bit 1) { j ^ bit; } j ^ bit; if (i j) { complex_t tmp x[i]; x[i] x[j]; x[j] tmp; } } }逻辑说明j 始终维护 i 的位反转值。内层循环把 j 里从高到低连续的 1 清零相当于二进制加法进位再把第一个 0 置 1。i j 的判断保证每对只交换一次避免换回去。这个写法比「逐位拼反转值」的版本少一层循环是常见做法。2.4 三层循环的蝴蝶主体// 原地基2 DIT-FFTx 长度 nw 为预计算旋转因子 void fft(complex_t *x, int n, const complex_t *w) { bit_reverse(x, n); // len: 当前合并的子变换长度, 从 2 到 n 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 k 0; k half; k) { complex_t u x[i k]; complex_t v c_mul(x[i k half], w[k * step]); x[i k] c_add(u, v); x[i k half] c_sub(u, v); } } } }参数说明len 是当前层的蝶形跨度每层翻倍half 是半跨度step n/len 是旋转因子抽取间隔因为第 len 点变换用的因子是 W_len^k W_n^{k·n/len}。u 是上半支v 是下半支乘完旋转因子的结果一加一减就是标准蝶形。整个函数原地操作不需要额外输出数组内存友好。2.5 一个能直接编译验证的 mainint main(void) { const int N 8; complex_t x[8] { {1,0},{2,0},{3,0},{4,0},{5,0},{6,0},{7,0},{8,0} }; complex_t w[4]; init_twiddle(w, N); fft(x, N, w); for (int i 0; i N; i) { printf(X[%d] %.4f %.4f j\n, i, x[i].re, x[i].im); } return 0; }编译命令gcc fft.c -o fft -lm注意-lm不能漏否则链接 cos/sin 会报错。验证方法输入 1 到 8 的实序列X[0] 应该等于 36直流分量即所有样本之和其余点按共轭对称分布。如果 X[0] 不对先查位倒序如果只有个别点不对查旋转因子步进 step。3. 从能跑到能用实数输入、窗函数与频谱换算3.1 实数序列别硬套复数 FFT工程里 90% 的输入是实数ADC 采样、音频 PCM直接当复数喂进去虚部全填 0等于白算一半。常见做法有两种一是用 N/2 点复数 FFT 打包两个实数序列适合批量处理二是利用实序列频谱的共轭对称性只算前 N/21 个点。对单路信号最省事的是照常跑 N 点复数 FFT然后只取前 N/2 个输出后一半是镜像扔掉即可。// 实数输入虚部清零后调用 fft只取前 N/21 个点 void fft_real(const double *in, complex_t *out, int n, const complex_t *w) { for (int i 0; i n; i) { out[i].re in[i]; out[i].im 0.0; } fft(out, n, w); // out[0..n/2] 有效out[n/21..n-1] 是共轭镜像 }3.2 加窗不加窗的频谱泄漏能骗你半天直接对截断信号做 FFT如果信号频率不是 bin 频率的整数倍能量会泄漏到邻近 bin谱线变胖、幅度偏低。这不是 FFT 的 bug是有限长度截断的必然结果。解决办法是加窗把两端压到接近零。窗类型主瓣宽度旁瓣衰减适用场景矩形窗最窄−13 dB整周期采样、瞬态分析Hann 窗中等−31 dB通用频谱分析Hamming 窗中等−43 dB语音、需要低旁瓣Blackman 窗宽−58 dB强干扰下的弱信号检测Hann 窗生成很简单w[n] 0.5·(1 − cos(2πn/(N−1)))。加窗后幅度要乘补偿系数Hann 约 2.0否则测出来的幅值偏小一半这是新手最常踩的坑之一。void apply_hann(double *x, int n) { for (int i 0; i n; i) { x[i] * 0.5 * (1.0 - cos(2.0 * M_PI * i / (n - 1))); } }3.3 幅度谱和 dB 换算FFT 输出是复数工程上要看的是幅度谱mag[k] sqrt(re² im²)。要还原真实物理幅值还得除以 N单边谱再乘 2直流和 Nyquist 点不乘。转 dB 用 20·log10(mag/ref)ref 取满量程或 1.0 看你标定方式。// 计算单边幅度谱n 点输入输出 n/21 个幅值 void magnitude_spectrum(const complex_t *spec, double *mag, int n) { mag[0] sqrt(spec[0].re * spec[0].re spec[0].im * spec[0].im) / n; for (int k 1; k n / 2; k) { double m sqrt(spec[k].re * spec[k].re spec[k].im * spec[k].im); mag[k] 2.0 * m / n; // 单边谱补偿 } mag[n / 2] sqrt(spec[n/2].re * spec[n/2].re spec[n/2].im * spec[n/2].im) / n; }频率轴换算f[k] k · fs / nfs 是采样率。比如 fs 8000 Hz、n 1024bin 间隔约 7.8 Hz。想提高频率分辨率只有两条路加长采样时间增大 n或降低采样率没有第三条。4. 避坑与排查FFT 调试里最常见的 5 个翻车现场4.1 输出全是 NaN 或极大值现象跑完 FFT打印出来一堆 nan 或者 1e300 量级的数。原因通常是旋转因子表没初始化就用了或者数组越界写坏了内存。解决在 fft 入口加断言检查 n 是否为 2 的幂(n (n-1)) 0确认 init_twiddle 在 fft 之前调用w 数组长度不小于 n/2。用 valgrind 或编译时开-fsanitizeaddress能快速定位越界。4.2 频谱左右不对称实数输入时现象实数信号做完 FFT前一半和后一半不是共轭对称。原因多半是位倒序写错或者旋转因子步进 step 算错。排查方法先用 N4 的极小案例手算对比输入 {1,0,0,0}正确输出应该全是 1。如果不对把 bit_reverse 单独拎出来打印下标映射表一眼就能看出问题。4.3 峰值频率差一个 bin现象输入 1 kHz 正弦频谱峰值却落在 992 Hz 或 1008 Hz。这是频谱泄漏加栅栏效应不是代码错。解决要么让采样点数是信号周期的整数倍同步采样要么加 Hann 窗后用插值法如三点抛物线插值修正峰值位置。指望 FFT 直接给出精确频率在非整周期采样下是不现实的。4.4 幅值随 N 变化现象同样的信号N 从 1024 改成 2048幅度谱数值变了。原因是你忘了归一化或者单边谱补偿系数用错。记住幅度谱必须除以 N单边谱非直流点乘 2。加窗后还要除以窗函数的相干增益Hann 为 0.5否则幅值系统性偏低。4.5 在 MCU 上跑一次要几十毫秒现象1024 点 FFT 在 72 MHz 的 M3 上耗时超过 50 ms实时性崩了。原因通常是用了 double、旋转因子现算、或者没开编译器优化。解决换 float、预计算旋转因子表放 Flash、编译开-O2能砍掉一大半时间。如果还不够考虑 Q15 定点或改用 CMSIS-DSP 这类针对内核优化过的库——但前提是你已经理解了自己手写版本的每一行否则出了问题连查都不知道从哪查。5. 进阶技巧用 N/2 点复数 FFT 一次算两路实数当你需要同时处理两路实数信号比如 I/Q 两路或者立体声左右声道有个经典技巧把两路实数打包成一个复数序列做一次 N 点复数 FFT再用对称性分离出两路各自的频谱。这样一次变换顶两次用计算量直接减半。设两路实信号为 a[n]、b[n]构造 z[n] a[n] j·b[n]对 z 做 N 点 FFT 得 Z[k]。利用 Z[k] 和 Z[N−k] 的共轭关系可以还原// 从 Z[k] 分离出 A[k] 和 B[k]n 点k 取 0..n/2 void split_real_spectra(const complex_t *Z, complex_t *A, complex_t *B, int n) { for (int k 0; k n / 2; k) { int mk (n - k) % n; // 镜像下标 complex_t Zk Z[k]; complex_t Zmk Z[mk]; // 共轭 Z[N-k] complex_t Zc { Zmk.re, -Zmk.im }; // A[k] (Z[k] conj(Z[N-k])) / 2 A[k].re 0.5 * (Zk.re Zc.re); A[k].im 0.5 * (Zk.im Zc.im); // B[k] (Z[k] - conj(Z[N-k])) / (2j) complex_t diff c_sub(Zk, Zc); B[k].re 0.5 * diff.im; B[k].im -0.5 * diff.re; } }逻辑说明Z[k] 的实部来自 a 的频谱、虚部来自 b 的频谱但两者在频域混叠。利用实序列频谱的共轭对称性 X[N−k] conj(X[k])联立方程就能解出 A[k] 和 B[k]。参数上 k 只需算到 n/2后面是镜像。这个技巧在雷达 I/Q 解调、双通道音频分析里非常常见省下的一次 FFT 在资源紧张的场合就是能不能实时跑的区别。验证方法造两路已知正弦 a[n]cos(2π·3n/N)、b[n]cos(2π·7n/N)跑一遍分离看 A 的峰值是否落在 bin 3、B 落在 bin 7。如果串了检查共轭那一步的虚部符号——这是最容易写反的地方我自己就因为这个符号调了一下午。最后一个习惯每次改完 FFT 代码先用 N8 的已知序列跑一遍和手算或 Python 的 numpy.fft 对一遍确认无误再上真实数据。FFT 的 bug 不会崩溃只会悄悄给你一条错的谱线等你在产品上发现的时候返工成本就不是一下午了。希望帮到你。本文还有配套的精品资源点击获取

C语言手写FFT：从DFT到O(N log N)的实战指南

相关推荐

易语言大漠模块封装与免注册调用实战指南

数字孪生工程落地：三层解耦与四步校验实战指南

从GcRoot出发：彻底搞懂Handler内存泄漏的完整引用链

更多资讯

ASP.NET MVC工作流源码深度解析与生产避坑指南

从95%到8%：国产AI芯片如何翻越CUDA生态大山

分布式存储系统设计：一致性哈希、副本放置与Raft实践

STM32驱动BH1750光照传感器实战：I²C寄存器配置与读数为0的六大根源

【C++】 priority_queue知识点总结

Intel HD630非DCH驱动26.20.100.6573稳定适配指南

今日更新

小波包能量谱原理与MATLAB实现：从频带分解到故障特征提取

基于LSTM的网易云音乐评论情感分析：数据清洗与模型调优实战

Alpha-Beta剪枝与极大极小值算法：博弈树搜索原理与实战优化

本周精选

交通拥堵预测毕设：基于Python+LSTM的车辆流量预测系统全流程解析

MATLAB/Simulink搭建10机39节点电力系统暂态稳定仿真实战

无人机航拍三维重建全流程：从SfM到网格生成的避坑指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？