1. 项目概述在STM32上玩转DSP库的FFT最近在做一个基于STM32F103的简易频谱分析仪项目核心需求是从ADC采集一段音频信号然后快速分析出里面的主要频率成分。这种需求在振动监测、音频处理、电力谐波分析等领域太常见了。一开始我琢磨着自己手搓一个FFT快速傅里叶变换算法但试了下发现在72MHz主频的Cortex-M3内核上用纯C语言实现一个256点的FFT计算耗时感人根本达不到实时性要求。于是我把目光投向了ST官方提供的DSP库。这个库里面封装了高度优化的汇编级FFT函数号称性能强劲。网上能找到的例程不少但大多只是把代码跑起来屏幕上显示个频谱图就完事了。对于“为什么采样率要这么设”、“FFT结果数组到底怎么解读”、“幅值计算里的65536又是哪来的”这些关键问题往往一笔带过。我在实际调试时就在这些地方栽了跟头频谱显示总是怪怪的。所以我决定结合ST官方那个经典的FFT演示工程就是那个FFT_Demo把整个从信号生成、FFT调用到结果处理的流程彻底拆解清楚。这篇文章不光会带你走通代码更重要的是我会把我调试过程中遇到的坑、对关键参数的理解以及如何将这套方法移植到你自己的实际项目比如接上真实的ADC中的经验毫无保留地分享出来。无论你是刚接触数字信号处理的嵌入式新手还是想优化现有算法性能的老鸟相信都能从中找到有用的东西。2. 核心思路与DSP库解析2.1 为什么选择STM32的DSP库在资源受限的MCU上做FFT我们无非几个选择一是自己用C写可读性好但速度慢二是用第三方优化库如CMSIS-DSP通用性强三就是直接用芯片厂商提供的库。ST的DSP库属于第三种它最大的优势是针对自家Cortex-M3内核指令集做了深度汇编优化。这个库里的FFT函数比如我们例程里用的cr4_fft_64_stm32其内部是用汇编编写的充分利用了STM32的硬件乘法和单周期指令速度比纯C实现快出一个数量级。对于追求实时性的应用这点性能提升至关重要。库函数处理的数据格式是Q15定点数范围-32768到32767这正好匹配STM32的ADC输出通常是12位0-4095也避免了浮点运算STM32F103没有FPU带来的巨大开销。注意这个DSP库是较老的版本现在ST主推的是集成在CubeMX和HAL库里的CMSIS-DSP库。但老库的代码结构更清晰对于理解底层原理非常有帮助且其核心的优化思想是相通的。学会这个再迁移到新的CMSIS-DSP会非常轻松。2.2 官方例程框架拆解提供的代码是一个完整的工程我们可以把它看成一次“离线仿真”。它没有连接真实的ADC而是通过软件生成一个包含特定频率的正弦波信号然后对这个信号做FFT最后把频谱结果显示到LCD上。这个流程完美模拟了真实场景信号源 - 采样 - FFT变换 - 频谱显示。整个工程的核心函数调用链非常清晰main()-MyDualSweep()-MygSin()生成信号 -cr4_fft_xx_stm32()执行FFT -powerMag()计算幅值。其中MyDualSweep函数做了两轮测试第一轮生成单个频率变化的信号第二轮生成两个频率叠加的信号。这其实就是在验证FFT的频率分辨能力和对多频率信号的分离能力是判断FFT算法是否工作正常的经典方法。3. 关键参数深度解读与配置3.1 采样定理与参数定义看懂这个例程第一步必须吃透代码开头那几个宏定义。它们不是随便写的数字每一个都对应着数字信号处理的核心理论。#define FreqSample 3200 // 采样频率 #define NPT 64 // 采样点数 #define ScalingFactor 20000 // 幅值系数采样频率 (FreqSample) 这里设为3200 Hz。根据奈奎斯特采样定理能无失真还原的最高信号频率 (Fmax) 必须小于采样频率 (Fs) 的一半即Fmax Fs/2。这里的Fs/2 1600 Hz意味着这个系统能分析的最高频率是1600Hz。如果你要分析50Hz的工频信号这个采样率绰绰有余但如果你想分析音频20Hz-20kHz那就得大幅提高FreqSample比如设为44.1kHz。采样点数 (NPT) 这里是64点。FFT点数直接影响两个关键指标频率分辨率 (Δf) 即频谱图上相邻两根谱线代表的频率差。计算公式是Δf Fs / NPT。本例中Δf 3200 Hz / 64 50 Hz。这意味着频谱图横轴上第0个点代表0Hz直流分量第1个点代表50Hz第2个点代表100Hz以此类推。如果你想区分两个靠得很近的频率比如51Hz和52Hz64点的分辨率50Hz是做不到的你必须增加NPT比如1024点来减小Δf。计算量和存储开销 FFT点数越多计算时间越长需要的RAM也越多输入、输出数组都要变大。STM32F103的RAM很小需要权衡。幅值系数 (ScalingFactor) 这个20000是用来调整生成的正弦波幅值的。因为FFT库函数使用Q15格式-32767到32767所以生成的信号幅值不能超过这个范围。ScalingFactor设为20000再乘以一个范围在[-1, 1]左右的正弦值可以确保信号幅值在合理范围内避免运算溢出。3.2 FFT函数调用与数据格式例程中FFT的调用非常简单根据NPT的值选择不同的函数#if (NPT64) cr4_fft_64_stm32(lBUFOUT, lBUFIN, NPT); #elif (NPT256) cr4_fft_256_stm32(lBUFOUT, lBUFIN, NPT); #endif这里有两个关键数组lBUFIN[NPT]:输入数组。注意它的数据类型是long32位。但FFT处理的是复数STM32 DSP库用一种巧妙的方式打包每个long型数据的高16位bits 31:16存放实部Re低16位bits 15:0存放虚部Im。对于我们从ADC采到的实信号虚部全部填0即可。这就是为什么在MygSin函数里你会看到lBUFIN[i] ((short)fZ) 16 ;这行代码它把生成的实部数据左移16位低16位自然就是0虚部。lBUFOUT[NPT]:输出数组。格式和输入数组完全一样也是每个元素的高16位是实部低16位是虚部。FFT的结果是复数我们需要从这个数组中提取出每个频率分量的幅度和相位信息。实操心得很多初学者在这里会困惑为什么输入输出用同一个数组不行因为大多数FFT算法包括这个是“原位运算”理论上输入输出可以共用内存。但ST的这个库函数设计成了输入输出分开可能是为了接口清晰或兼容性。在实际应用时如果内存紧张可以仔细研究函数说明看是否支持原位计算。4. 从复数结果到实际频谱幅值计算详解4.1powerMag函数逐行解析FFT输出是复数但我们通常只关心幅度谱。powerMag函数就是干这个的。这是整个工程里最容易出错也最需要理解的一环。void powerMag(long nfill, char* strPara) { static int32_t lX,lY; uint32_t i; for (i0; i nfill/2; i) { // 注意只循环前NPT/2个点 lX (lBUFOUT[i] 16)16; /* 取低16bit实部这里注释有误 */ lY (lBUFOUT[i] 16); /* 取高16bit虚部这里注释有误 */ { float X NPT * ((float)lX) / 32768; float Y NPT * ((float)lY) / 32768; float Mag sqrt(X*X Y*Y) / nfill; lBUFMAG[i] (uint32_t)(Mag * 65536); } } }首先纠正一个关键误区原代码注释把lX和lY说反了。根据我们之前讲的数据格式高16位实部低16位虚部(lBUFOUT[i] 16)16 先左移16位丢掉高16位再右移16位符号扩展回来。这个操作取的是低16位即虚部(Im)。(lBUFOUT[i] 16) 直接右移16位取的是高16位即实部(Re)。 所以lX是虚部(Im)lY是实部(Re)。这点一定要记清楚否则后面计算全错。然后我们看计算过程float X NPT * ((float)lX) / 32768; 将取出的Q15格式的虚部值范围-32768~32767转换为浮点数并除以32768归一化到[-1, 1)区间再乘以点数NPT。为什么乘以NPT这是因为FFT的常用定义中变换后的幅度会放大NPT倍。除以NPT是后续步骤。float Y ... 对实部做同样的处理。float Mag sqrt(X*X Y*Y) / nfill; 计算复数的模幅度。sqrt(Re^2 Im^2)得到的是该频率分量的“能量”大小。最后除以nfill即NPT才是正确的幅值这是很多初学者忽略的一步如果不除得到的幅值会非常大。lBUFMAG[i] (uint32_t)(Mag * 65536); 将计算出的浮点幅值乘以65536再转成整型。这其实是将结果转换为了Q16格式即小数点在bit16之前方便后续的定点显示或传输。如果你用LCD直接显示浮点数这一步可以不用。为什么循环只到nfill/2因为对于实信号输入其FFT结果具有共轭对称性。即第k个点和第N-k个点是共轭的幅度相同相位相反。所以只有前NPT/2个点从0到NPT/2-1包含独立的频率信息后一半是前一半的镜像没有实际意义。第0个点是直流分量频率0Hz第NPT/2个点是奈奎斯特频率分量本例中1600Hz。4.2 单边谱与双边谱代码中有一个未使用的onesided函数和传入的2SIDED参数这涉及到单边谱和双边谱的概念。双边谱 显示了从-Fs/2到Fs/2的整个频谱共轭对称。我们上面计算出来的lBUFMAG数组长度NPT/2对应的是从0Hz到Fs/2的部分这已经是单边谱了。单边谱 为了能量守恒在绘制单边谱时除了直流分量0点和奈奎斯特频率点NPT/2点其他点的幅值需要乘以2。因为真实信号的能量被平均分到了正负频率上。onesided函数试图做这个事但例程里没调用它。对于大多数实际应用如显示频谱我们使用计算出的lBUFMAG即单边谱的前半部分就足够了。如果你想得到物理上正确的幅值对于i从1到NPT/2-1的点需要将lBUFMAG[i]乘以2。5. 移植到真实ADC采集项目官方例程用函数生成信号是理想情况。真正有挑战的是把ADC采集的实时数据灌进去。下面是我的移植步骤和踩坑记录。5.1 ADC配置与数据缓冲假设我们用STM32F103的ADC1规则通道以固定的采样率Fs采集信号。定时器触发ADC 这是保证等间隔采样的关键。配置一个定时器如TIM2使其溢出频率等于你想要的采样频率Fs。然后在定时器更新事件中触发ADC转换。// 示例系统时钟72MHz定时器预分频72则计数器时钟为1MHz。 // 若要Fs3200Hz则自动重载值ARR 1MHz / 3200Hz 312.5取整312或313会略有误差。 TIM_TimeBaseInitTypeDef TIM_InitStructure; TIM_InitStructure.TIM_Period 312; // 自动重载值 TIM_InitStructure.TIM_Prescaler 71; // 预分频值72-1 // ... 其他配置 TIM_SelectOutputTrigger(TIM2, TIM_TRGOSource_Update); // 设置更新事件为触发输出 ADC_ExternalTrigConvCmd(ADC1, ENABLE); ADC_ExternalTrigConvConfig(ADC1, ADC_ExternalTrigConv_T2_TRGO); // ADC由TIM2触发DMA双缓冲循环采集 这是实现连续无缝采集的标配。配置DMA将ADC的转换结果寄存器ADC1-DR自动搬运到内存中的一个缓冲区。设置双缓冲区Ping-Pong Buffer每个缓冲区大小为NPTFFT点数。当一个缓冲区填满时触发DMA半传输/传输完成中断在中断里启动FFT计算同时DMA继续向另一个缓冲区写入数据。#define NPT 256 uint16_t adc_buffer[2][NPT]; // 双缓冲区 volatile uint8_t buffer_ready 0; // 缓冲区就绪标志0:无1:缓冲区0满2:缓冲区1满 // DMA中断服务函数 void DMA1_Channel1_IRQHandler(void) { if(DMA_GetITStatus(DMA1_IT_TC1)) { // 传输完成第二个缓冲区满 DMA_ClearITPendingBit(DMA1_IT_TC1); buffer_ready 2; // 标记缓冲区1数据就绪 } else if(DMA_GetITStatus(DMA1_IT_HT1)) { // 半传输完成第一个缓冲区满 DMA_ClearITPendingBit(DMA1_IT_HT1); buffer_ready 1; // 标记缓冲区0数据就绪 } }5.2 数据预处理与FFT调用在主循环或一个专门的任务中检查buffer_ready标志然后处理数据。数据搬移与格式转换 ADC的结果是12位无符号整数0-4095。需要先转换成有符号的Q15格式并打包成FFT库要求的格式。void ProcessFFTBuffer(uint8_t buf_id) { uint16_t *src_buf adc_buffer[buf_id]; for(int i0; iNPT; i) { // 1. 将12位ADC值0-4095转换为有符号16位-32768~32767 // 先减去直流偏置假设是中间值2048再缩放。 int32_t temp (int32_t)src_buf[i] - 2048; // 去直流 temp temp * 8; // 粗略缩放使动态范围接近Q15满量程。可根据实际信号调整。 // 确保不溢出 if(temp 32767) temp 32767; if(temp -32768) temp -32768; // 2. 打包高16位为实部(有符号值)低16位为0虚部 lBUFIN[i] ((int16_t)temp) 16; } // 3. 执行FFT #if (NPT256) cr4_fft_256_stm32(lBUFOUT, lBUFIN, NPT); #endif // 4. 计算幅值 powerMag(NPT, 2SIDED); // 5. 后续显示频谱、查找峰值频率等 DisplaySpectrum(lBUFMAG, NPT/2); }去直流与加窗 上面的代码简单减去了一个固定值2048来去直流。但实际中信号的直流分量可能会漂移。更稳健的做法是计算一个缓冲区数据的平均值然后每个点都减去这个平均值。此外如果采样不是整周期会发生“频谱泄漏”主峰会扩散到旁边的频率点。为了抑制泄漏需要对采样数据加窗如汉宁窗Hamming、汉明窗Hanning。加窗函数需要预先计算好一个长度为NPT的数组在数据打包到lBUFIN之前先乘以窗函数系数。// 预先计算汉宁窗系数Q15格式 int16_t hanning_window[NPT]; for(int i0; iNPT; i) { hanning_window[i] (int16_t)(32767 * (0.5 - 0.5 * cos(2*PI*i/(NPT-1)))); } // 在数据转换时加窗 temp ((int32_t)src_buf[i] - dc_offset) * 8; // dc_offset是计算出的直流平均值 temp (temp * hanning_window[i]) 15; // Q15乘法右移15位得到结果 lBUFIN[i] ((int16_t)temp) 16;5.3 频谱显示与峰值查找计算出lBUFMAG数组后就可以显示频谱了。在LCD上通常用柱状图来表示。幅值缩放lBUFMAG是Q16格式的数值可能很大。需要将其缩放到LCD屏幕的高度范围内。#define LCD_HEIGHT 240 uint32_t max_mag 0; // 首先找到本次频谱中的最大值 for(int i1; iNPT/2; i) { // 通常忽略直流分量(i0) if(lBUFMAG[i] max_mag) max_mag lBUFMAG[i]; } // 绘制频谱柱状图 for(int i1; iNPT/2; i) { int bar_height (int)((float)lBUFMAG[i] / max_mag * (LCD_HEIGHT - 20)); if(bar_height LCD_HEIGHT-20) bar_height LCD_HEIGHT-20; LCD_DrawLine(i*2, LCD_HEIGHT, i*2, LCD_HEIGHT - bar_height, RED); // 简单绘图 }峰值频率查找 除了显示我们常常需要知道主频率是多少。最简单的就是找lBUFMAG数组中忽略前几个可能包含直流泄漏的点最大值对应的索引i。uint32_t peak_index 0; uint32_t peak_value 0; for(int i5; iNPT/2; i) { // 从第5个点开始找避开直流和低频干扰 if(lBUFMAG[i] peak_value) { peak_value lBUFMAG[i]; peak_index i; } } float peak_freq (float)peak_index * FreqSample / NPT; // 计算实际频率 printf(Peak Frequency: %.2f Hz\n, peak_freq);更高级的算法可以插值如重心法来获得比Δf更精确的频率估计。6. 调试心得与常见问题排查6.1 频谱结果不对一步步排查问题频谱全是噪声没有明显峰值。检查ADC数据 首先别急着做FFT。把lBUFIN数组里的实部数据高16位通过串口打印出来或者用IDE的Memory/Graph功能画出来看。确保你采集到的是一个正常的、幅度合适的正弦波或你想要的信号。如果这里信号就不对问题在ADC和采样环节。检查数据格式 确认你打包到lBUFIN的格式是否正确。高16位是实部有符号Q15低16位是0。可以用一个简单测试手动给lBUFIN填充一个已知频率的正弦波序列像例程的MygSin那样如果这样FFT结果正确但用ADC数据就不行那肯定是数据转换或ADC环节的问题。问题峰值频率和预期对不上。检查采样频率(Fs)和点数(NPT) 重新核算Δf Fs / NPT。你预期的频率f应该大致等于peak_index * Δf。如果偏差是固定的倍数可能是Fs计算错了定时器配置错误。如果偏差是随机的可能是信号本身频率不稳或者存在频谱泄漏导致峰值扩散。验证FFT函数 用例程的MygSin函数生成一个50Hz的信号Fs3200, NPT64则Δf50Hz理论上峰值应该出现在lBUFMAG[1]。用这个标准测试来验证你的powerMag计算过程是否正确。问题频谱幅值看起来不合理太大或太小。检查powerMag函数 重点看sqrt(X*X Y*Y) / nfill这一步除以nfill即NPT了吗这是最常见的错误。如果不除幅值会大NPT倍。检查缩放系数 从ADC原始值到Q15的缩放系数上面代码中的* 8需要根据你信号的实际情况调整。先用一个已知幅度的信号比如用信号发生器输入一个1Vpp的正弦波进行校准调整这个系数使得FFT计算出的幅值与实际物理幅值成比例。问题存在很大的直流分量lBUFMAG[0]值特别大。去直流处理 这说明你的信号有很强的直流偏置。务必在FFT前减去整个缓冲区的平均值dc_offset而不是一个固定值。硬件耦合 检查硬件电路如果是交流信号是否使用了隔直电容6.2 性能优化与进阶技巧整数开方优化powerMag函数中的sqrt浮点运算在M3上比较慢。如果不需要非常精确的幅值或者只需要找峰值可以用绝对值求和法近似Mag (abs(X) abs(Y)) / nfill。或者使用平方和查表法预先计算好X*XY*Y与Mag的对应关系。使用实数FFTRFFT ST的DSP库和CMSIS-DSP库都提供了专门的实数FFT函数如arm_rfft_fast_f32。对于实信号输入RFFT的计算量比复数FFT大约少一半而且输出格式更友好前半部分就是复数FFT结果的前N/21个点。这是未来升级到CMSIS-DSP后的首选。动态调整采样率 如果你要分析的频率范围很宽可以采用多采样率策略。对于低频高精度分析用高点数、低采样率对于高频分析用低点数、高采样率。这需要在不同定时器配置间切换。使用DSP库的定点数函数 如果切换到CMSIS-DSP它提供了丰富的定点数Q31, Q15FFT函数如arm_cfft_q15。这些函数同样高度优化并且输入输出格式更标准独立的实部数组和虚部数组避免了手动打包的麻烦。移植这个FFT例程到实际项目最磨人的不是代码本身而是对信号处理概念的理解和细节的把握。我最开始就是卡在数据格式和幅值计算上出来的频谱怎么看怎么怪。后来是把每一步的中间数据都打印出来和MATLAB计算的结果一点点对比才最终定位到问题。嵌入式DSP调试一定要有“分而治之”的思路先确保信号源和ADC是对的再验证FFT函数本身最后处理幅值转换和显示。当你第一次在小小的LCD屏幕上看到清晰的、随着输入信号频率变化而跳动的频谱柱时那种成就感绝对是熬夜调试的最佳回报。
STM32 DSP库FFT实战：从原理到ADC频谱分析移植
发布时间：2026/6/7 15:19:55
1. 项目概述在STM32上玩转DSP库的FFT最近在做一个基于STM32F103的简易频谱分析仪项目核心需求是从ADC采集一段音频信号然后快速分析出里面的主要频率成分。这种需求在振动监测、音频处理、电力谐波分析等领域太常见了。一开始我琢磨着自己手搓一个FFT快速傅里叶变换算法但试了下发现在72MHz主频的Cortex-M3内核上用纯C语言实现一个256点的FFT计算耗时感人根本达不到实时性要求。于是我把目光投向了ST官方提供的DSP库。这个库里面封装了高度优化的汇编级FFT函数号称性能强劲。网上能找到的例程不少但大多只是把代码跑起来屏幕上显示个频谱图就完事了。对于“为什么采样率要这么设”、“FFT结果数组到底怎么解读”、“幅值计算里的65536又是哪来的”这些关键问题往往一笔带过。我在实际调试时就在这些地方栽了跟头频谱显示总是怪怪的。所以我决定结合ST官方那个经典的FFT演示工程就是那个FFT_Demo把整个从信号生成、FFT调用到结果处理的流程彻底拆解清楚。这篇文章不光会带你走通代码更重要的是我会把我调试过程中遇到的坑、对关键参数的理解以及如何将这套方法移植到你自己的实际项目比如接上真实的ADC中的经验毫无保留地分享出来。无论你是刚接触数字信号处理的嵌入式新手还是想优化现有算法性能的老鸟相信都能从中找到有用的东西。2. 核心思路与DSP库解析2.1 为什么选择STM32的DSP库在资源受限的MCU上做FFT我们无非几个选择一是自己用C写可读性好但速度慢二是用第三方优化库如CMSIS-DSP通用性强三就是直接用芯片厂商提供的库。ST的DSP库属于第三种它最大的优势是针对自家Cortex-M3内核指令集做了深度汇编优化。这个库里的FFT函数比如我们例程里用的cr4_fft_64_stm32其内部是用汇编编写的充分利用了STM32的硬件乘法和单周期指令速度比纯C实现快出一个数量级。对于追求实时性的应用这点性能提升至关重要。库函数处理的数据格式是Q15定点数范围-32768到32767这正好匹配STM32的ADC输出通常是12位0-4095也避免了浮点运算STM32F103没有FPU带来的巨大开销。注意这个DSP库是较老的版本现在ST主推的是集成在CubeMX和HAL库里的CMSIS-DSP库。但老库的代码结构更清晰对于理解底层原理非常有帮助且其核心的优化思想是相通的。学会这个再迁移到新的CMSIS-DSP会非常轻松。2.2 官方例程框架拆解提供的代码是一个完整的工程我们可以把它看成一次“离线仿真”。它没有连接真实的ADC而是通过软件生成一个包含特定频率的正弦波信号然后对这个信号做FFT最后把频谱结果显示到LCD上。这个流程完美模拟了真实场景信号源 - 采样 - FFT变换 - 频谱显示。整个工程的核心函数调用链非常清晰main()-MyDualSweep()-MygSin()生成信号 -cr4_fft_xx_stm32()执行FFT -powerMag()计算幅值。其中MyDualSweep函数做了两轮测试第一轮生成单个频率变化的信号第二轮生成两个频率叠加的信号。这其实就是在验证FFT的频率分辨能力和对多频率信号的分离能力是判断FFT算法是否工作正常的经典方法。3. 关键参数深度解读与配置3.1 采样定理与参数定义看懂这个例程第一步必须吃透代码开头那几个宏定义。它们不是随便写的数字每一个都对应着数字信号处理的核心理论。#define FreqSample 3200 // 采样频率 #define NPT 64 // 采样点数 #define ScalingFactor 20000 // 幅值系数采样频率 (FreqSample) 这里设为3200 Hz。根据奈奎斯特采样定理能无失真还原的最高信号频率 (Fmax) 必须小于采样频率 (Fs) 的一半即Fmax Fs/2。这里的Fs/2 1600 Hz意味着这个系统能分析的最高频率是1600Hz。如果你要分析50Hz的工频信号这个采样率绰绰有余但如果你想分析音频20Hz-20kHz那就得大幅提高FreqSample比如设为44.1kHz。采样点数 (NPT) 这里是64点。FFT点数直接影响两个关键指标频率分辨率 (Δf) 即频谱图上相邻两根谱线代表的频率差。计算公式是Δf Fs / NPT。本例中Δf 3200 Hz / 64 50 Hz。这意味着频谱图横轴上第0个点代表0Hz直流分量第1个点代表50Hz第2个点代表100Hz以此类推。如果你想区分两个靠得很近的频率比如51Hz和52Hz64点的分辨率50Hz是做不到的你必须增加NPT比如1024点来减小Δf。计算量和存储开销 FFT点数越多计算时间越长需要的RAM也越多输入、输出数组都要变大。STM32F103的RAM很小需要权衡。幅值系数 (ScalingFactor) 这个20000是用来调整生成的正弦波幅值的。因为FFT库函数使用Q15格式-32767到32767所以生成的信号幅值不能超过这个范围。ScalingFactor设为20000再乘以一个范围在[-1, 1]左右的正弦值可以确保信号幅值在合理范围内避免运算溢出。3.2 FFT函数调用与数据格式例程中FFT的调用非常简单根据NPT的值选择不同的函数#if (NPT64) cr4_fft_64_stm32(lBUFOUT, lBUFIN, NPT); #elif (NPT256) cr4_fft_256_stm32(lBUFOUT, lBUFIN, NPT); #endif这里有两个关键数组lBUFIN[NPT]:输入数组。注意它的数据类型是long32位。但FFT处理的是复数STM32 DSP库用一种巧妙的方式打包每个long型数据的高16位bits 31:16存放实部Re低16位bits 15:0存放虚部Im。对于我们从ADC采到的实信号虚部全部填0即可。这就是为什么在MygSin函数里你会看到lBUFIN[i] ((short)fZ) 16 ;这行代码它把生成的实部数据左移16位低16位自然就是0虚部。lBUFOUT[NPT]:输出数组。格式和输入数组完全一样也是每个元素的高16位是实部低16位是虚部。FFT的结果是复数我们需要从这个数组中提取出每个频率分量的幅度和相位信息。实操心得很多初学者在这里会困惑为什么输入输出用同一个数组不行因为大多数FFT算法包括这个是“原位运算”理论上输入输出可以共用内存。但ST的这个库函数设计成了输入输出分开可能是为了接口清晰或兼容性。在实际应用时如果内存紧张可以仔细研究函数说明看是否支持原位计算。4. 从复数结果到实际频谱幅值计算详解4.1powerMag函数逐行解析FFT输出是复数但我们通常只关心幅度谱。powerMag函数就是干这个的。这是整个工程里最容易出错也最需要理解的一环。void powerMag(long nfill, char* strPara) { static int32_t lX,lY; uint32_t i; for (i0; i nfill/2; i) { // 注意只循环前NPT/2个点 lX (lBUFOUT[i] 16)16; /* 取低16bit实部这里注释有误 */ lY (lBUFOUT[i] 16); /* 取高16bit虚部这里注释有误 */ { float X NPT * ((float)lX) / 32768; float Y NPT * ((float)lY) / 32768; float Mag sqrt(X*X Y*Y) / nfill; lBUFMAG[i] (uint32_t)(Mag * 65536); } } }首先纠正一个关键误区原代码注释把lX和lY说反了。根据我们之前讲的数据格式高16位实部低16位虚部(lBUFOUT[i] 16)16 先左移16位丢掉高16位再右移16位符号扩展回来。这个操作取的是低16位即虚部(Im)。(lBUFOUT[i] 16) 直接右移16位取的是高16位即实部(Re)。 所以lX是虚部(Im)lY是实部(Re)。这点一定要记清楚否则后面计算全错。然后我们看计算过程float X NPT * ((float)lX) / 32768; 将取出的Q15格式的虚部值范围-32768~32767转换为浮点数并除以32768归一化到[-1, 1)区间再乘以点数NPT。为什么乘以NPT这是因为FFT的常用定义中变换后的幅度会放大NPT倍。除以NPT是后续步骤。float Y ... 对实部做同样的处理。float Mag sqrt(X*X Y*Y) / nfill; 计算复数的模幅度。sqrt(Re^2 Im^2)得到的是该频率分量的“能量”大小。最后除以nfill即NPT才是正确的幅值这是很多初学者忽略的一步如果不除得到的幅值会非常大。lBUFMAG[i] (uint32_t)(Mag * 65536); 将计算出的浮点幅值乘以65536再转成整型。这其实是将结果转换为了Q16格式即小数点在bit16之前方便后续的定点显示或传输。如果你用LCD直接显示浮点数这一步可以不用。为什么循环只到nfill/2因为对于实信号输入其FFT结果具有共轭对称性。即第k个点和第N-k个点是共轭的幅度相同相位相反。所以只有前NPT/2个点从0到NPT/2-1包含独立的频率信息后一半是前一半的镜像没有实际意义。第0个点是直流分量频率0Hz第NPT/2个点是奈奎斯特频率分量本例中1600Hz。4.2 单边谱与双边谱代码中有一个未使用的onesided函数和传入的2SIDED参数这涉及到单边谱和双边谱的概念。双边谱 显示了从-Fs/2到Fs/2的整个频谱共轭对称。我们上面计算出来的lBUFMAG数组长度NPT/2对应的是从0Hz到Fs/2的部分这已经是单边谱了。单边谱 为了能量守恒在绘制单边谱时除了直流分量0点和奈奎斯特频率点NPT/2点其他点的幅值需要乘以2。因为真实信号的能量被平均分到了正负频率上。onesided函数试图做这个事但例程里没调用它。对于大多数实际应用如显示频谱我们使用计算出的lBUFMAG即单边谱的前半部分就足够了。如果你想得到物理上正确的幅值对于i从1到NPT/2-1的点需要将lBUFMAG[i]乘以2。5. 移植到真实ADC采集项目官方例程用函数生成信号是理想情况。真正有挑战的是把ADC采集的实时数据灌进去。下面是我的移植步骤和踩坑记录。5.1 ADC配置与数据缓冲假设我们用STM32F103的ADC1规则通道以固定的采样率Fs采集信号。定时器触发ADC 这是保证等间隔采样的关键。配置一个定时器如TIM2使其溢出频率等于你想要的采样频率Fs。然后在定时器更新事件中触发ADC转换。// 示例系统时钟72MHz定时器预分频72则计数器时钟为1MHz。 // 若要Fs3200Hz则自动重载值ARR 1MHz / 3200Hz 312.5取整312或313会略有误差。 TIM_TimeBaseInitTypeDef TIM_InitStructure; TIM_InitStructure.TIM_Period 312; // 自动重载值 TIM_InitStructure.TIM_Prescaler 71; // 预分频值72-1 // ... 其他配置 TIM_SelectOutputTrigger(TIM2, TIM_TRGOSource_Update); // 设置更新事件为触发输出 ADC_ExternalTrigConvCmd(ADC1, ENABLE); ADC_ExternalTrigConvConfig(ADC1, ADC_ExternalTrigConv_T2_TRGO); // ADC由TIM2触发DMA双缓冲循环采集 这是实现连续无缝采集的标配。配置DMA将ADC的转换结果寄存器ADC1-DR自动搬运到内存中的一个缓冲区。设置双缓冲区Ping-Pong Buffer每个缓冲区大小为NPTFFT点数。当一个缓冲区填满时触发DMA半传输/传输完成中断在中断里启动FFT计算同时DMA继续向另一个缓冲区写入数据。#define NPT 256 uint16_t adc_buffer[2][NPT]; // 双缓冲区 volatile uint8_t buffer_ready 0; // 缓冲区就绪标志0:无1:缓冲区0满2:缓冲区1满 // DMA中断服务函数 void DMA1_Channel1_IRQHandler(void) { if(DMA_GetITStatus(DMA1_IT_TC1)) { // 传输完成第二个缓冲区满 DMA_ClearITPendingBit(DMA1_IT_TC1); buffer_ready 2; // 标记缓冲区1数据就绪 } else if(DMA_GetITStatus(DMA1_IT_HT1)) { // 半传输完成第一个缓冲区满 DMA_ClearITPendingBit(DMA1_IT_HT1); buffer_ready 1; // 标记缓冲区0数据就绪 } }5.2 数据预处理与FFT调用在主循环或一个专门的任务中检查buffer_ready标志然后处理数据。数据搬移与格式转换 ADC的结果是12位无符号整数0-4095。需要先转换成有符号的Q15格式并打包成FFT库要求的格式。void ProcessFFTBuffer(uint8_t buf_id) { uint16_t *src_buf adc_buffer[buf_id]; for(int i0; iNPT; i) { // 1. 将12位ADC值0-4095转换为有符号16位-32768~32767 // 先减去直流偏置假设是中间值2048再缩放。 int32_t temp (int32_t)src_buf[i] - 2048; // 去直流 temp temp * 8; // 粗略缩放使动态范围接近Q15满量程。可根据实际信号调整。 // 确保不溢出 if(temp 32767) temp 32767; if(temp -32768) temp -32768; // 2. 打包高16位为实部(有符号值)低16位为0虚部 lBUFIN[i] ((int16_t)temp) 16; } // 3. 执行FFT #if (NPT256) cr4_fft_256_stm32(lBUFOUT, lBUFIN, NPT); #endif // 4. 计算幅值 powerMag(NPT, 2SIDED); // 5. 后续显示频谱、查找峰值频率等 DisplaySpectrum(lBUFMAG, NPT/2); }去直流与加窗 上面的代码简单减去了一个固定值2048来去直流。但实际中信号的直流分量可能会漂移。更稳健的做法是计算一个缓冲区数据的平均值然后每个点都减去这个平均值。此外如果采样不是整周期会发生“频谱泄漏”主峰会扩散到旁边的频率点。为了抑制泄漏需要对采样数据加窗如汉宁窗Hamming、汉明窗Hanning。加窗函数需要预先计算好一个长度为NPT的数组在数据打包到lBUFIN之前先乘以窗函数系数。// 预先计算汉宁窗系数Q15格式 int16_t hanning_window[NPT]; for(int i0; iNPT; i) { hanning_window[i] (int16_t)(32767 * (0.5 - 0.5 * cos(2*PI*i/(NPT-1)))); } // 在数据转换时加窗 temp ((int32_t)src_buf[i] - dc_offset) * 8; // dc_offset是计算出的直流平均值 temp (temp * hanning_window[i]) 15; // Q15乘法右移15位得到结果 lBUFIN[i] ((int16_t)temp) 16;5.3 频谱显示与峰值查找计算出lBUFMAG数组后就可以显示频谱了。在LCD上通常用柱状图来表示。幅值缩放lBUFMAG是Q16格式的数值可能很大。需要将其缩放到LCD屏幕的高度范围内。#define LCD_HEIGHT 240 uint32_t max_mag 0; // 首先找到本次频谱中的最大值 for(int i1; iNPT/2; i) { // 通常忽略直流分量(i0) if(lBUFMAG[i] max_mag) max_mag lBUFMAG[i]; } // 绘制频谱柱状图 for(int i1; iNPT/2; i) { int bar_height (int)((float)lBUFMAG[i] / max_mag * (LCD_HEIGHT - 20)); if(bar_height LCD_HEIGHT-20) bar_height LCD_HEIGHT-20; LCD_DrawLine(i*2, LCD_HEIGHT, i*2, LCD_HEIGHT - bar_height, RED); // 简单绘图 }峰值频率查找 除了显示我们常常需要知道主频率是多少。最简单的就是找lBUFMAG数组中忽略前几个可能包含直流泄漏的点最大值对应的索引i。uint32_t peak_index 0; uint32_t peak_value 0; for(int i5; iNPT/2; i) { // 从第5个点开始找避开直流和低频干扰 if(lBUFMAG[i] peak_value) { peak_value lBUFMAG[i]; peak_index i; } } float peak_freq (float)peak_index * FreqSample / NPT; // 计算实际频率 printf(Peak Frequency: %.2f Hz\n, peak_freq);更高级的算法可以插值如重心法来获得比Δf更精确的频率估计。6. 调试心得与常见问题排查6.1 频谱结果不对一步步排查问题频谱全是噪声没有明显峰值。检查ADC数据 首先别急着做FFT。把lBUFIN数组里的实部数据高16位通过串口打印出来或者用IDE的Memory/Graph功能画出来看。确保你采集到的是一个正常的、幅度合适的正弦波或你想要的信号。如果这里信号就不对问题在ADC和采样环节。检查数据格式 确认你打包到lBUFIN的格式是否正确。高16位是实部有符号Q15低16位是0。可以用一个简单测试手动给lBUFIN填充一个已知频率的正弦波序列像例程的MygSin那样如果这样FFT结果正确但用ADC数据就不行那肯定是数据转换或ADC环节的问题。问题峰值频率和预期对不上。检查采样频率(Fs)和点数(NPT) 重新核算Δf Fs / NPT。你预期的频率f应该大致等于peak_index * Δf。如果偏差是固定的倍数可能是Fs计算错了定时器配置错误。如果偏差是随机的可能是信号本身频率不稳或者存在频谱泄漏导致峰值扩散。验证FFT函数 用例程的MygSin函数生成一个50Hz的信号Fs3200, NPT64则Δf50Hz理论上峰值应该出现在lBUFMAG[1]。用这个标准测试来验证你的powerMag计算过程是否正确。问题频谱幅值看起来不合理太大或太小。检查powerMag函数 重点看sqrt(X*X Y*Y) / nfill这一步除以nfill即NPT了吗这是最常见的错误。如果不除幅值会大NPT倍。检查缩放系数 从ADC原始值到Q15的缩放系数上面代码中的* 8需要根据你信号的实际情况调整。先用一个已知幅度的信号比如用信号发生器输入一个1Vpp的正弦波进行校准调整这个系数使得FFT计算出的幅值与实际物理幅值成比例。问题存在很大的直流分量lBUFMAG[0]值特别大。去直流处理 这说明你的信号有很强的直流偏置。务必在FFT前减去整个缓冲区的平均值dc_offset而不是一个固定值。硬件耦合 检查硬件电路如果是交流信号是否使用了隔直电容6.2 性能优化与进阶技巧整数开方优化powerMag函数中的sqrt浮点运算在M3上比较慢。如果不需要非常精确的幅值或者只需要找峰值可以用绝对值求和法近似Mag (abs(X) abs(Y)) / nfill。或者使用平方和查表法预先计算好X*XY*Y与Mag的对应关系。使用实数FFTRFFT ST的DSP库和CMSIS-DSP库都提供了专门的实数FFT函数如arm_rfft_fast_f32。对于实信号输入RFFT的计算量比复数FFT大约少一半而且输出格式更友好前半部分就是复数FFT结果的前N/21个点。这是未来升级到CMSIS-DSP后的首选。动态调整采样率 如果你要分析的频率范围很宽可以采用多采样率策略。对于低频高精度分析用高点数、低采样率对于高频分析用低点数、高采样率。这需要在不同定时器配置间切换。使用DSP库的定点数函数 如果切换到CMSIS-DSP它提供了丰富的定点数Q31, Q15FFT函数如arm_cfft_q15。这些函数同样高度优化并且输入输出格式更标准独立的实部数组和虚部数组避免了手动打包的麻烦。移植这个FFT例程到实际项目最磨人的不是代码本身而是对信号处理概念的理解和细节的把握。我最开始就是卡在数据格式和幅值计算上出来的频谱怎么看怎么怪。后来是把每一步的中间数据都打印出来和MATLAB计算的结果一点点对比才最终定位到问题。嵌入式DSP调试一定要有“分而治之”的思路先确保信号源和ADC是对的再验证FFT函数本身最后处理幅值转换和显示。当你第一次在小小的LCD屏幕上看到清晰的、随着输入信号频率变化而跳动的频谱柱时那种成就感绝对是熬夜调试的最佳回报。
