简介一份基于C语言的FFT快速傅里叶变换算法实现与DSP工程示例面向数字信号处理初学者、嵌入式开发者及需要掌握基2 FFT与蝶形运算原理的学习者。压缩包内共9个文件以CCS工程文件.pjt、.cmd、.lkf和C源码.c为核心配合日志、配置及说明文档构成完整的DSP实验工程可帮助理解从算法推导到实际编译调试的全过程。包体仅5KB结构紧凑方便下载后直接查看代码与工程配置。已有203人学习下载。通过学习这份资源读者能获得可运行的FFT C语言代码、DSP开发环境下的工程组织方式以及位反转、蝶形操作、时间抽取等关键步骤的代码化呈现适合对照教材进行实验验证也可作为自行编写频谱分析、滤波器设计等应用的起点。1. FFT算法C语言实现为什么今天还要自己写一遍拿到一个叫 FFT.rar 的压缩包里面通常是 FFT 算法的 C 语言源码、测试数据和一份简短的说明文档。这类资源在网络上一搜一大把但真正打开后能直接用在项目里的很少有的依赖特定编译器有的只实现了基 2 时间抽取但没说清楚输入序列该怎么排有的旋转因子算出来精度不够在 1024 点以上跑出来的频谱和 MATLAB 对不上。问题不在算法本身而在于 C 语言实现里那些容易被忽略的细节——位反转、旋转因子预计算、定点和浮点的取舍。FFT 不是新东西但 FFT 算法 C 这个检索组合说明很多人正在做的是把算法落到具体硬件或工程里嵌入式 DSP、音频处理、振动分析、电力谐波检测。用 Python 调库谁都会可一旦要移植到 STM32、移植到没有数学库的裸机环境或者要在大循环里反复调用几千次就绕不开 C 语言实现。写这篇文章的目的就是给出一个可直接复现的基 2 FFT 实现讲清楚每一处设计的理由和参数边界同时把工程里最常见的几个坑——复数数组布局、旋转因子累乘误差、输入输出正反变换的差异——逐一拆开。2. 先从 DFT 到 FFTC语言实现前必须明确的三个数学前提2.1 为什么快速傅里叶变换能拆成蝶形运算离散傅里叶变换的定义是X[k] Σ x[n] * W_N^(nk)直接计算 N 点 DFT 需要 N^2 次复数乘法和 N(N-1) 次复数加法。当 N1024 时这个量级在单核 MCU 上不可接受。FFT 的核心思路是利用旋转因子W_N^(nk)的周期性和对称性把 N 点 DFT 拆解成两个 N/2 点 DFT再递归下去。最终每一级只需要 N/2 次复数乘法总计算量降为 N/2 * log2(N) 次。这就是基 2 时间抽取DIT算法的数学基础。// 蝶形运算的核心一次复数乘加 // 输入 a, b 为复数w 为旋转因子 // 输出 a a w*b, b a - w*b这一组运算是所有 FFT 实现的最小单元。C 语言里不直接支持复数类型除非用 C99 的complex.h所以工程上常见的做法是用两个浮点数组分别存实部和虚部或者定义一个struct { float re; float im; }。前者在缓存利用上更友好后者在代码可读性上更直观。之后展开的实现采用实部虚部分离的双数组方案原因在于它可以方便地扩展为定点整数实现——只需要把float换成int32_t并调整旋转因子的定点格式。2.2 位反转排序输入序列必须先重排基 2 DIT-FFT 要求输入序列按照二进制位反转的顺序排列。例如 N8 时自然序 0,1,2,3,4,5,6,7 对应的位反转序为 0,4,2,6,1,5,3,7。如果跳过这一步输出频谱的频率索引会完全错乱。// 位反转将 idx 按 bit_num 位宽反转 unsigned int bit_reverse(unsigned int idx, unsigned int bit_num) { unsigned int rev 0; for (unsigned int i 0; i bit_num; i) { rev (rev 1) | (idx 1); idx 1; } return rev; }参数bit_num由 N 决定等于log2(N)。N 必须是 2 的整数次幂否则该函数没有意义。在 ARM Cortex-M 上这个循环可以用RBIT指令一条搞定但可移植的 C 代码就只能用循环实现。若 N 固定建议把位反转表预先算好存入常量数组省去每次变换前的计算开销。2.3 旋转因子怎么算才能保证精度旋转因子W_N^k cos(2πk/N) - j*sin(2πk/N)。C 语言里标准做法是用cos()和sin()逐项计算但每次调用数学库函数耗时较大。常见优化是只计算前 N/2 个角度利用对称性补全其余值更进一步的做法是查表将浮点结果转成float或int16_t存表。我一般会区分两种情况如果 FFT 点数固定且变换频率很高例如音频频谱显示就预计算整张旋转因子表如果点数可变则采用半表加对称映射的方法每级只取需要的项。注意旋转因子表必须是const存储放在 Flash 里而非 RAM嵌入式环境下这一点直接决定能否跑起来 4096 点变换。方案内存占用 (N1024)适用场景全表 float1024 * 8 8KB点数固定、追求速度半表 float512 * 8 4KB点数可变、折中全表 int16 定点1024 * 4 4KB无 FPU 的 MCU实时计算0但每次调 sin/cos单次变换、不频繁旋转因子精度不足会在多级级联后产生误差累积具体表现是频谱峰值幅度偏低、非峰值处出现本底噪声抬高。后文会在精度分析部分专门讨论。3. 用 C 语言写出可复用的 FFT 核心代码3.1 复数运算的最小封装C 语言实现 FFT不建复杂的抽象层只做一组静态内联函数处理复数乘法避免重复代码。一个复数乘法需要四次实数乘法和两次加减法。// 复数乘法: (ar ai*j) * (br bi*j) // 常规实现耗时 4 次乘法 2 次加法 static inline void complex_mul(float ar, float ai, float br, float bi, float *rr, float *ri) { *rr ar * br - ai * bi; *ri ar * bi ai * br; }这段代码本身不复杂但需要留意它被调用多少次。N1024 时蝶形总数为N/2 * log2(N) 5120每个蝶形一次复数乘法五次解析下来要调用上万次。如果编译器开了-O2这种static inline写法会被完全展开消除函数调用开销。工程中谨慎将复数类型定义成结构体并频繁传值因为结构体传参会引入额外的内存拷贝在 ARM 硬浮点环境下让性能打折扣。3.2 按时间抽取的迭代式实现标准基 2 DIT-FFT 有三个嵌套循环外层是级数stage中间层是每一级内的蝶形组block内层是组内的蝶形k。每级步长step 1 stage每块大小block_size step 1。void fft_c(float *re, float *im, unsigned int n) { // 前提: n 为 2 的整数次幂 unsigned int bit_num 0; while ((1U bit_num) n) bit_num; // Step 1: 位反转重排输入 for (unsigned int i 0; i n; i) { unsigned int j bit_reverse(i, bit_num); if (j i) { float t re[i]; re[i] re[j]; re[j] t; t im[i]; im[i] im[j]; im[j] t; } } // Step 2: 多级蝶形运算 for (unsigned int stage 1; stage bit_num; stage) { unsigned int step 1U stage; // 当前级蝶形跨度 unsigned int half step 1; // 蝶形两点距离 for (unsigned int group 0; group n; group step) { for (unsigned int k 0; k half; k) { // 旋转因子: W_N^k, 但需按级缩放角度 float angle -2.0f * PI * k / (float)step; float wr cosf(angle); float wi sinf(angle); unsigned int idx_a group k; unsigned int idx_b idx_a half; float ar re[idx_a], ai im[idx_a]; float br re[idx_b], bi im[idx_b]; // 蝶形计算 float tr wr * br - wi * bi; float ti wr * bi wi * br; re[idx_a] ar tr; im[idx_a] ai ti; re[idx_b] ar - tr; im[idx_b] ai - ti; } } } }逻辑拆解如下。最外层的stage循环从 1 跑到bit_num每轮处理当前级的所有蝶形。step的含义是本级两个输入点之间的距离half是蝶形两个输出端的间隔。内层k循环计算旋转因子时以step为周期而不是以总点数 N——这一点是初学最容易写错的地方。idx_a与idx_b相差half对应蝶形图里的上下两支。这段代码可以工作但旋转因子在每级内反复调用cosf和sinf性能较差。下一节会换成查表法。3.3 预计算旋转因子表的改进版本为消除三角函数调用将旋转因子按 N 点完整预计算一次后续每级按索引跳跃取用。基 2 DIT 的第stage级需要的旋转因子索引是k * (N stage)。static float *g_wr, *g_wi; // 长度为 n/2 的旋转因子表 void fft_init(unsigned int n) { // n/2 个表项覆盖 W_N^0 到 W_N^(n/2-1) g_wr (float *)malloc(sizeof(float) * (n / 2)); g_wi (float *)malloc(sizeof(float) * (n / 2)); for (unsigned int i 0; i n / 2; i) { float angle -2.0f * PI * i / (float)n; g_wr[i] cosf(angle); g_wi[i] sinf(angle); } } // 蝶形内取旋转因子改为查表 // 第 stage 级 (1-based) 的步长为 step2^stage // 索引 k 对应的表项下标为 k * (n stage)这里的关键参数是表长。N 点 FFT 的旋转因子具有对称性完整表只需要 N/2 项即可覆盖所有蝶形。每级实际用到的索引是均匀取点——第 1 级只用W_N^0第 2 级用W_N^0和W_N^(N/4)第 3 级用 4 个点依此类推。这使得越靠前的级扫描的表项越稀疏但查表的开销远小于函数调用。实测一段 1024 点 FFT 从实时计算旋转因子改为查表后在无 FPU 的 MCU 上加速比可达 3 倍以上在桌面 CPU 上也有 20% 左右的提升。3.4 和 MATLAB/NumPy 输出对不齐怎么排查C 语言的 FFT 输出和 MATLAB 对不齐最常见的原因有三个缩放因子MATLAB 的fft()默认不除以 N而一些 C 语言库如某些 DSP 库会在输出时自动除以 N或提供正反变换归一化选项。旋转因子符号DFT 定义里指数项是负号即W_N e^(-j2π/N)。若实现中旋转因子角度写成了正号频谱会呈镜像翻转幅值不变但相位相反。位反转遗漏或反转方式错误DIT 需要输入重排DIF频率抽取需要输出重排。混用两种方式的重排逻辑会得出完全错乱的频谱。排查技巧是输入一个单频正弦波频率设为F_s / N的整数倍比如采样率 1024 Hz、N1024 时输入 2 Hz 正弦波。理论上输出频谱应在 k2 处有一个峰值。如果峰值出现在 k1022说明旋转因子符号反了如果频谱看起来像噪声散布几乎可以断定位反转环节有问题。4. FFT算法在频谱分析中的实操与参数边界4.1 从时域采到频域幅值谱的完整流程工程上做频谱分析不是直接喂原始采样数据给 FFT 就算完。流程是采样 → 去直流 → 加窗 → FFT → 幅值修正 → 频率轴映射。每一步处理不好得到的频谱都无法反映真实信号。// 以 256 点为例的完整流程 #define N 256 float sample_re[N], sample_im[N]; // 1. 采集信号填入 sample_re虚部清零 // 2. 去直流: 减去均值 float mean 0.0f; for (int i 0; i N; i) mean sample_re[i]; mean / N; for (int i 0; i N; i) sample_re[i] - mean; // 3. 加汉宁窗非矩形窗 for (int i 0; i N; i) { float w 0.5f * (1.0f - cosf(2.0f * PI * i / (N - 1))); sample_re[i] * w; } // 4. FFT fft_c(sample_re, sample_im, N); // 5. 幅值谱修正: 单边谱除 N非矩形窗还需除窗增益 float amp_spectrum[N/2]; float coh_gain 0.0f; // 窗幅度增益 (矩形窗为 1.0, 汉宁窗为 0.5) for (int i 0; i N; i) coh_gain 0.5f * (1.0f - cosf(2.0f * PI * i / (N - 1))); coh_gain / N; for (int k 0; k N/2; k) { float mag sqrtf(sample_re[k]*sample_re[k] sample_im[k]*sample_im[k]); amp_spectrum[k] mag / N / coh_gain * 2.0f; // 排除直流和奈奎斯特项时要小心 }参数说明coh_gain是窗函数的相干增益汉宁窗约为 0.5矩形窗为 1.0。除以coh_gain是为了补偿加窗造成的能量衰减。乘以 2.0 是单边谱换算因为正负频率分量对称只取一边时需要把幅值加倍。直流分量k0和奈奎斯特频率kN/2不翻倍否则幅值会虚高一倍。这些细节决定最终算出的幅值谱和真实幅度之间的误差能否控制在 1% 以内。4.2 常用窗函数的选择与参数对照FFT 点数固定时窗函数决定的是频谱泄漏和主瓣宽度的取舍。矩形窗主瓣最窄、频率分辨能力最好但旁瓣只衰减约 13 dB遇到近距双频信号会把弱信号掩盖。汉宁窗旁瓣衰减约 31 dB是音频和振动测试中最常用的窗。汉明窗在近旁瓣表现略好但远旁瓣衰减不如汉宁窗。平顶窗用于幅值精度要求高的场景代价是主瓣展宽到矩形窗的约 3.7 倍。窗类型主瓣宽度(归一化)旁瓣衰减幅值精度典型用途矩形1.0-13 dB差瞬态信号、整周期采样汉宁2.0-31 dB好音频、振动、一般频谱汉明2.0-41 dB(近)中语音处理平顶3.7-70 dB(远)极好幅值校准FFT 算法 C 语言实现本身不依赖任何窗加窗发生在变换前的时域乘法和变换后的幅值修正环节。如果只做频域峰值检测而不需要精确幅值可以省略除coh_gain这一步。4.3 点数不足时的补零操作与频率分辨率陷阱补零是指在尾部填充零值使序列长度达到 FFT 要求的 2 的幂。补零不增加真实频率分辨率频率分辨率由采样时长T N_orig / F_s决定补零只做插值让频谱曲线更平滑、峰值定位更精细。// 原始 300 点补零到 512 点 #define N_ORIG 300 #define N_FFT 512 float fft_in_re[N_FFT], fft_in_im[N_FFT]; for (int i 0; i N_FFT; i) { fft_in_re[i] (i N_ORIG) ? raw_data[i] : 0.0f; fft_in_im[i] 0.0f; }一个容易忽略的点是补零前应不应对原始数据加窗做法是先加窗再补零。因为补零等效于对加窗序列做更长区间的周期延拓如果先补零再加窗零值区域也参与加权会造成窗浪费有效数据段。反过来先加窗再补零窗只作用于实际数据段频谱形状才是预期的窗谱。4.4 实信号用实 FFTrfft还是复 FFT绝大多数工程输入是实信号。直接调用复 FFT 会把虚部全置零浪费一半运算量。优化方案是利用实序列频谱的共轭对称性把 N 点实序列打包成 N/2 点复序列做复 FFT再拆出 N/2 点正频率分量。这种实 FFTrfft在 N1024 时能节省约 40% 的计算时间。// 实 FFT 打包技巧示意: 将 even 放实部odd 放虚部 // 做 N/2 点复 FFT 后通过对称关系重组出 N 点实序列频谱 // X[k] 0.5 * (Z[k] conj(Z[N/2 - k])) // - 0.5j * (Z[k] - conj(Z[N/2 - k])) * e^(-j2πk/N)这个重组过程的推导依赖共轭对称性编码时容易在符号上出错。对于 N 不大的场景直接用复 FFT 省心得多N 大于 2048 且性能紧张时再考虑 rfft。5. 逆变换、帧重叠与工程落地的三个关键技巧5.1 IFFT 用 FFT 同函数实现的小技巧逆变换的计算公式与正变换只差一个共轭和 1/N 缩放。C 语言代码里可以复用正变换函数先把输入取共轭调用fft_c再对输出取共轭并除以 N。void ifft_c(float *re, float *im, unsigned int n) { // Step 1: 取共轭 for (unsigned int i 0; i n; i) im[i] -im[i]; // Step 2: 复用正变换 fft_c(re, im, n); // Step 3: 取共轭并缩放 float inv_n 1.0f / (float)n; for (unsigned int i 0; i n; i) { im[i] -im[i]; re[i] * inv_n; im[i] * inv_n; } }这种做法的好处是不需要维护两套旋向逻辑代码量小、易验证。代价是比专用 IFFT 实现多两次共轭操作但共轭仅改变符号位标志不是浮点运算开销几乎可忽略。5.2 流式处理中的分帧重叠与频谱更新策略FFT 算法 C 语言实现的典型应用是实时频谱显示器。直接每隔 N 个采样点做一次 N 点 FFT频谱会以帧为单位跳动且窗口边缘的不连续会导致频谱泄漏。重叠处理overlap能有效缓解。常见参数是重叠 50% 或 75%每帧取 N 个点但相邻帧起点仅移动 N/2 或 N/4 个采样点。// 帧移动步长: N/2 (50% 重叠) #define FRAME_N 512 #define HOP_LEN 256 static float ring_buf[FRAME_N * 2]; static int ring_pos 0; // 每收到 HOP_LEN 个新样本拼一帧做 FFT for (int i 0; i HOP_LEN; i) { // 将新样本推入环形缓冲buffer[frame_idx FRAME_N] 缓存下一帧起点数据 } // 加窗 → fft → 幅值谱更新显示使用重叠处理后帧率提高一倍50% 重叠频谱的时间平滑性更好但计算量也翻倍。选择重叠率要在 CPU 占用与视觉/听觉平滑度之间平衡。另外重叠帧之间可以做幅值平均或峰值保持前者适合观测稳态信号后者适合捕捉瞬态事件。这些后续处理已经不属于 FFT 本身但决定最终用户体验。5.3 验证 FFT 结果正确性的三组测试向量拿到任何人提供的 FFT 算法 C 代码第一件事是跑测试而不是直接集成。带三组测试向量进代码覆盖率基本足够第一组是直流信号输入x[n] 1.0N8。输出应为X[0]8其余全部为 0。如果虚部有微小非零值是浮点误差量级应在1e-6以下。第二组是整周期正弦波N1024采样率 1024 Hz输入x[n] sin(2π * 2 * n / 1024)。输出X[2]幅值应为 512X[1022]应为 512共轭对称其余接近 0。第三组是单位脉冲x[0]1其余为 0。输出幅值谱应全为 1相位谱全为 0。这一组用来测试旋转因子符号是否正确——符号反了相位谱会变成线性递增一眼就能看出。这三组向量跑完正确性基本可以确认。最后一步检查性能瓶颈如果点数超过 4096用-O2编译选项、开启硬浮点指令、确保旋转因子查表而不是实时计算这三条做到位C 语言 FFT 的性能已经接近理论峰值。本文还有配套的精品资源点击获取
