简介：一套基于C语言的DSP FFT语音频谱分析示例工程，面向数字信号处理初学者及嵌入式开发者，旨在帮助理解快速傅里叶变换原理，并实践语音信号的频域分析方法。压缩包内含59个文件，体积约95KB，主要包含C源码、DSP工程配置文件（.pjt、.cmd、.lkf）、头文件、目标文件以及多段正弦测试音频（MP3格式），覆盖1kHz至16kHz多个频率点，便于验证频谱计算结果的正确性。资源目前已吸引268人学习，适合在CCS等DSP开发环境中结合语音处理项目参考。通过阅读源码与工程结构，读者可掌握FFT库调用、采样数据预处理及频谱幅度计算的基本流程，并借助附带的多频测试音频对比输入与输出，从而快速上手语音频谱分析相关开发。
1. 为什么做语音频谱的人最终都会回到 FFT 和 C 语言
在嵌入式设备上做语音频谱显示，最不容易绕开的链路是：麦克风拾音、ADC 采样、帧切分、频谱变换和可视化。时域波形很难直接告诉你一段人声里有没有 800Hz 的共振峰，而 FFT 能在几十毫秒内把这帧信号的能量按频率拆开。对用 C 语言做 DSP 开发的工程师来说，光会调用一个 FFT 函数不够，还得清楚点数怎么选、窗函数加不加、输出数组里每个下标对应多少赫兹，这些参数直接决定频谱图是可用还是只是几条乱跳的柱子。下文按算法原理、DSP 平台实现、语音频谱完整流程、验证与调优四段展开，适合正在用 STM32、常见 MCU 或定点 DSP 做音频分析的开发者，也适合准备从控制类程序转音频处理的 C 语言开发。
2. 从 DFT 到 FFT：C 语言实现频谱计算的最小闭环
2.1 为什么不直接套 DFT 公式
离散傅里叶变换的教科书定义是 X[k] = Σ x[n]·e^(-j2πkn/N)，每计算一个频点需要 N 次复数乘法，N 个频点就是 N² 次。当 N 取 1024，DFT 要做一百多万次复数乘法，在大多数没有专门复数乘法部件的处理器上要消耗毫秒级时间；语音一帧通常只有 20~30ms，留给频谱变换的预算本就不多，后面还要接取模、显示、降噪等步骤，DFT 几乎没有余量。
FFT 走的是另一条路：利用旋转因子 W_N^(kn) 的周期性和对称性，把长序列逐级拆成短序列，每一级只做 N/2 次蝶形运算，总计算量变成 (N/2)·log2 N。同样是 1024 点，FFT 只需要 5120 次蝶形运算，比 DFT 少两个数量级。后面要实现的基 2 FFT 就是这条思路最典型的落地结构，也是大多数 C 语言和 DSP 库实现的基础。两类实现的乘法量差别很直观：
| FFT 点数 | DFT 复数乘法次数 | FFT 蝶形次数 | 加速比（约） |
|---|---|---|---|
| 256 | 65536 | 1024 | 64 |
| 1024 | 1048576 | 5120 | 205 |
2.2 基 2 FFT 的 C 语言实现：位逆序加蝶形运算
常见的基 2 时间抽取 FFT 分两步：先把输入序列按二进制位逆序重排，再做 log2 N 级蝶形迭代。下面这段 C 语言代码可以在 PC 上直接编译验证，也方便后续移植到嵌入式环境。
#include <math.h> #include <stdint.h> #define FFT_N 256 #define PI 3.14159265358979323846f typedef struct { float real; float imag; } complex_t; // 位逆序重排：把 x[i] 换到 bit-reverse(i) 对应的位置 static void bit_reverse(complex_t *x, uint16_t n) { for (uint16_t i = 1, j = 0; i < n; i++) { uint16_t bit = n >> 1; for (; j & bit; bit >>= 1) { j ^= bit; } j ^= bit; if (i < j) { complex_t tmp = x[i]; x[i] = x[j]; x[j] = tmp; } } } // 基 2 时间抽取 FFT，点数 n 必须是 2 的幂 void fft_radix2(complex_t *x, uint16_t n) { bit_reverse(x, n); for (uint16_t len = 2; len <= n; len <<= 1) { float ang = -2.0f * PI / (float)len; complex_t w = {cosf(ang), sinf(ang)}; // 本级旋转因子 for (uint16_t i = 0; i < n; i += len) { complex_t wn = {1.0f, 0.0f}; for (uint16_t j = 0; j < len / 2; j++) { complex_t u = x[i + j]; complex_t v = { wn.real * x[i + j + len / 2].real - wn.imag * x[i + j + len / 2].imag, wn.real * x[i + j + len / 2].imag + wn.imag * x[i + j + len / 2].real }; x[i + j].real = u.real + v.real; x[i + j].imag = u.imag + v.imag; x[i + j + len / 2].real = u.real - v.real; x[i + j + len / 2].imag = u.imag - v.imag; float wreal = wn.real * w.real - wn.imag * w.imag; float wimag = wn.real * w.imag + wn.imag * w.real; wn.real = wreal; wn.imag = wimag; } } } }代码逻辑分三层循环：len表示当前蝶形的跨度，从 N=2 开始逐级翻倍，直到覆盖整个序列；内层i枚举每一组蝶形的起点；最内层j完成一组里上下两个点的复数乘加，也就是经典的“蝶形”操作。wn从 1 开始，每处理完一个点就乘以 w，相当于用迭代代替重复计算三角函数，代价是乘出来的累积误差，对大多数 8kHz 语音应用，float 精度足够；不建议改成 double，否则在无 FPU 平台上成倍增加开销。
需要注意两点：bit_reverse里的嵌套循环是位运算实现下标翻转的标准写法，这对 i 从 1 开始循环，j 通过按位与和异或模拟翻位过程，运行前要求 n 为 2 的幂，否则位掩码的移位逻辑会出错，这是基 2 FFT 的硬性前提。调用时把实数序列装进complex_t数组，虚部置 0，执行fft_radix2(x, FFT_N)后，x[k].real和x[k].imag就是第 k 个频点的实部与虚部。
2.3 从 FFT 输出数组换算成频率和幅值
变换完的数组长度还是 N，但只有前半部分是有效频谱。以 fs=8000Hz、N=256 为例，频率分辨率 Δf = fs / N = 31.25Hz，第 k 个频点对应的频率是 k × 31.25Hz；由于采样定理的限制，实际可用的频率范围是 0 到 fs/2=4000Hz，也就是 k 从 0 到 128 共 129 个有效频点，k 大于 128 的部分与前半部分镜像对称，在语音频谱显示中通常直接丢弃。
幅值计算最常用的是幅度谱：mag[k] = sqrt(re² + im²)。但在显示端几乎没有人直接画线性幅度，因为语音信号低频能量往往比高频高几个数量级，线性绘制会把高频细节压成一条贴着零轴的线。常规做法是先取 20·log10(mag[k]) 转成 dB 单位，再映射到屏幕高度，这样人耳感知的音量级与显示高度才近似线性。
3. DSP 平台上的 FFT 实现：定点、查表与资源权衡
3.1 浮点和定点：DSP 上先决定用什么数域
大多数人先在 PC 上用 float 写通 FFT，然后直接烧到 MCU 或 DSP 上跑，结果发现耗时翻了几倍甚至几十倍。原因无非两种：处理器没有硬件浮点单元，或浮点需要软件模拟。于是 DSP 平台上的常见做法是切换到 Q15 定点表示：一个 16 位整数，最高位为符号位，数值范围落在 [-1, 1)，1 对应 0x7FFF，-1 对应 0x8000。如果 ADC 输出是偏移二进制格式，先减去 0x8000 再进入 FFT；如果是二进制补码格式，则可以直接按 Q15 使用，省去一次格式转换。
定点 FFT 与浮点的差异集中在乘法上。两个 Q15 数相乘，结果为 Q30 格式，也就是 30 个小数位，必须算术右移 15 位还原成 Q15，否则数据范围对不上；而加减运算结果可能超出 [-1, 1)，这时要么依赖 DSP 的饱和指令，要么手动限幅。在浮点代码中一个简单的中间变量，到定点实现里就变成“先移位、再饱和”三步操作，这也是很多工程师从 C 语言浮点转 DSP 定点后第一周的主要调试内容。
3.2 用旋转因子查表法替代运行时三角函数
浮点实现里的cosf和sinf在定点 DSP 上是非常奢侈的调用，库函数可能展开成多项式逼近，一次调用几百个周期。更常见的做法是把旋转因子提前算好，固化在 const 数组中，运行时用查表替代三角函数计算。基 2 FFT 里旋转因子 W = e^(-j2πk/len)，最完整的一张表需要存 N 个点；实际利用对称性只需存 N/4 组并映射四个象限，可以显著减少 Flash 占用。
// twiddle_table.h 由上位机脚本生成，包含 cos_table 与 sin_table // 数组元素为 Q15 定点值，长度 FFT_N / 4 #include "twiddle_table.h" // 用位与取模，等价于 angle % FFT_N，要求 FFT_N 为 2 的幂 static uint16_t angle_to_index(uint16_t angle) { return angle & (FFT_N - 1); } // 查表取旋转因子：angle 为 0 到 N-1 的索引，wr/wi 为输出 Q15 值 static void twiddle_load(uint16_t angle, int16_t *wr, int16_t *wi) { uint16_t idx = angle_to_index(angle); if (idx < FFT_N / 4) { *wr = cos_table[idx]; *wi = sin_table[idx]; } else if (idx < FFT_N / 2) { uint16_t k = FFT_N / 2 - idx; *wr = -cos_table[k]; *wi = sin_table[k]; } else if (idx < 3 * FFT_N / 4) { uint16_t k = idx - FFT_N / 2; *wr = -cos_table[k]; *wi = -sin_table[k]; } else { uint16_t k = FFT_N - idx; *wr = cos_table[k]; *wi = -sin_table[k]; } }这段代码的关键在象限折叠：angle取 0 到 N-1 之间的索引，通过位与操作把角度映射到第一象限对应的查表下标，再根据所在象限对余弦和正弦取正负。这样只需要 N/4 组表项就能覆盖完整旋转因子。表本身由 PC 端脚本生成后复制进 C 文件，脚本用 Python 或 MATLAB 都行，生成的是 Q15 整型常量，DSP 上运行时只做查表和符号翻转，不再碰三角函数。
在真正的定点蝶形运算里，查表得到的wr、wi要和复数乘法的结果一样走 Q15 乘法流程：乘积右移 15 位，检测溢出再做饱和。这里最常见的错误是对中途结果移位不足，导致频谱输出近似正确但某些频点出现明显跳变；排查方法也比较直接，用一组全 0 输入跑一次，再输入只在 0 号位置为 32767 的冲激信号，输出应该严格等于表值，任何偏差都指向移位或溢出问题。
3.3 不同 FFT 点数用存储和实时性换什么
语音频谱分析中 FFT 点数是需要拍板的第一参数。点数太少，频率分辨率粗，相邻频带区分不开；点数太多，帧长变长，显示刷新率和实时性下降。下面给出 8kHz 采样率下的典型配置参考。
| FFT 点数 | 帧长 ms | 频率分辨率 Hz | 有效频谱范围 Hz | 复数数组占用（8 字节/复数） | 蝶形级数 |
|---|---|---|---|---|---|
| 128 | 16 | 62.5 | 0~4000 | 1KB | 7 |
| 256 | 32 | 31.25 | 0~4000 | 2KB | 8 |
| 512 | 64 | 15.625 | 0~4000 | 4KB | 9 |
| 1024 | 128 | 7.8125 | 0~4000 | 8KB | 10 |
这里的帧长是 N/fs。语音清音段大约几十毫秒，浊音段周期约几毫秒到十几毫秒，256 点帧长 32ms 正好覆盖一个完整音调周期，所以多数语音频谱显示用 256 点起步。如果做音乐类、需要区分相邻半个音（约 3%）的场景，则要 512 或 1024 点。还要注意表格里只算了原始复数数组，实际加上窗函数、旋转因子表和输出功率谱数组，Flash 和 RAM 占用大概是表的 2~3 倍，选片时留出余量。
4. 语音频谱分析的完整流程：从 ADC 采样到频谱显示
4.1 先定帧参数：采样率、帧长与重叠
语音频谱分析不是把麦克风数据源源不断丢进 FFT 就能出结果，而是要先把连续语音切成“帧”。帧长由采样率和 FFT 点数共同决定：采样率 8kHz、FFT 点数 256 时，一帧时长为 32ms，与分析点要求一致。相邻帧之间通常还有一定重叠，比如帧移 16ms 即 50% 重叠，这样频谱随时间变化能更平滑，减少帧边界造成的跳变。
如果采样率提高，比如 16kHz，而 FFT 点数还是 256，那频谱范围变成 0~8kHz，分辨率反而降到 62.5Hz。换句话讲，FFT 点数决定的是分辨率，采样率决定的是最高可分析频率，两者不能混为一谈。常见做法是先用最低满足需求的采样率，再按最小可接受的分辨率选点数，避免无谓增加计算负担。
4.2 加窗为什么不可省：矩形窗的频谱泄漏
直接截取一帧相当于给无限长语音信号乘了一个矩形窗，频谱上表现为主瓣很窄但旁瓣很高，频率成分会“泄漏”到相邻频点。比如一个正好不落在整数频点上的正弦信号，矩形窗会让它在好几个 bin 都有能量，峰值位置看起来模糊不清。语音频谱分析通常使用 Hamming 窗或 Hann 窗，它们旁瓣衰减更大，代价是主瓣略宽、分辨率略有损失。
Hamming 窗系数 w[n] = 0.54 - 0.46·cos(2πn/(N-1))。这组系数可以先算好放在数组里，加上窗时逐点与输入信号相乘，每帧都要执行一次。因为 FFT 本身要求周期延拓，而语音是非周期信号，不加窗的频谱会出现明显的拖尾，加窗之后频谱更“干净”，峰值更好定位。
4.3 从 ADC 数据到幅度谱的 C 代码实现
下面这段代码把采样、加窗、FFT 和幅度谱串在一起。它假设 FFT 函数已经是前面实现的fft_radix2，或者后续替换为 DSP 厂商库的接口。
#include <stdint.h> #include <math.h> #define PI 3.14159265358979323846f #define SAMPLES 256 // 每帧采样点数 #define FS 8000 // 采样率 Hz #define FRAME_SHIFT 128 // 50% 重叠，帧移 = 16ms // Hamming 窗系数，常驻 RAM 或 flash static float hamming_window[SAMPLES]; static void build_hamming(void) { for (int n = 0; n < SAMPLES; n++) { hamming_window[n] = 0.54f - 0.46f * cosf(2.0f * PI * n / (SAMPLES - 1)); } } // 处理一帧：输入原始 ADC 序列，输出有效频点的幅度谱（单位 dBFS） // fft_points 必须小于等于 SAMPLES，且为 2 的幂 void process_frame(const int16_t *adc_buf, float *mag_db, int fft_points) { static complex_t spec[SAMPLES]; // 专用数组，避免每次进入重新分配 for (int i = 0; i < fft_points; i++) { spec[i].real = (float)adc_buf[i] * hamming_window[i]; spec[i].imag = 0.0f; } fft_radix2(spec, fft_points); for (int k = 0; k <= fft_points / 2; k++) { float re = spec[k].real; float im = spec[k].imag; float mag = sqrtf(re * re + im * im); // 转换成 dBFS，参照满幅 32768 做归一化 float norm = mag / (32768.0f * fft_points / 2.0f); mag_db[k] = 20.0f * log10f(norm + 1e-9f); } }参数说明：adc_buf是 DMA 或中断采集来的原始 16 位数据，虚部直接清零。process_frame里第一个循环完成加窗，第二个循环调用 FFT，第三个循环取有效频点幅度并转成 dBFS。1e-9f是防止 log10 收到 0 而出现负无穷，这个技巧在对数频谱里几乎必须用。归一化系数32768 * N / 2是因为输入满幅正弦的 FFT 峰值约为幅值×N/2，除以它才能把 0dBFS 对齐到满幅信号，否则相同的输入在不同 N 下显示电平不一样。Hamming 窗会把信号能量压缩到约 0.54 倍，如果要求更严格的 0dBFS 对齐，可以在归一化系数里再除以窗系数的平均幅度，实际显示中影响不大。
调用时需要考虑帧移。缓冲区用环形队列保存最新 256 个点，每攒够FRAME_SHIFT个新样本就取出一帧调用process_frame，频谱数组刷新率就是 8000/128=62.5Hz，大约每 16ms 更新一次显示，视觉上足够平滑。如果某些显示模块刷新慢，可以提高帧移到一帧长度，代价是频谱对语音变化的跟随性变差。
4.4 结果怎么读：频点、峰值与语音特征
帧处理完后的mag_db数组，下标 k 对应的频率是 k×31.25Hz。假如语音里有明显的 250Hz 基频，那么峰值会出现在 k≈8 的位置；在 250Hz 整数倍的位置还会看到谐波峰。找到最大峰值的常用方法是遍历 k=1 到 128（跳过 k=0 的直流分量），记录最大幅度及其下标。谱峰位置对应音高，谱峰高度对比对应音强，而高频区能量占比可以用来判断清音与浊音的相对强弱。
5. 用已知信号验证 FFT 代码，再做三处进阶处理
5.1 合成 1kHz 正弦的验证方法
代码移植到新平台后，第一件事不是接麦克风，而是用合成信号验证链路。生成一段 1kHz 正弦，采样率 8kHz，FFT 点数 256，理想情况下频率分辨率 31.25Hz，真实频率 1000/31.25=32，也就是第 32 个频点应该出现一个明显的峰值，其他频点接近零。用下面这段代码可以直接在 PC 上定位峰值下标。
#include <stdio.h> #include <math.h> int main(void) { complex_t x[FFT_N]; for (int n = 0; n < FFT_N; n++) { x[n].real = 32767.0f * sinf(2.0f * PI * 1000.0f * n / 8000.0f); x[n].imag = 0.0f; } fft_radix2(x, FFT_N); int peak_idx = 0; float peak_val = 0.0f; for (int k = 1; k <= FFT_N / 2; k++) { float mag = sqrtf(x[k].real * x[k].real + x[k].imag * x[k].imag); if (mag > peak_val) { peak_val = mag; peak_idx = k; } } printf("peak idx = %d, freq = %.2f Hz\n", peak_idx, peak_idx * 8000.0f / FFT_N); return 0; }输出应该是 32 附近。如果差一到两个 bin，检查是不是忘了加窗或者采样率宏定义错了；如果差得离谱，多半是位逆序函数或蝶形因子符号写反，可以输入仅 0 号非零的冲激信号，FFT 输出应该全部是同一个实数值。
5.2 频率分辨率不够时：补零和插值
语音频谱里两个相邻频率分量相差不到 31.25Hz 时，直接用 256 点 FFT 区分不开。简单办法是把 256 点数据后面补 256 个零做 512 点 FFT，计算量变大，但频谱曲线更平滑，峰值定位精度有所提高；要注意补零并不能让两个真实紧邻的频率真正分离，它只是让频谱“插值变细”。想要真正提高分辨率就必须增加实际采样点数，或改用频谱插值算法在峰值附近用抛物线拟合出更精确的频率估计，后者在音高检测场景很常用。
5.3 用峰值保持让语音频谱显示更稳定
直接逐帧刷新幅度谱时，人眼会看到柱子抖动很厉害，因为语音短时波动完全映射到了显示上。常见做法是给每个频点做时间方向的平滑或峰值保持。峰值保持的做法是记录每一帧各频点的最大值，并按一个随时间衰减的遗忘因子缓慢降低，只有当新帧的幅度超过当前显示值时才更新为更高的值。这样既保留瞬态冲击的峰值痕迹，又不会让显示每个帧都剧烈跳动。频率轴和幅度轴映射到屏幕时，频率轴建议用对数间隔显示，与人的听觉特性更接近。所有参数调通后再接入实时麦克风数据，用手机放一段人声或单音吉他，频谱中和基频与泛音相关的峰就会清晰可见。
本文还有配套的精品资源，点击获取
