单片机FFT实战：从ADC采样到频谱显示的全流程嵌入式实现

简介本资源是一套面向嵌入式初学者与单片机开发者的FFT频谱分析实践方案聚焦于在资源受限的单片机平台上实现快速傅里叶变换并实时显示频谱结果。项目完整覆盖信号采样、复数运算优化、蝶形算法移植、LCD液晶驱动显示等关键环节特别适合电子类课程设计、智能仪器开发及音频/振动信号分析类实训场景。压缩包共15个文件含核心C源码fft.c、Keil工程配置文件uv2、opt、m51、编译中间产物lst、obj、hex及备份文件bak总大小仅64KB轻量紧凑便于快速导入与调试。目前已有265人下载学习资源结构清晰、模块分工明确——主程序负责ADC采样与FFT计算液晶显示模块直观呈现频域幅值分布配套汇编启动文件startup.a51和链接脚本lnp保障底层兼容性是理解嵌入式信号处理全流程的典型入门范例。1. 单片机上跑 FFT 不是“移植个库就完事”采样、定点、内存、实时性四重关卡必须逐个击破很多人看到“单片机 FFT”第一反应是网上一搜抄段 C 代码接个 ADC串口打印几个数字——结果频谱图毛刺飞溅、主频峰偏移、幅值严重失真甚至 FFT 输出全为零。根本原因在于FFT 在 PC 上是浮点密集计算在单片机里却是资源受限环境下的精密工程ADC 采样率没对齐频率轴就全错用 float 算 1024 点 FFTSTM32F4 的 RAM 直接爆掉定点缩放系数设错一位整个频域动态范围塌缩 6dB更别说中断嵌套导致采样丢点、DMA 传输未对齐引发数据错位。本文聚焦真实嵌入式场景——以 STM32F407 和 51 单片机为双主线讲清从 ADC 配置到频谱显示的完整链路不依赖 GUI 库、不调用现成 DSP 库除非明确说明、所有代码可直接在 Keil/MDK 或 SDCC 下编译运行。适合正在做音频监测、电机振动分析、电力谐波检测或毕业设计需要实测频谱的同学尤其适合那些已经烧过几块板子、发现“理论上能跑通实际上频谱完全不对”的开发者。2. 从 ADC 采样开始采样率、触发方式与数据对齐决定频谱精度上限FFT 结果的物理意义完全由采样参数定义。在单片机上这一步出错后续所有计算都是空中楼阁。常见误区是把“能采集到数据”等同于“采样合格”而实际需同时满足三要素等间隔、无混叠、相位连续。2.1 采样率设置不是越高越好而是要匹配目标频带与 FFT 点数假设你要分析 0–2 kHz 的信号如电机轴承故障特征频带根据奈奎斯特采样定理最低采样率需 ≥4 kHz。但实际中必须留余量若用 STM32F4 的 ADC1最高采样率约 2.4 MSPS单通道但受 GPIO 带宽和 DMA 吞吐限制稳定运行建议 ≤1 MSPS若用 STC8H 系列 51 单片机内部 ADC 最高约 500 kSPS推荐工作在 200 kSPS 以下以保证精度关键约束采样率 × 采集时间 FFT 点数 N。例如N1024目标分析时长 100 ms则采样率必须严格为 10240 Hz10.24 kHz。若设为 10 kHz采集 100 ms 得到 1000 点不足 1024强行补零会导致频谱泄漏加剧。提示不要用SystemCoreClock / prescaler粗略估算 ADC 时钟。务必查阅对应芯片手册中 ADCCLK 与采样周期Sampling Time的关系表。例如 STM32F407 的 ADC_SMPR1 寄存器中对 12-bit 模式采样时间需设为 480 个 ADCCLK 周期才能达到 12-bit 精度此时若 ADCCLK30 MHz单次采样耗时 ≈16 μs理论最大采样率 ≈62.5 kSPS。2.2 触发与同步用定时器 TRGO 触发 ADC杜绝软件延时抖动手写while(ADC_GetFlagStatus(ADC1, ADC_FLAG_EOC) RESET);这类轮询方式会导致采样间隔严重不均——尤其在中断频繁时。正确做法是配置通用定时器如 TIM2为向上计数模式ARR采样周期倒数如采样率 10 kHz → ARRSystemCoreClock/10000−1将 TIM2 的 TRGO 信号连接至 ADC1 的外部触发源EXTSEL0b101对应 TIM2_TRGO开启 ADC 连续转换模式CONT1与 DMA 请求DMA1启动 TIM2 后ADC 自动按精确周期启动转换DMA 将结果搬入指定缓冲区。// STM32F4 示例TIM2 触发 ADC1 DMA RCC_APB1PeriphClockCmd(RCC_APB1PERIPH_TIM2, ENABLE); RCC_APB2PeriphClockCmd(RCC_APB2PERIPH_ADC1 | RCC_APB2PERIPH_GPIOA, ENABLE); GPIO_InitTypeDef GPIO_InitStructure; GPIO_InitStructure.GPIO_Pin GPIO_Pin_0; // PA0 为 ADC1_IN0 GPIO_InitStructure.GPIO_Mode GPIO_Mode_AIN; GPIO_Init(GPIOA, GPIO_InitStructure); ADC_CommonInitTypeDef ADC_CommonInitStructure; ADC_CommonInitStructure.ADC_Prescaler ADC_Prescaler_Div4; ADC_CommonInit(ADC_CommonInitStructure); ADC_InitTypeDef ADC_InitStructure; ADC_InitStructure.ADC_Resolution ADC_Resolution_12b; ADC_InitStructure.ADC_ScanConvMode DISABLE; ADC_InitStructure.ADC_ContinuousConvMode ENABLE; ADC_InitStructure.ADC_ExternalTrigConv ADC_ExternalTrigConv_T2_TRGO; // 关键 ADC_InitStructure.ADC_DataAlign ADC_DataAlign_Right; ADC_InitStructure.ADC_NbrOfConversion 1; ADC_Init(ADC1, ADC_InitStructure); // DMA 配置省略初始化代码 DMA_InitTypeDef DMA_InitStructure; DMA_InitStructure.DMA_BufferSize 1024; // 与 FFT 点数一致 DMA_InitStructure.DMA_MemoryInc DMA_MemoryInc_Enable; DMA_InitStructure.DMA_PeripheralDataSize DMA_PeripheralDataSize_HalfWord; DMA_InitStructure.DMA_MemoryDataSize DMA_MemoryDataSize_HalfWord; DMA_InitStructure.DMA_DIR DMA_DIR_PeripheralToMemory; DMA_InitStructure.DMA_Mode DMA_Mode_Circular; // 循环模式持续采集 DMA_Init(DMA2_Stream0, DMA_InitStructure); ADC_DMACmd(ADC1, ENABLE); ADC_Cmd(ADC1, ENABLE); TIM_Cmd(TIM2, ENABLE);这段代码确保每 100 μs10 kHz产生一次 ADC 转换DMA 自动填满adc_buffer[1024]无 CPU 干预。注意DMA_Mode_Circular是关键避免缓冲区溢出后停止采集。2.3 数据预处理直流偏置消除与窗函数选择原始 ADC 数据常含显著直流分量如传感器零点漂移直接 FFT 会在 0 Hz 处出现巨大尖峰掩盖邻近低频成分。必须先减去均值uint32_t sum 0; for(uint16_t i 0; i N; i) sum adc_buffer[i]; uint32_t mean sum / N; for(uint16_t i 0; i N; i) { int16_t x (int16_t)adc_buffer[i] - (int16_t)mean; input_real[i] x; // 输入 FFT 的实部数组 }窗函数用于抑制频谱泄漏。单片机资源有限汉宁窗Hanning是性价比最优解计算仅需乘加系数可查表存储。其公式为$$ w(n) 0.5 - 0.5 \cos\left(\frac{2\pi n}{N-1}\right), \quad n0,1,\dots,N-1 $$生成 1024 点汉宁窗系数Q15 定点格式范围 −32768~32767const int16_t hanning_1024[1024] { 0, 6, 24, 54, 96, 150, 216, 294, /* ... 省略中间 1008 项 ... */, 294, 216, 150, 96, 54, 24, 6 }; // 实际使用时input_real[i] (int32_t)input_real[i] * hanning_1024[i] 15;注意窗函数会衰减信号有效幅度需在幅值计算后补偿如汉宁窗能量补偿因子为 1.5。若不做补偿相同信号在不同窗下幅值不可比。3. 定点 FFT 实现选型、缩放与内存布局决定能否在 64KB RAM 内跑通 1024 点在 STM32F4 上用arm_math.h的arm_cfft_radix4_q15()是最稳妥方案但需理解其底层约束在 51 单片机上则必须手写精简版基 2 DIT-FFT并严控中间变量生命周期。3.1 STM32F4ARM CMSIS-DSP 库的 Q15 模式深度解析CMSIS-DSP 提供arm_cfft_radix4_q15()要求输入为复数格式实部虚部交替且N 必须为 4 的整数幂128/256/512/1024。关键参数有三参数说明典型值影响*SFFT 实例结构体指针arm_cfft_sR_q15_len1024决定蝶形运算次数与 twiddle 因子表*p1输入/输出缓冲区复数格式fft_buffer[2048]长度 2×N单位为 int16_tifftFlag是否执行 IFFT0FFT设为 1 时执行逆变换bitReverseFlag是否位反转输出1输出按自然序排列否则需额外位反转初始化步骤arm_cfft_instance_q15 S; arm_cfft_init_q15(S, 1024); // 初始化 1024 点 FFT 实例 // 构造复数输入实部存偶地址虚部存奇地址 int16_t fft_buffer[2048]; // 1024 个复数 2048 个 int16_t for(uint16_t i 0; i 1024; i) { fft_buffer[2*i] input_real[i]; // 实部 fft_buffer[2*i1] 0; // 虚部置 0实信号 } arm_cfft_q15(S, fft_buffer, 0, 1); // 执行 FFT位反转输出缩放是核心难点Q15 格式数值范围为 −1.0 ~ 0.99997而 FFT 中间结果极易溢出。CMSIS 默认采用级联缩放stage-by-stage scaling每级蝶形后右移 1 位即除以 2最终结果需左移 log₂(N) 位恢复。对 1024 点需左移 10 位// 计算幅值谱|X[k]| sqrt(Re² Im²)Q15 输入 → Q30 中间 → Q15 输出 for(uint16_t k 0; k 1024; k) { int32_t re (int32_t)fft_buffer[2*k]; int32_t im (int32_t)fft_buffer[2*k1]; int32_t mag_sq re*re im*im; // Q30 × Q30 Q60需截断 int16_t mag (int16_t)(sqrtf((float)mag_sq / 65536.0f) * 32767.0f); // 归一化到 Q15 magnitude[k] mag; }提示sqrtf()在 FPU 开启时效率尚可但若追求极致速度可用查表法256 点平方根表 插值或改用arm_sqrt_q15()需将 mag_sq 右移 15 位转为 Q15。3.2 51 单片机手写基 2 DIT-FFT 的内存与循环优化策略STC8H 或 AT89C52 等 51 单片机 RAM 通常仅 256–1024 字节无法容纳 1024 点复数缓冲区。务实方案是64 点或 128 点 FFT并采用原位计算in-place与位反转索引预计算。核心优化点复数结构体压缩不用struct {int16_t r,i;}4 字节/点改用两个独立数组real[128]和imag[128]各 256 字节减少指针运算开销twiddle 因子查表预先计算 cos/sin 值存入 code 区ROM避免 runtime 浮点运算蝶形运算内联将Wn cos(2πn/N) j·sin(2πn/N)展开为实部/虚部计算消除复数乘法循环展开对每一级stage手动展开内层循环减少for判断开销。64 点 FFT 的蝶形层级数为 log₂64 6第stage级的跨度stride 2^stage旋转因子索引w_index k * (N / (2*stride))。关键代码片段SDCC 编译// 64 点基2 DIT-FFT输入 real[64], imag[64] void fft_64(int16_t real[], int16_t imag[]) { uint8_t stage, k, j, stride, w_index; int16_t tr, ti, ur, ui; // 位反转重排预计算好 bitrev[64] 数组 for(k 0; k 64; k) { j bitrev64[k]; if(j k) { tr real[k]; real[k] real[j]; real[j] tr; ti imag[k]; imag[k] imag[j]; imag[j] ti; } } // 6 级蝶形 for(stage 1; stage 6; stage) { stride 1 (stage-1); for(k 0; k 64; k 2*stride) { for(j 0; j stride; j) { w_index j * (64 stage); // twiddle index ur (int16_t)((long)cos_table[w_index] * real[kjstride] 15) - (int16_t)((long)sin_table[w_index] * imag[kjstride] 15); ui (int16_t)((long)sin_table[w_index] * real[kjstride] 15) (int16_t)((long)cos_table[w_index] * imag[kjstride] 15); tr real[kj] - ur; ti imag[kj] - ui; real[kj] ur; imag[kj] ui; real[kjstride] tr; imag[kjstride] ti; } } } }cos_table[]和sin_table[]为 Q15 格式预计算表共 32 项占用约 128 字节 ROM。此实现可在 STC8H3K64S264KB Flash/8KB RAM上稳定运行全程无动态内存分配。4. 频谱映射与显示从 FFT 输出到可读频率-幅值关系的硬核转换FFT 输出的索引k不是物理频率必须通过采样率fs和点数N映射为 Hz。更易被忽略的是幅值需归一化、对数压缩、频带合并否则串口打印的数字毫无工程意义。4.1 频率轴计算线性映射与奈奎斯特上限的硬约束FFT 输出X[k]对应的物理频率为$$ f_k k \times \frac{f_s}{N}, \quad k 0,1,\dots,N-1 $$但因输入为实信号X[k]关于N/2共轭对称有效频带仅为k0到kN/2含对应0到f_s/2奈奎斯特频率。例如fs10 kHz, N1024则k0→ 0 HzDCk1→ 9.766 Hzk512→ 5000 Hzfs/2注意k512是最后一个有效点k513至k1023是镜像应丢弃。4.2 幅值归一化补偿窗函数与 FFT 缩放的双重衰减原始 ADC 数据经汉宁窗后总能量衰减约 1/3理论值 0.375。CMSIS-DSP 的 Q15 FFT 在 1024 点下默认每级缩放最终幅值需乘以N/2 512才接近理论值。综合补偿因子为$$ \text{ScaleFactor} \frac{N}{2} \times \frac{1}{\text{WindowEnergy}} $$汉宁窗能量因子为sum(w²)/N 0.375故ScaleFactor ≈ 512 / 0.375 ≈ 1365。实际代码中为避免溢出常采用分步缩放// 计算 0~512 点的幅值谱k0 to 512 for(uint16_t k 0; k 512; k) { int32_t re (int32_t)fft_buffer[2*k]; int32_t im (int32_t)fft_buffer[2*k1]; int32_t mag_sq re*re im*im; // Q15 输入 → Q30 mag_sq → 右移 15 位得 Q15 幅值平方 int16_t mag_sq_q15 (int16_t)(mag_sq 15); // 开方得 Q15 幅值再乘以 1365Q15 格式为 136515? 不用 Q1.14 表示 1365/16384 int32_t scaled_mag ((int32_t)arm_sqrt_q15(mag_sq_q15) * 1365) 10; // 粗略缩放 magnitude_db[k] (scaled_mag 0) ? (int16_t)(20 * log10f((float)scaled_mag / 32767.0f)) : -120; }注意log10f()在无 FPU 的单片机上极慢。生产环境应改用查表线性插值或直接用arm_log10_q15()需将 scaled_mag 转为 Q15。4.3 分段频谱图按工程需求合并频带降低显示复杂度直接显示 513 个点的频谱对单片机 UI 不现实。典型做法是按倍频程octave或三分之一倍频程1/3 octave分组。例如电力谐波分析常用 50 Hz 基频的整数倍频带频带编号频率范围 (Hz)对应 k 范围fs10kHz, N1024合并方式10–50k0–5取最大值250–100k6–10取最大值3100–150k11–15取最大值…………1004950–5000k509–512取最大值代码实现#define BAND_NUM 100 int16_t band_max[BAND_NUM]; uint16_t k_start[BAND_NUM] {0,6,11,...}; // 预计算每带起始 k uint16_t k_end[BAND_NUM] {5,10,15,...}; // 预计算每带结束 k for(uint8_t b 0; b BAND_NUM; b) { int16_t max_val 0; for(uint16_t k k_start[b]; k k_end[b]; k) { if(magnitude_db[k] max_val) max_val magnitude_db[k]; } band_max[b] max_val; } // 此时 band_max[0] 对应 0–50 Hz 带的最大 dB 值可直接送 LCD 或串口这种分段大幅降低数据量100 点 vs 513 点且符合人耳听感和工业标准如 ISO 2631 振动评价。5. 实时性验证与抗干扰技巧用示波器看中断响应用已知信号标定系统误差写完代码只是起点真正落地需用硬件工具验证时序与精度。很多“频谱不准”的问题根源不在算法而在信号链路的隐性失真。5.1 中断响应时间测量确认 ADC 采样无丢点用示波器探头接 STM32 的某个 GPIO如 PC13在 ADC 转换完成中断ADC_IRQn服务函数开头置高结尾置低void ADC_IRQHandler(void) { GPIO_SetBits(GPIOC, GPIO_Pin_13); // 处理 DMA 传输完成标志非 EOC if(DMA_GetFlagStatus(DMA2_Stream0, DMA_FLAG_TCIF0) ! RESET) { DMA_ClearFlag(DMA2_Stream0, DMA_FLAG_TCIF0); // 触发 FFT 计算 fft_ready 1; } GPIO_ResetBits(GPIOC, GPIO_Pin_13); }观察 PC13 的脉冲宽度应稳定在 1–2 μsCortex-M4 主频 168 MHz 下。若脉冲宽度跳变剧烈说明中断被更高优先级任务阻塞需检查 SysTick 或其他外设中断配置。5.2 信号源标定用函数发生器注入正弦波验证频谱峰值位置这是最可靠的验证手段。步骤函数发生器输出 1 kHz 正弦波幅值调至 ADC 满量程的 70%如 2.5 Vpp → 1.25 V 峰值单片机以fs10.24 kHz采集 1024 点FFT 后查找magnitude_db[k]最大值位置k_max计算f_calc k_max × fs / N应与 1000 Hz 误差 0.5%若偏差大检查ADC 时钟是否准确用示波器测 PA0 波形周期、DMA 缓冲区是否被意外覆盖、FFT 输入数组是否被其他任务修改。5.3 抗干扰实战技巧电源滤波、PCB 布局与软件滤波协同电源噪声在 ADC 电源引脚VDDA就近放置 100 nF 10 μF 陶瓷电解电容避免开关电源纹波耦合PCB 布局模拟地AGND与数字地DGND单点连接于 ADC 旁走线远离高速数字线如 USB、SPI软件滤波对band_max[]数组做滑动平均如 5 帧消除突发干扰static int16_t band_hist[BAND_NUM][5]; // 每带保存最近 5 帧 for(uint8_t b 0; b BAND_NUM; b) { // 移动历史帧 for(uint8_t i 4; i 0; i--) band_hist[b][i] band_hist[b][i-1]; band_hist[b][0] band_max[b]; // 计算平均 int32_t sum 0; for(uint8_t i 0; i 5; i) sum band_hist[b][i]; band_smoothed[b] (int16_t)(sum / 5); }这些措施组合使用可使单片机频谱分析系统在工业现场稳定运行而非仅限实验室理想环境。本文还有配套的精品资源点击获取