STM32定点FFT精度实战：采样链路与幅值频率相位测量

简介面向嵌入式开发者的STM32信号处理实践资源聚焦在STM32F407平台借助标准库实现FFT把时域正弦波采样转换为频域分析结果从而精准测量幅值、频率与相位差。内容覆盖ADC采样、窗函数抑制频谱泄露、FFT变换、频谱分析与结果输出等完整链路适合已具备C语言与STM32外设基础、希望把快速傅里叶变换落到工业控制、电机控制、音频分析或无线通信等实时场景的中级学习者。压缩包共283个文件约9.98MB包含46个c与60个h源文件以及Keil工程文件、库文件、hex与map等分别对应算法实现、外设驱动、编译配置与最终产物便于直接编译验证与二次移植。已有1668人学习下载可据此复现从采集到频域计算的完整流程理解FFT结果与真实信号参数的对应关系并借鉴其工程组织方式排查采样、窗函数与频谱换算中的常见偏差。1. 1024 点 FFT 跑通的那一刻精度问题才刚刚开始很多人第一次在 STM32 上跑通 FFT看到串口吐出一串频域数据、峰值点正好落在某个 bin 上就觉得这事已经成了。可把信号源调到准确的 1 kHz读数却跳出 1003.9 Hz把输入幅值从 1.0 V 压到 0.5 V算出来的幅值没有线性减半两路同相的正弦波相位差却稳定偏了十几度。快速傅里叶变换本身没算错错的是从模拟前端到频谱读数之间那几段没人讲透的链路采样时刻、窗函数增益、栅栏效应、双通道同步。搞懂这套机制的价值在于你能分清哪些误差是算法固有的、哪些是自己接线的锅而不是对着一个 0.3% 的频率偏差反复重写蝶形运算。2. 从 DFT 到蝶形O(N²) 是怎么被旋转因子对称性砍成 O(N log N) 的离散傅里叶变换的定义本身很短N 点变换要做 N² 次复数乘加。1024 点就是一百多万次复数乘法在 72 MHz 的 Cortex-M 上按软件浮点跑单次测量要几十毫秒根本谈不上实时。FFT 不是新算法它只是把 DFT 里大量重复的乘法提前算好、复用起来。理解这一层才能明白为什么定点实现里旋转因子的位宽、查表大小会直接影响最终精度。2.1 旋转因子的周期性与共轭对称DFT 的核心因子是 W_N^k e^(-j2πk/N)。它有两个性质周期性 W_N^(kN) W_N^k共轭对称性 W_N^(kN/2) -W_N^k。这意味着 N 点变换里真正独立的因子只有前 N/2 个后一半只是前一半取负。基 2 分解正是靠这条性质把一次 N 点变换拆成两个 N/2 点变换再把两个 N/2 点变换的结果用一对加减组合起来组合时只需要乘一次旋转因子。递归下去乘法次数从 N² 降到 (N/2)·log₂N。落到实现上旋转因子通常预生成成一张表长度 N/2元素是 Q15 格式的 cos 和 sin。表越长精度越高但占的 Flash 也越多1024 点、两个 q15_t 数组一共 4 KB对多数 STM32 型号是可以接受的。常见的实现思路是把表放在 Flash 里用 const 修饰运行时只读不写省掉 RAM。2.2 基 2 蝶形与位反转重排时域抽取的基 2 FFT 有一个前置动作把输入序列按二进制位反转的顺序重排。原因是每一级分解都把偶数下标和奇数下标分到两组分解到底之后原始下标的顺序被打乱了重排回去才能保证输出是自然顺序。位反转本身不涉及乘法但它是原地运算的前提跳过它直接算出来的频谱是乱序的峰值位置对不上任何物理频率。蝶形运算是整个算法的原子操作取一对输入 (u, v)做 uv·W 和 u−v·W两个输出写回原来的两个位置。每一级的所有蝶形互相独立可以原地进行这是 FFT 能省内存的关键。#define FFT_N 1024 /* 预生成旋转因子表cos(-2*pi*k/N)、sin(-2*pi*k/N)Q15 定标 */ extern const q15_t tw_cos[FFT_N / 2]; extern const q15_t tw_sin[FFT_N / 2]; /* 基 2 时域抽取 FFT原地运算进入前必须完成位反转重排 */ void fft_radix2_q15(q15_t *xr, q15_t *xi) { uint16_t len, i, j, k; for (len 2; len FFT_N; len 1) { /* 逐级扩大蝶形跨度 */ uint16_t half len 1; uint16_t step FFT_N / len; /* 旋转因子表索引步进 */ for (i 0; i FFT_N; i len) { /* 遍历本级所有蝶形块 */ k 0; for (j i; j i half; j) { /* 块内每一对 */ q15_t ur xr[j], ui xi[j]; q15_t vr xr[j half], vi xi[j half]; /* 复数乘法(vr j*vi) * (cos j*sin)32 位承接防溢出 */ int32_t tr (int32_t)vr * tw_cos[k] - (int32_t)vi * tw_sin[k]; int32_t ti (int32_t)vr * tw_sin[k] (int32_t)vi * tw_cos[k]; q15_t wr (q15_t)(tr 15); /* Q30 右移回 Q15 */ q15_t wi (q15_t)(ti 15); xr[j] (q15_t)(ur wr); xi[j] (q15_t)(ui wi); xr[j half] (q15_t)(ur - wr); xi[j half] (q15_t)(ui - wi); k step; /* 跳到下一个旋转因子 */ } } } }这段代码有几处值得展开。step FFT_N / len是旋转因子索引的步进跨度越大块内需要的因子越稀疏索引每次跳step正好对应 W_N^(j·N/len)。中间乘积用int32_t承接是因为两个 Q15 相乘会得到 Q30直接存回 q15_t 会截断符号位。 15完成 Q30 到 Q15 的定标回退但要注意这是算术右移负数才能保持正确。蝶形最后的加减法是最容易翻车的地方|ur| 和 |wr| 各自最大接近 32767相加轻松超过 int16 范围所以输入信号必须做功率归一化或者整体右移一位否则高电平输入下频谱会出鬼影。2.3 定点 Q15 下的动态范围取舍Q15 表示 [-1, 1) 区间的定点小数分辨率 1/32768。ADC 是 12 位原始码值 0~4095转成 Q15 时要先减去直流偏置 2048再左移四位相当于乘 16这样满量程 2048 映射到 32768 附近。注意左移四位后最大值 2048×16 32768 已经溢出 int16 正半轴实际实现里通常只左移三位或者减偏置后乘 15。代价是动态范围变窄。一个 1 mV 的小信号叠在 1 V 的直流上减完偏置后可能只剩几十个码值量化噪声占比很高。这是定点 FFT 的固有问题不是代码写错。如果测量对象是小信号要么在前端加交流耦合和放大要么换成浮点 FFT——后者在带 FPU 的型号上单次 1024 点大约 1~2 ms比定点慢但仍可用。3. 采样链路为什么定时器触发 ADC DMA 是标准库下的稳妥组合FFT 的结果只对“等间隔采样”这个前提负责。一旦采样时刻抖动频谱就会糊掉再精确的蝶形运算也救不回来。标准库时代没有 HAL 的便利封装定时器 TRGO 触发 ADC、DMA 搬运、半满中断处理这套组合几乎是唯一能把抖动压到纳秒级的做法。3.1 定时器触发决定了频率刻度的准确性用软件循环里调ADC_SoftwareStartConvCmd采一个点两次转换之间的间隔受中断响应、其他任务抢占影响抖动轻松到几微秒。对 10 kHz 采样率来说这就是百分之几的周期误差频率读数直接跟着漂。正确做法是让定时器在更新事件时输出 TRGOADC 配置为外部触发上升沿启动硬件自动完成“到点就采”CPU 完全不参与。采样率 fs 由定时器周期决定fs TIM_CLK / (PSC1) / (ARR1)。72 MHz 时钟、预分频到 1 MHz、ARR 设成 99得到 10 kHz。这个 fs 必须和代码里传给 FFT 的频率刻度保持一致任何一处改了指另一处没改频率读数就会整体缩放——这是排查时第一个要核对的地方。#define ADC_BUF_LEN 2048 /* 双通道交织每通道 1024 点 */ volatile uint16_t adc_buf[ADC_BUF_LEN]; void sample_chain_init(uint32_t fs_hz) { /* 1. TIM2 更新事件作为 TRGO决定采样节拍 */ TIM_TimeBaseInitTypeDef tb; RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE); TIM_TimeBaseStructInit(tb); tb.TIM_Prescaler 72 - 1; /* 72 MHz - 1 MHz */ tb.TIM_Period (1000000 / fs_hz) - 1; /* 决定 fs */ tb.TIM_CounterMode TIM_CounterMode_Up; TIM_TimeBaseInit(TIM2, tb); TIM_SelectOutputTrigger(TIM2, TIM_TRGO_Update); TIM_Cmd(TIM2, ENABLE); /* 2. ADC1 规则组外部触发源选 TIM2_TRGO */ ADC_InitTypeDef adc; RCC_APB2PeriphClockCmd(RCC_APB2Periph_ADC1, ENABLE); ADC_StructInit(adc); adc.ADC_ContinuousConvMode DISABLE; /* 单次由触发驱动 */ adc.ADC_ExternalTrigConv ADC_ExternalTrigConv_T2_TRGO; adc.ADC_NbrOfChannel 2; /* 双通道交替扫描 */ ADC_Init(ADC1, adc); ADC_RegularChannelConfig(ADC1, ADC_Channel_0, 1, ADC_SampleTime_28Cycles); ADC_RegularChannelConfig(ADC1, ADC_Channel_1, 2, ADC_SampleTime_28Cycles); ADC_ExternalTrigConvCmd(ADC1, ENABLE); /* 3. DMA1 通道1 循环搬运半满/全满各触发一次中断 */ DMA_InitTypeDef dma; RCC_AHBPeriphClockCmd(RCC_AHBPeriph_DMA1, ENABLE); DMA_StructInit(dma); dma.DMA_PeripheralBaseAddr (uint32_t)ADC1-DR; dma.DMA_MemoryBaseAddr (uint32_t)adc_buf; dma.DMA_DIR DMA_DIR_PeripheralSRC; dma.DMA_BufferSize ADC_BUF_LEN; dma.DMA_Mode DMA_Mode_Circular; dma.DMA_IT_TC | dma.DMA_IT_HT; DMA_Init(DMA1_Channel1, dma); DMA_Cmd(DMA1_Channel1, ENABLE); ADC_DMACmd(ADC1, ENABLE); ADC_Cmd(ADC1, ENABLE); }关键参数逐个说。ADC_ExternalTrigConv_T2_TRGO让转换起点由硬件给出CPU 只负责事后收数据。ADC_NbrOfChannel 2配成扫描模式两个通道在同一触发下依次转换但它们并非严格同时刻采样——这个细节在第 5 章会引发相位差问题。DMA_Mode_Circular让缓冲区循环复用配合半满和全满两个中断就能实现“算前一半时采后一半”的双缓冲避免读写同一块内存。3.2 双缓冲与数据就绪的边界单缓冲最常见的写法是等 DMA 全满中断进去把 1024 点拷出来再调 FFT。问题是从中断触发到拷贝完成这段时间里DMA 已经开始覆盖缓冲区头部拷出来的数据前半段是上一轮、后半段是下一轮频谱直接裂开。双缓冲把 2048 个 uint16_t 分成两半半满中断处理前半、全满中断处理后半处理期间另一半正在被填充互不干扰。中断服务函数里只做一件事置标志位。FFT 和后续计算都放到主循环否则中断执行时间过长会拖慢采样节拍甚至丢触发。常见的实现思路是维护一个volatile uint8_t half_ready状态机主循环轮询到标志就取走对应半区取完清标志。3.3 从 12 位码值到 Q15 的映射ADC 输出是无符号 12 位FFT 输入要求有符号定点。转换分三步减直流偏置、转有符号、定标。偏置取 2048半量程减完得到 [-2048, 2047]再乘一个定标系数搬到 Q15 的高位。/* 把交织缓冲区里指定通道的 1024 点转成 Q15 实序列 */ void pack_channel_q15(const uint16_t *src, int ch, q15_t *dst, uint16_t n) { for (uint16_t i 0; i n; i) { int32_t v (int32_t)src[i * 2 ch] - 2048; /* 去直流偏置 */ v 3; /* 2048*8 16384留半量程余量 */ if (v 32767) v 32767; if (v -32768) v -32768; dst[i] (q15_t)v; } }左移三位而不是四位是给蝶形运算的中间累加留出 2 倍余量。如果输入信号本来就只有满量程的一半可以改成左移四位提升信噪比代价是大幅值信号会削顶削顶在频域表现为谐波和底噪抬升比量化噪声更难排查。4. 幅值、频率、相位差三个测量量各自的血泪细节FFT 输出的是复数数组索引 k 对应频率 k·fs/N。从数组里读一个峰值就能得到频率、幅值、相位看起来是三个mag[k]和atan2的事。实际测量里这三个量各自的误差来源完全不同必须分开处理。4.1 幅值归一化系数与窗函数增益补偿单频信号的频谱是一根谱线但只有信号频率正好落在 k·fs/N 上时才成立。不落在整数 bin 上时能量会泄漏到邻近若干 bin这叫频谱泄露。为了压制泄露通常会加窗汉宁窗最常用代价是主瓣变宽峰值幅值也被窗函数的相干增益压低。幅值计算分两步。先算复数的模mag sqrt(re² im²)。再乘归一化系数单边谱要乘 2加窗后还要除以窗函数的相干增益。汉宁窗的相干增益是 0.5所以总系数是 4矩形窗是 1总系数是 2。这个系数如果漏乘幅值会恒定偏小一半而且不随频率变化很容易被误判成前端增益问题。/* 从 1024 点频谱中提取单频分量的幅值、频率、相位 */ void analyze_tone(const q15_t *xr, const q15_t *xi, float fs, tone_t *out) { float mag[FFT_N / 2]; uint16_t kmax 1; for (uint16_t k 1; k FFT_N / 2; k) { float re xr[k] / 32768.0f; float im xi[k] / 32768.0f; mag[k] sqrtf(re * re im * im); if (mag[k] mag[kmax]) kmax k; } /* 抛物线插值修正栅栏效应delta 落在 [-0.5, 0.5] */ float a mag[kmax - 1], b mag[kmax], c mag[kmax 1]; float delta 0.5f * (a - c) / (a - 2.0f * b c); out-freq_hz (kmax delta) * fs / FFT_N; /* 单边谱 x2汉宁窗再除以 0.5 相干增益 */ float scale 2.0f / 0.5f; out-amp_vpp mag[kmax] * scale; /* 相位取插值后主峰处的 atan2 */ out-phase_rad atan2f(xi[kmax], xr[kmax]); }scale这个系数是本段的重点除以 0.5 就是补偿汉宁窗的相干增益乘 2 是把双边谱折成单边谱的物理幅值。窗口类型一换这个系数必须跟着换否则幅值读数会成比例失真。有些实现图省事直接乘 4凑巧对上汉宁窗换窗后就全错。4.2 频率栅栏效应与插值校正频率分辨率是 fs/N。10 kHz 采样、1024 点每个 bin 是 9.77 Hz。如果信号是 1000 Hz它落在第 102.4 个 bin 上峰值会出现在第 102 和第 103 个 bin 之间直接读 kmax 得到的是 996.1 Hz 或 1005.9 Hz误差最大可达半个 bin。抛物线插值用峰点及左右邻点拟合一条抛物线把峰位偏移 delta 估出来精度能提升一个数量级。前提是频谱确实像抛物线——加窗不当或者信号含两个相近频率时插值反而会引入偏差。更严谨的做法是 Jacobsen 估计或者基于相位差的频率估计但计算量更大在 MCU 上要算一笔账。还有一种“免费”提高分辨率的方法补零。把 1024 点补到 4096 点再做 FFTbin 间隔变成 fs/4096。但补零不增加信息量只是把频谱插值得更密主瓣宽度不变对实际频率精度的提升有限主要作用是让插值更容易定位峰位。4.3 相位差两路 FFT 的相位相减为什么不准测两路同频信号的相位差直觉做法是分别做两次 1024 点 FFT各自取峰值 bin 的atan2再相减。这个做法在仿真里很准在实际硬件上经常偏出一个固定角度而且这个偏差随频率变化。根因在采样时刻。标准库下配置 ADC 双通道扫描时通道 1 和通道 2 是同一个触发沿启动但它们之间隔了一次完整的转换时间采样 转换几十个 ADC 时钟周期。这个时间差对 1 kHz 信号来说很小等于相位偏移2π·f·Δt但对 10 kHz 甚至更高频率偏移就不可忽略。频率越高固定时间差引入的相位误差越大这正是“偏差随频率变化”的来源。解决方向有两条。一是改用双 ADC 同步规则模式两个 ADC 在同一触发下各自采样同一时刻的不同通道时间差压到纳秒级。二是承认这个时间差在软件里按2π·f·(n·T_adc)做补偿其中 T_adc 是单次转换周期。前者硬件成本高一点但一劳永逸后者只在频率已知且不高的场景下可靠。另外峰值 bin 的相位对频率偏移很敏感信号频率不完全落在 bin 上时相邻 bin 的相位会围绕真实值旋转插值位置稍有偏差相位就跟着偏。更稳的办法是用插值后的频率和复数谱在峰点附近做最小二乘拟相位或者用两路信号的互谱X1·conj(X2)取相位互谱法对幅度不敏感只依赖相位一致性。5. 避坑与排查五条把结果带偏的常见问题5.1 频率读数整体偏移、幅值随频率变化现象输入 1 kHz读数稳定在 980 Hz 附近换 5 kHz 后误差更大幅值读数在低频段正常高频段明显偏低。原因分两类。频率整体偏移通常是 fs 标称值和实际值不符——定时器预分频或 ARR 算错或者晶振本身偏差。幅值随频率下降是前端抗混叠滤波器的滚降造成的滤波器在通带边缘已经衰减了几个 dB但代码里仍按理想通带归一化。解决先用已知频率的信号源反推实际 fs即fs_real fs_nominal × f_measured_known / f_read核对定时器参数幅值补偿则要在标定环节测出每个频点的衰减系数建一张小表做修正或者把测量范围限制在滤波器平坦段。5.2 两路相位差恒定偏十几度现象两路同相正弦输入相位差读数稳定在 12° 左右改变信号频率后偏角跟着变。原因ADC 双通道扫描不是同时采样通道之间存在固定的转换时间差如 4.3 节所述这个时间差乘上信号角频率就是相位误差。另外如果两路走线长度差异大也会叠加一个固定延时。解决换双 ADC 同步模式或在软件里按已知 Δt 做相位补偿Δφ 2π·f·Δt。补偿前先确认 Δt 是“采样保持时间 转换时间”的总和只算转换时间会补偿不足。5.3 偶发数据整段为 0 或频谱乱跳现象串口打印的 1024 点里有整段 0或者某次 FFT 结果完全对不上重跑一次又正常。原因DMA 缓冲区被覆盖。用了单缓冲主循环取数据时 DMA 已经回卷到头部或者中断里置了标志但主循环处理太慢下一轮中断到来时数据已经被改写。解决改成双缓冲 半满/全满中断主循环处理一块时另一块在填充同时把 FFT 耗时测量出来确保“处理一块的时间 填充一块的时间”否则要降采样率或减点数。5.4 频谱上出现不该有的谐波和底噪现象输入是纯正弦频谱里却能看到 2 次、3 次谐波底噪也偏高。原因ADC 输入过载削顶或者输入阻抗不匹配导致波形失真也可能是 Q15 定标时左移太多蝶形运算中间结果溢出产生回绕。解决示波器看 ADC 引脚波形有没有顶部削平把 Q15 定标系数降低一档观察谐波是否消失确认蝶形加减法之前有没有做整体右移保护。5.5 直流附近有一根巨大的谱线现象k0 处幅值远大于信号峰且随输入直流偏置变化。原因去偏置没做或偏置值不准。ADC 的实际零点可能不是 2048受基准电压和前端失调影响。解决采一段静态输入求平均值作为真实偏置替代固定 2048同时在 FFT 前对数据做一次均值扣除双重保险。注意扣除均值会轻微改变频谱形状对极低频测量要谨慎。6. 设计取舍定点还是浮点、点数选多少、精度换什么把这条链路走完会发现每一个参数都是拿一种资源换另一种资源。点数 N 决定分辨率和耗时。1024 点、10 kHz 采样分辨率 9.77 Hz频率精度靠插值能到 1 Hz 量级换成 4096 点分辨率提到 2.44 Hz但内存占用翻四倍FFT 耗时大约变成 4.5 倍N log N 的比例。在 SRAM 只有几十 KB 的型号上4096 点 Q15 复数数组就要 32 KB加上旋转因子表和双缓冲很容易顶到天花板。我一般会先把 N 定在 1024 或 2048用插值补精度而不是盲目加点数。定点还是浮点本质是速度和动态范围的交换。Q15 定点在无 FPU 的型号上比软件浮点快三到五倍但小信号叠大直流时量化噪声压不住带 FPU 的型号跑浮点 FFT1024 点约 1~2 ms动态范围直接拉满代价是代码体积和功耗略高。如果测量对象是固定幅度的工频或音频信号定点足够如果是宽动态范围的传感器输出浮点省心得多。窗函数和插值是另一组交换。矩形窗主瓣最窄、频率分辨率最高但泄露严重幅值读数对频偏极其敏感汉宁窗泄露小、幅值稳定但主瓣变宽两个相近频率分不开。测单一正弦、要求幅值准选汉宁窗加相干增益补偿要求频率分辨率、信号又是整周期采样矩形窗反而更好。插值校正能补频率精度但它假设主瓣形状规整多频信号下会失效这时候要么加点数要么改用更复杂的频率估计算法。相位差这一环最容易被低估。软件补偿能救急但补偿系数依赖对 ADC 时序的准确建模换一款芯片或改一次采样时间就得重标定。我现在的习惯是只要涉及两路相位优先确认硬件是否支持同步采样不支持就在设计阶段把测量频率上限压下来让固定时间差引入的相位误差落在可接受范围内而不是事后靠公式硬掰。从那以后我每次搭这套链路都强制走一遍“已知信号源标定 → 核对 fs → 查削顶 → 验相位”的流程任何一步跳过后面调参都是盲调。希望这套拆解能帮你在自己的板子上少绕几圈——文末附带的工程里把采样链路、定点 FFT、插值校正和双缓冲都拆成了独立文件可以直接对照着改参数验证。本文还有配套的精品资源点击获取