简介面向嵌入式信号处理初学者与需要做频域测量的工程师这份资料围绕STM32F407平台用标准库实现FFT并完成正弦波幅值、频率与相位差的测量。内容涉及ADC采样、窗函数处理、FFT变换、频谱分析与结果输出适合电机控制、音频分析和无线通信等实时信号处理场景的学习与验证。压缩包共283个文件约9.98MB包含C/H源文件、Keil工程配置、编译中间文件及可执行输出等便于直接打开工程、对照源码理解算法流程与工程结构。已有1668人学习下载可从中获得一套较完整的STM32 FFT测量实现参考并据此梳理采样、频谱计算和结果输出之间的衔接关系。1. 一块最小系统板加两路正弦波输入就能得到这样的读数一块 STM32 最小系统板两路同频正弦波分别从 PA0、PA1 输入串口每隔 100 ms 输出一行结果A1.002V F1000.31Hz dPhi44.97deg——这是我把 FFT 真正跑在单片机上之后一直关注的东西。很多人对快速傅里叶变换的印象还停留在上位机里那条漂亮的频谱曲线可一旦搬到 STM32 上难点根本不在蝶形公式本身而在采样率怎么定、直流偏置怎么去掉、频率没有正好落在某个频率间隔上时幅值为什么忽大忽小、两个通道不是同一时刻采样时相位差为什么会偏。这套基于标准外设库的 FFT 测参方案做的就是整条链路定时器触发 ADC 采 1024 点一次基 2 变换把频域里的峰值还原成幅值、频率和双通道相位差。它适合手上已经有 F103 或 F407 板子、想把示波器那点测量能力加入自己产品的开发者也适合刚啃完 DSP 教材、想找个能上板验证的练手项目的人。2. 信号链与采样前端把连续正弦波变成 FFT 能处理的数组采样前端这一层最容易被轻视。算法写得再漂亮进来的数组本身就是乱的后面全是白费功夫。这一章先把参数算清楚再把定时器、ADC、DMA 这套配置搭出来最后做码值到浮点数组的预处理。2.1 采样率、点数和分辨率三者互相制约给定采样率fs和点数N频率分辨率Δf fs / N一帧的时间长度T N / fs能无混叠测量的最高频率是fs / 2。测 1 kHz 正弦取fs 20 kHz、N 1024得到Δf ≈ 19.53 Hz、T 51.2 ms。这 51.2 ms 意味着一帧内 1 kHz 信号转了 51 个周期这一点很关键周期数越多主瓣越窄幅值受频谱泄漏影响越小。反过来只采 5 个周期就算加了窗幅值读数也会随信号频率抖动这类精度问题很容易被忽略。参数典型取值决定什么FFT_N1024频率分辨率与运算量1024 点基 2 是 5120 次蝶形FS_HZ20000可测上限 10 kHz触发间隔 50 µsΔf19.53 Hz频率的量化单位T51.2 ms数据刷新率的上限ADC 时钟12 MHzF103 上不得超过 14 MHz采样率不能只满足奈奎斯特实际工程里我一般会留出 510 倍余量理由是两个一是抗混叠滤波器不可能做成理想的矩形滚降需要过渡带二是每周期点数太少时ADC 的孔径抖动和触发抖动会被放大成相位噪声。20 kHz 采 1 kHz 得到每周期 20 点配上 1024 点一帧算是精度和运算量之间比较舒服的一档。抗混叠那一级千万别省。PA0 前面串一个 1 kΩ 电阻、对地挂 33 nF截止频率约 4.8 kHz正好压在信号带宽之上、fs/2之下把 10 kHz 以上会折叠回来的成分先衰减掉。没有这一级示波器探头夹上去看到的“干干净净的正弦”在频谱上照样能冒出一堆莫名其妙的谱线。注意ADC 采样时间要跟信号源内阻匹配。41.5 个采样周期在 12 MHz 下是 3.46 µs源阻抗几 kΩ 以内都够用如果前面加了很大的分压电阻采样电容充不满幅值会系统性偏小而且偏得很“稳定”特别容易误判成算法问题。2.2 定时器触发 ADC 扫描 DMA 搬运的标准库写法采样必须由硬件定时器触发绝不能在主循环里手动启动转换否则采样间隔抖动直接变成相位误差。标准外设库这套初始化比较长但每一步都有明确目的照着改宏定义就能用。#define FFT_N 1024u #define FS_HZ 20000.0f #define VREF 3.3f #define ADC_FULL 4095.0f uint16_t adc_raw[FFT_N * 2]; /* 两通道交替存放CH0,CH1,CH0,CH1... */ volatile uint8_t frame_ready 0; void sample_init(void) { GPIO_InitTypeDef gpio; ADC_InitTypeDef adc; DMA_InitTypeDef dma; TIM_TimeBaseInitTypeDef tim; RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOA | RCC_APB2Periph_ADC1, ENABLE); RCC_AHBPeriphClockCmd(RCC_AHBPeriph_DMA1, ENABLE); RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE); RCC_ADCCLKConfig(RCC_PCLK2_Div6); /* 72MHz/6 12MHz别超 14MHz */ gpio.GPIO_Pin GPIO_Pin_0 | GPIO_Pin_1; gpio.GPIO_Mode GPIO_Mode_AIN; GPIO_Init(GPIOA, gpio); ADC_DeInit(ADC1); adc.ADC_Mode ADC_Mode_Independent; adc.ADC_ScanConvMode ENABLE; /* 两通道轮流转换 */ adc.ADC_ContinuousConvMode DISABLE; /* 由 TIM2 触发不能连续 */ adc.ADC_ExternalTrigConv ADC_ExternalTrigConv_T2_TRGO; adc.ADC_DataAlign ADC_DataAlign_Right; adc.ADC_NbrOfChannel 2; ADC_Init(ADC1, adc); ADC_RegularChannelConfig(ADC1, ADC_Channel_0, 1, ADC_SampleTime_41Cycles5); ADC_RegularChannelConfig(ADC1, ADC_Channel_1, 2, ADC_SampleTime_41Cycles5); ADC_DMACmd(ADC1, ENABLE); ADC_ExternalTrigConvCmd(ADC1, ENABLE); DMA_DeInit(DMA1_Channel1); dma.DMA_PeripheralBaseAddr (uint32_t)ADC1-DR; dma.DMA_MemoryBaseAddr (uint32_t)adc_raw; dma.DMA_DIR DMA_DIR_PeripheralSRC; dma.DMA_BufferSize FFT_N * 2; /* 半字为单位 2*N */ dma.DMA_PeripheralInc DMA_PeripheralInc_Disable; dma.DMA_MemoryInc DMA_MemoryInc_Enable; dma.DMA_PeripheralDataSize DMA_PeripheralDataSize_HalfWord; dma.DMA_MemoryDataSize DMA_MemoryDataSize_HalfWord; dma.DMA_Mode DMA_Mode_Circular; dma.DMA_Priority DMA_Priority_High; DMA_Init(DMA1_Channel1, dma); DMA_ITConfig(DMA1_Channel1, DMA_IT_TC, ENABLE); DMA_Cmd(DMA1_Channel1, ENABLE); /* 72MHz / 20000Hz 3600 个计数向上计数更新事件作 TRGO */ tim.TIM_Period 3600 - 1; tim.TIM_Prescaler 0; tim.TIM_ClockDivision TIM_CKD_DIV1; tim.TIM_CounterMode TIM_CounterMode_Up; TIM_TimeBaseInit(TIM2, tim); TIM_SelectOutputTrigger(TIM2, TIM_TRGO_Update); TIM_Cmd(TIM2, ENABLE); ADC_Cmd(ADC1, ENABLE); ADC_ResetCalibration(ADC1); while (ADC_GetResetCalibrationStatus(ADC1)); ADC_StartCalibration(ADC1); while (ADC_GetCalibrationStatus(ADC1)); ADC_SoftwareStartConvCmd(ADC1, ENABLE); /* 首次启动后续靠 TRGO */ }几个参数值得单独说。TIM_Period 3600 - 1是 72 MHz 除以 20 kHz 的结果改采样率就改这个数改完记得同步FS_HZ宏两处不一致时频率读数会整体缩放这种错误很典型。ADC_ExternalTrigConv_T2_TRGO配上TIM_SelectOutputTrigger(TIM2, TIM_TRGO_Update)才是完整的一对只写一半的话 ADC 会一直等不到触发现象是 DMA 一个数都不搬。DMA_BufferSize的单位是传输元素个数而不是字节配置成半字的时候就是2 * FFT_N写成4 * FFT_N会多搬一倍数组越界却不报错。还有一点扫描模式下 CH0 和 CH1 不是同时采的CH1 比 CH0 晚了“采样时间 转换时间”这么久——41.5 12.5 54 个 ADC 周期在 12 MHz 下正好 4.5 µs。这个固定偏斜后面算相位差时必须补回来第 4 章会用到。2.3 码值到浮点数组去直流比标定更要紧转换这一步做两件事把 12 位码值映射成电压把直流分量减掉。第二件事如果忘了0 号频点的能量会大到把自动量程、峰值搜索全部带跑偏。static void raw_to_float(const uint16_t *raw, float *x0, float *x1) { float sum0 0.0f, sum1 0.0f; for (uint16_t i 0; i FFT_N; i) { x0[i] raw[2 * i] * (VREF / ADC_FULL); /* CH0 */ x1[i] raw[2 * i 1] * (VREF / ADC_FULL); /* CH1 */ sum0 x0[i]; sum1 x1[i]; } float dc0 sum0 / (float)FFT_N; /* 均值即直流分量 */ float dc1 sum1 / (float)FFT_N; for (uint16_t i 0; i FFT_N; i) { x0[i] - dc0; x1[i] - dc1; } }用均值当直流分量是最省事的做法前提是一帧里包含整数个信号周期否则均值本身会带一点残余。测 1 kHz、采 51.2 ms 这种整周期关系时误差可以忽略信号频率未知或者不是整周期就该换成硬件上的中点偏置——两路信号先经过 1.65 V 的偏置网络再进 ADC把 0 V 抬到量程中央然后减固定常数VREF/2比减均值稳定得多。硬件上还得多一句STM32 的 ADC 输入不能承受负电压正弦波必须抬到 03.3 V 之间。信号源直接输出 ±1 V 就进管脚轻则波形下半周削平、频谱上出现一堆奇次谐波重则把管脚打坏。我一般用一只 10 kΩ 电阻把信号拉到 1.65 V 分压点串 100 nF 隔直简单可靠。3. FFT 核心在 STM32 上跑 1024 点基 2 快速傅里叶变换前端交出两个干净的浮点数组接下来是纯算法部分。1024 点基 2 是 10 级、每级 512 个蝶形一共 5120 次蝶形运算。手写一遍的价值在于你能确切知道每个下标在干什么出了问题也能定位到具体哪一级。3.1 位反转置换与旋转因子表基 2 的时域抽取DIT要求输入按位反转顺序排列输出才是自然顺序。位反转这段代码短但下标逻辑很容易写错写成死循环或者只对一部分元素生效都不报错只是频谱结果完全不对。static float fft_re[FFT_N], fft_im[FFT_N]; static float tw_re[FFT_N / 2], tw_im[FFT_N / 2]; static void fft_init_twiddle(void) { for (uint16_t i 0; i FFT_N / 2; i) { float ang -2.0f * 3.14159265358979f * (float)i / (float)FFT_N; tw_re[i] cosf(ang); tw_im[i] sinf(ang); } } static void bit_reverse(float *re, float *im, uint16_t n) { uint16_t j 0; for (uint16_t i 1; i n; i) { uint16_t bit n 1; while (j bit) { /* 找到 j 的下一个反转值 */ j ^ bit; bit 1; } j ^ bit; if (i j) { /* 只在前者小于后者时交换保证每对换一次 */ float t re[i]; re[i] re[j]; re[j] t; t im[i]; im[i] im[j]; im[j] t; } } }if (i j)这个判断是必须的去掉它会把每对元素交换两次等于什么都没做而数组里只有实部非零、虚部全 0 时错误又不会崩溃属于典型的“静默错误”。旋转因子我放在const大数组里烧到 Flash 更省 RAM1024 点需要 512 个float共 2 KBF103C8 那点 20 KB RAM 里塞得下但也不宽裕可以按需只存第一象限再用对称性推出其余三个象限。3.2 三层循环的蝶形运算蝶形那段是整个算法最容易写错下标的地方。以 DIT 为例外层按级数增长块长度len中间层遍历所有块内层在块内做一半的蝶形。void fft_radix2(void) { bit_reverse(fft_re, fft_im, FFT_N); for (uint16_t len 2; len FFT_N; len 1) { uint16_t half len 1; uint16_t step FFT_N / len; /* 旋转因子表的下标步长 */ for (uint16_t i 0; i FFT_N; i len) { for (uint16_t j 0; j half; j) { uint16_t k i j; /* 上半个蝶形 */ uint16_t m k half; /* 下半个蝶形 */ uint16_t ti j * step; /* 取表下标不是 j */ float wr tw_re[ti], wi tw_im[ti]; float mr fft_re[m], mi fft_im[m]; float tr mr * wr - mi * wi; /* 复数乘 W */ float tv mr * wi mi * wr; fft_re[m] fft_re[k] - tr; fft_im[m] fft_im[k] - tv; fft_re[k] tr; fft_im[k] tv; } } } }ti j * step是最关键的一行。每一级的旋转因子只用到表里的一部分第一级len 2只用W[0]最后一级len N要用完整 512 个step就是用来按级抽取的。写成ti j的话第一级就会取到W[1]、W[2]结果看着像频谱、频率也对得上但幅值和相位全都慢慢偏掉——这类错误最难排查因为波形“看起来是对的”。先做减法再赋值也需要注意顺序fft_re[k]被更新之前必须先算出tr上面代码已经先存了中间变量如果把fft_re[k] tr写在前面fft_re[m]那一行用的就是脏数据。RAM 方面三个 1024 点浮点数组实部、虚部加上输入的两路就是 12 KB 以上F103C8 的 20 KB 一下就去掉大半。真要在小容量芯片上跑输入可以直接原地变换两路信号分时处理反正 FFT 是线性变换两路之间没有耦合。3.3 浮点、Q15 与 CMSIS-DSP 的三条路同等正确性下摆在面前的选择有三条。第一条是上面这种手写单精度浮点代码透明改参数方便缺点是没有 FPU 的 Cortex-M3 上全靠软件浮点库指令数大约是整数的三到五倍。我在 72 MHz 的 F103 上跑 1024 点单精度量级落在几十毫秒做 10 Hz 级的慢速监测完全够做实时音频就吃力了。带 FPU 的 F407 是另一个数量级1 ms 以内能出结果代价是芯片贵一截、功耗高一点。第二条是 Q15 定点。所有数据缩放到 ±1.0 对应的 ±32767蝶形每级的增益是 2 倍所以每做一级就要右移一位防止溢出1024 点做完 10 级整体相当于右移了 10 位小信号的信噪比掉得厉害。想要精度就得用块浮点每级统一判断是否溢出再决定移几位复杂度立刻上一个台阶。第三条是直接用 CMSIS-DSP这也是我最常用的方案。它是 ARM 的官方 DSP 库和标准外设库并不冲突Keil 工程里把arm_cortexM3l_math.libM4F 用arm_cortexM4lf_math.lib加进来就行。#include arm_math.h static arm_rfft_fast_instance_f32 rfft; static float fft_in[FFT_N], fft_out[FFT_N]; void dsp_fft_prepare(void) { arm_rfft_fast_init_f32(rfft, FFT_N); /* 只需初始化一次 */ } /* 输入实数序列输出打包格式的频谱 */ void dsp_fft_run(const float *x, float *mag, uint16_t half_n) { memcpy(fft_in, x, sizeof(float) * FFT_N); arm_rfft_fast_f32(rfft, fft_in, fft_out, 0); /* 0 正变换 */ /* 打包格式out[0]DCout[1]Nyquist 之后 out[2k]Re[k]out[2k1]Im[k]k 从 1 到 N/2-1 */ mag[0] fabsf(fft_out[0]); mag[half_n] fabsf(fft_out[1]); for (uint16_t k 1; k half_n; k) { float re fft_out[2 * k]; float im fft_out[2 * k 1]; mag[k] sqrtf(re * re im * im); } }arm_rfft_fast_f32用的是实数 FFT比起把虚部填零再跑复数 FFT运算量直接减半这个便宜不占白不占。它输出的打包格式是刚接触时最容易踩坑的地方不是[Re0, Im0, Re1, Im1, ...]那种规整排列直流和奈奎斯特两个点被单独拎出来塞在开头。照着这个顺序去取sqrt(re²im²)出来的幅值会整体错位读数完全对不上。提示arm_rfft_fast_f32只支持 32 位单精度。M3 上没有 FPU 的时候arm_rfft_q15配合arm_cfft_q15会快很多但要把缩放因子一路管到底调试期间我一般不这么干。如果只是教学或者想弄清楚每一步手写那一版更好如果项目要交付、要跑得快直接用 CMSIS-DSP。这两条路我都留着因为前者用来对照后者能很快判定问题出在算法还是出在调用。4. 三参数提取从频谱到幅值、频率、相位差频谱算出来只是起点。真正要交给上层的是三个数幅值多少伏、频率多少赫兹、两路之间相位差多少度。这三件事里只有频率是自然落在频率间隔上的另外两个都需要做额外处理。4.1 峰值搜索与单边谱幅值还原实数信号的 FFT 是共轭对称的能量在正负频率各占一半所以正频率那一侧找到的幅值要乘 2 才是真实幅值。加上窗函数之后还要除掉窗的相干增益汉宁窗是 0.5。typedef struct { float amp; /* 幅值单位 V */ float freq; /* 频率单位 Hz */ float phase; /* 初相位单位 rad */ } wave_t; void spectrum_measure(const float *x, wave_t *out) { float win_sum 0.0f; for (uint16_t i 0; i FFT_N; i) { /* 汉宁窗N-1 做分母保证首尾对称 */ float w 0.5f - 0.5f * cosf(2.0f * 3.14159265f * (float)i / (float)(FFT_N - 1)); fft_re[i] x[i] * w; fft_im[i] 0.0f; win_sum w; } fft_radix2(); uint16_t kmax 1; float mmax 0.0f; for (uint16_t k 1; k FFT_N / 2; k) { /* 跳过 DC 与奈奎斯特 */ float m sqrtf(fft_re[k] * fft_re[k] fft_im[k] * fft_im[k]); if (m mmax) { mmax m; kmax k; } } /* 单边谱 *2再除以窗相干增益 win_sum/N汉宁窗约 0.5 */ out-amp 2.0f * mmax / win_sum; out-freq (float)kmax * FS_HZ / (float)FFT_N; out-phase atan2f(fft_im[kmax], fft_re[kmax]); }2.0f * mmax / win_sum这一句把窗增益归一化合并进去了比硬写mmax / (FFT_N * 0.5f)更通用——换布莱克曼窗或者平顶窗的时候只要改窗函数不用回头改幅值公式。平顶窗的主瓣很宽频率分辨能力差但幅值精度高得多专门用来做幅值表的场合我会换成它。峰值搜索的范围同理从k 1开始、到N/2 - 1结束。包含 0 号或者N/2那个点的话只要信号里有一点残余直流峰值就会被吸到 0 号上去频率直接读成 0 Hz。4.2 频率细化插值法和相位差分法kmax * fs / N只能分辨到 19.53 Hz 这一档。想再细得用频谱插值把频率间隔内部的位置估出来。最简单的是对数幅值的抛物线插值/* 在峰值周围三点做抛物线拟合求小数偏移 delta ∈ [-1, 1] */ static float freq_delta(uint16_t k, const float *mag) { float a mag[k - 1], b mag[k], c mag[k 1]; float denom a - 2.0f * b c; if (fabsf(denom) 1e-12f) return 0.0f; /* 三点等高无法拟合 */ return 0.5f * (a - c) / denom; } float f_refined ((float)kmax freq_delta(kmax, mag)) * FS_HZ / (float)FFT_N;插值法简单、单帧就能算但精度依赖信噪比通常能做到Δf的 1/10 到 1/20也就是 12 Hz。要更好的效果就用相位差分法隔M个采样点取两帧同一频点k的相位差Δθ和频率偏移δ满足Δθ 2πδM/N于是/* theta1、theta2 是相隔 M 个采样点两帧在同一 bin 上的相位 */ float dtheta theta2 - theta1; while (dtheta 3.14159265f) dtheta - 2.0f * 3.14159265f; while (dtheta -3.14159265f) dtheta 2.0f * 3.14159265f; float delta dtheta * (float)FFT_N / (2.0f * 3.14159265f * (float)M); float f_fine ((float)kmax delta) * FS_HZ / (float)FFT_N;代价是要采两帧刷新率降一半换来的是在中等信噪比下能把误差压到 0.1 Hz 量级。M越大越灵敏但2πδM/N超过 π 就会产生相位模糊M的上限大约是N/(2|δ|max)实际取N/4比较稳妥。4.3 双通道相位差atan2、频率偏移补偿与采样偏斜两路各自用atan2f(im, re)算出峰值处的相位相减就是相位差。但这里有两处系统误差必须补。第一处来自频率没有落在频率间隔正中。设真实频率是(k δ)·fs/N那么这个 bin 上测到的相位并不是信号的真实初相而是被旋转了一个角度。矩形窗和汉宁窗在这一点上关系一致/* 把 bin k 上的相位换算回信号真实初相 */ static float phase_correct(float phi_bin, uint16_t k, float delta) { return phi_bin - 3.14159265f * delta * (float)(FFT_N - 1) / (float)FFT_N; }π · δ · (N-1)/N这个量在N 1024、δ 0.5时大约是 1.57 rad将近 90°不补的话相位差读数会随信号频率漂移看着像“相位不稳定”实际是计算漏了一步。第二处来自两通道不是同时采样。2.2 节算过扫描模式下 CH1 比 CH0 晚 4.5 µs这等价于 CH1 的相位被额外超前了 2πf·本文还有配套的精品资源点击获取
