1. 项目概述在STM32上实现高精度相位测量最近在做一个电机控制相关的项目需要精确测量两路正弦信号的相位差。用传统的过零检测法精度受限于采样率和噪声在低频时误差很大。后来想到既然MCU已经采集了完整的波形数据为什么不直接用FFT快速傅里叶变换把信号的频谱和相位信息直接“算”出来呢这比硬件比较要灵活和精确得多。STM32的Cortex-M系列内核特别是M3/M4/M7都支持DSP指令集ST官方也提供了完善的DSP库。其中就包含了高度优化的实数FFT函数我们完全可以利用它来做一个高精度的“软件相位计”。这个思路不仅适用于电机控制在电力谐波分析、音频处理、振动监测等需要分析信号间时序关系的场景下都非常有用。今天我就把自己在STM32F103C8T6也就是常说的“蓝桥杯”核心板上利用标准外设库和DSP库实现FFT测相位的完整过程、踩过的坑以及关键技巧详细地梳理一遍。目标是让你看完就能在自己的板子上复现并且理解每一个参数背后的意义。2. 核心思路与方案选型为什么是FFT以及如何选型2.1 从过零检测到频谱分析的思维转变传统的相位测量尤其是对于工频50Hz信号很多工程师的第一反应是过零比较。方法很简单用比较器或者ADC采样后软件判断捕捉两个信号从负到正或从正到负跳变的时刻计算时间差再换算成相位差。这个方法直观但有几个硬伤抗噪能力差信号上的毛刺或轻微的失真都可能导致错误的过零点判断。依赖波形对称性如果信号含有直流偏移或偶次谐波过零点会偏移。低频测量精度低对于低频信号一个采样点的误差就会带来很大的相位误差。而FFT方法是一种“全局”分析方法。它不对波形上的单个点做判断而是对一整段采样数据比如256个点进行数学变换得到这段数据中各个频率分量的幅度和相位信息。对于纯净的正弦波FFT后会在其频率对应的“频点”上出现一个峰值这个峰值对应的复数结果实部虚部就包含了该频率信号的幅度和初始相位。通过计算两个信号在同一频率点上的初始相位再做差就能得到相位差。这种方法本质上是对一段时间内的信号整体特征进行提取抗噪声和干扰的能力强得多精度也由采样率和FFT点数决定可以做到很高。2.2 STM32 DSP库的优势与限制ST为Cortex-M4/M7等提供了CMSIS-DSP库这是一个由ARM主导的、高度优化的DSP函数库。对于Cortex-M3如STM32F103ST提供了一个较早的STM32_DSP_Lib库其中也包含了FFT函数。我们项目里用的就是这个。为什么选择它速度极快库函数用汇编和SIMD指令高度优化比我们自己用C语言写的FFT算法快一个数量级以上。这对于实时性要求高的应用至关重要。接口稳定函数接口简单专注于输入输出方便集成。资源明确官方示例给出了不同FFT点数64 256 1024所需的内存RAM大小方便我们根据手头芯片的资源配置进行选择。需要注意的限制定点运算这个老版本的DSP库主要处理Q15或Q31格式的定点数。这意味着我们需要把ADC采样的整数结果转换成库函数要求的定点格式FFT计算完成后再转换回我们能理解的幅度和相位值。这个过程涉及到**缩放因子Scaling Factor**的理解是新手最容易出错的地方。仅支持2的整数次幂点数库函数只支持64 256 1024等点数。点数越多频率分辨率越高但计算量也越大所需RAM也越多。实数FFT库函数cr4_fft_xxx_stm32执行的是实数FFT。它输入N个实数采样点输出N/2个复数频点包含正频率部分。这正好符合我们对实信号分析的需求比复数FFT效率更高。2.3 关键参数设计采样率、点数与频率分辨率这是一个核心的工程设计环节参数选不好结果要么不准要么算不出来。FFT点数NPT我们选择256点。为什么对于STM32F103C8T620K RAM1024点FFT需要约22K RAM如代码注释警告很可能导致内存溢出。64点分辨率又太低。256点是一个在精度和资源消耗之间很好的平衡点。采样率Fs代码中定义为FreqSample其值为50*NPT 50*256 12800 Hz。这里的设计非常巧妙目标频率我们想观测50Hz的信号工频。频率分辨率FFT的频率分辨率 Δf Fs / NPT。如果我们想让50Hz信号恰好落在某个频点上避免频谱泄漏就需要满足50Hz k * Δf k * (Fs / NPT)其中k是整数。推导将上式变形为 Fs 50 * NPT / k。为了计算方便通常取k1即一个信号周期内恰好采样NPT个点。但这要求Fs50*NPT12800Hz对于很多应用来说采样率过高。更常见的做法是取k为NPT的约数例如kNPT/2则Fs100Hz但这需要信号是严格的周期性且同步采样。代码中的设计代码采用了Fs 50 * NPT这意味着它模拟了一个“理想采样”场景采样率是信号频率的整数倍256倍并且采样起始点与信号相位是同步的这由MygSin函数中固定的初始相位196.0*PI/180.0决定。在实际项目中我们通常无法控制采样与信号的起始相位同步因此这个高采样率主要是为了演示计算的方便。在实际应用中Fs应根据奈奎斯特定律Fs 2*信号最高频率和所需频率分辨率Δf Fs/NPT来综合确定。缩放因子ScalingFactor设置为20000。这是因为库函数使用的Q15格式表示范围是[-1, 1)对应整数[-32768, 32767)。我们的正弦波幅度想设为接近满量程的1那么就需要乘以32767。这里取20000是为了留出余量防止运算溢出同时也模拟了ADC采样值不会总是满量程的情况。注意这些参数是演示代码设定的“理想情况”。在实际硬件ADC采样时Fs由定时器触发ADC的频率决定ScalingFactor就是ADC的实际采样值如12位ADC值范围0-4095需转换为有符号数用于FFT。理解这个对应关系是项目成功的关键。3. 代码深度解析与实操要点让我们跳出代码顺序抓住几个最关键的、容易混淆的函数进行解剖。3.1 信号生成函数MygSin理解模拟采样的构造这个函数用于生成测试用的正弦波数据模拟了ADC采样的过程。理解它就理解了如何准备FFT的输入数据。void MygSin(long nfill, long Fs, long Freq1, long Freq2, long Ampli) { uint32_t i; float fFs, fFreq1, fFreq2, fAmpli; float fZ,fY; fFs (float) Fs; fFreq1 (float) Freq1; fFreq2 (float) Freq2; fAmpli (float) Ampli; for (i0; i nfill; i) { // 核心公式生成正弦波样本 fY PERCENT_DC/2 // 一个很小的直流分量 sin((196.0*PI/180.0) PI2 * i * (fFreq1/fFs)); if (fFreq2) fY 0.2*sin((293.0*PI/180.0) PI2 * i * (fFreq2/fFs)); fZ fAmpli * fY; // 乘以缩放因子模拟ADC量化 // 关键操作将实数样本放入复数格式的输入数组 lBUFIN[i] ((short)fZ) 16; } }关键点解析初始相位196.0*PI/180.0和293.0*PI/180.0是直接设定的初始相位弧度制。这告诉我们FFT计算出的相位是相对于这段数据采样窗口起始时刻的初始相位。复数输入格式lBUFIN数组每个元素是long型32位。库函数要求输入是复数数组即使我们的信号是实数。这里采用了一种紧凑格式高16位存放实部低16位存放虚部。对于纯实数输入虚部为0。所以((short)fZ) 16就是将实数样本fZ作为实部放在高16位低16位默认为0。直流分量PERCENT_DC/2添加了一个很小的直流偏移。在FFT结果中直流分量会体现在第0个频点lBUFMAG[0]上。代码注释也提到了这一点。实操要点当你用真实的ADC数据替换这个函数时需要将ADC值比如0-4095转换为有符号整数。通常做法是int16_t adc_val (raw_adc - 2048);将其零中心化然后像代码一样lBUFIN[i] ((int16_t)adc_val) 16;。确保你的采样值在[-32768, 32767]范围内否则在后续定点运算中会溢出。3.2 FFT核心调用与幅值计算powerMag这是整个流程的核心。先调用FFT再计算每个频点的幅度。void powerMag(long nfill, char* strPara) { int32_t lX,lY; uint32_t i; for (i0; i nfill/2; i) { // 1. 从复数输出中提取实部(lX)和虚部(lY) lX (lBUFOUT[i] 16) 16; // 取低16位对应cos实部这里需要辨析 lY (lBUFOUT[i] 16); // 取高16位对应sin虚部 // 2. 转换为浮点数并进行缩放 float X (float)lX * (NPT / 32768.0); float Y (float)lY * (NPT / 32768.0); // 3. 计算幅度谱 float Mag sqrt(X*X Y*Y) / nfill; lBUFMAG[i] (uint32_t)(Mag * 65536); } }这里有一个巨大的坑也是很多初学者直接抄代码会出错的地方。注意看提取lX和lY的代码lX (lBUFOUT[i] 16) 16;这是一个先左移16位再右移16位的操作其效果是取低16位并做符号扩展。根据注释“sine_cosine -- cos”它被认为是实部。lY (lBUFOUT[i] 16);直接右移16位取高16位。注释认为是虚部。但是这与MygSin函数中存入数据的方式((short)fZ) 16矛盾了在MygSin中我们把数据放在了高16位实部低16位是0虚部。那么按照这个逻辑在powerMag中我们应该用lY (lBUFOUT[i] 16)来取实部用lX (lBUFOUT[i] 16) 16来取虚部才对。然而代码的注释写的是sine_cosine -- cos和sine_cosine -- sin。在数字信号处理中一个复数通常表示为a bj其中a是实部对应cos分量b是虚部对应sin分量。FFT的输出也是这种格式。经过对STM32 DSP库源码的查阅和实际测试可以确认库函数cr4_fft_xxx_stm32的输出格式是每个long型输出数据中高16位存放的是实部Re低16位存放的是虚部Im。因此代码中powerMag函数里的提取操作是正确的但注释可能引起了误解。正确的理解是lX (lBUFOUT[i] 16) 16;// 提取低16位即虚部(Im)。lY (lBUFOUT[i] 16);// 提取高16位即实部(Re)。这导致了后续X和Y的赋值与变量名相反。在计算Mag sqrt(X*X Y*Y)时因为平方和与顺序无关所以幅度计算是正确的。但在计算相位时如果弄反了实部和虚部结果将是错误的实操心得对于ST的DSP库记住这个约定输出数组lBUFOUT中每个元素的高16位是实部(Re)低16位是虚部(Im)。自己写代码时变量名最好用Re和Im避免用X/Y或cos/sin造成混淆。缩放因子NPT / 32768.032768对应Q15格式的1.0。NPT是FFT点数这个缩放是为了将定点数结果转换为实际的幅度值。最后除以nfill是FFT的常规归一化操作。3.3 相位计算函数power_Phase_Radians核心中的核心这是计算相位的关键函数也是错误的重灾区。void power_Phase_Radians(long nfill) { int32_t lX,lY; uint32_t i; for (i0; i nfill/2; i) { // 提取操作与powerMag相同 lX (lBUFOUT[i] 16) 16; // 虚部 Im lY (lBUFOUT[i] 16); // 实部 Re { // 缩放与powerMag一致 float X NPT * ((float)lX) / 32768; // 注意这里X对应的是虚部Im float Y NPT * ((float)lY) / 32768; // 这里Y对应的是实部Re // 计算相位角 atan(Y/X) 注意这里用的是Y和X即 atan(Re / Im) float phase atan(Y / X); // 象限判断与修正 if (Y 0) { // 实部Re 0 if (X 0) { // 虚部Im 0 ; // 第一象限phase已在[-π/2, π/2]内atan结果在[0, π/2]正确。 } else { // 虚部Im 0 phase PI; // 第二象限atan结果为负需加π转到[π/2, π] } } else { // 实部Re 0 if (X 0) { // 虚部Im 0 phase PI2; // 第四象限atan结果为负需加2π转到[3π/2, 2π) } else { // 虚部Im 0 phase PI; // 第三象限atan结果为正需加π转到[π, 3π/2] } } // 转换为角度制存储 lBUFPHASE[i] phase * 180.0 / PI; } } }现在我们来分析最关键的问题相位计算公式和象限修正。标准的复数相位计算公式是φ atan2(Im, Re)即atan2(虚部 实部)。C语言标准库atan2(y, x)函数返回的是y/x的反正切其值域是(-π, π]它已经自动处理了象限。而代码中计算的是phase atan(Y/X);其中YRe,XIm。这等价于atan(Re / Im)也就是atan(实部/虚部)。这与标准的atan2(虚部实部)是不同的。因此后面的象限判断逻辑if (Y0)...是基于YRe和XIm来写的它是在手动实现atan2(Re, Im)的功能而不是atan2(Im, Re)。这会导致什么结果计算出的相位phase与理论值会有一个90度π/2的固定偏差。因为atan(Re/Im)与atan(Im/Re)的关系是互余的相差90度再加上象限修正逻辑的差异最终结果可能不是简单的偏移。正确的做法应该是使用标准库函数atan2f直接避免象限判断的麻烦和错误。确保传入的参数顺序是(虚部 实部)。修正后的核心代码段应为// 正确提取实部(Re)和虚部(Im) int32_t lRe (lBUFOUT[i] 16); // 高16位是实部 int32_t lIm (lBUFOUT[i] 16) 16; // 低16位是虚部 // 转换为浮点数并缩放 float Re NPT * ((float)lRe) / 32768.0f; float Im NPT * ((float)lIm) / 32768.0f; // 使用atan2f计算相位参数顺序为(Im, Re) float phase_rad atan2f(Im, Re); // 返回值范围 [-π, π] // 如果需要转换为0到360度 if (phase_rad 0) { phase_rad 2 * PI; } float phase_deg phase_rad * 180.0f / PI;实操要点与避坑指南永远使用atan2f(Im, Re)这是计算复数相位的标准方法安全可靠。理解相位含义FFT计算出的相位是该频率分量正弦波在采样窗口起始时刻t0的初始相位。如果你采样窗口的开始时间相对于信号是随机的那么这个初始相位也是随机的。但是两个信号做FFT后在同一频率点上计算出的初始相位之差就是它们之间的相位差这个值与采样窗口的起始时刻无关是稳定的。相位缠绕Phase Wrappingatan2f返回的范围是[-π, π]。当你计算相位差时如果结果接近±180°可能会因为噪声从179°跳变到-181°实际上它们只差2°。在后续处理中比如用于锁相环可能需要使用phase_unwrap函数进行解缠绕。4. 完整实操流程与移植指南现在我们将理论转化为在你自己板子上可运行的步骤。假设你使用的是STM32F103C8T6核心板并已经配置好了基本的工程模板包含系统时钟、延时、串口等。4.1 步骤一获取并添加DSP库获取库文件从ST官网搜索“STM32 DSP Lib”或在你使用的STM32CubeIDE/Keil MDK的安装目录里寻找。通常文件名是STM32_DSP_Lib_Vx.x.xx.rar。解压后我们主要需要Lib文件夹下的.a或.lib库文件以及Examples文件夹下的头文件和示例源码。添加库到工程在工程目录下新建一个文件夹例如Drivers/CMSIS/DSP。将库文件如arm_cortexM3l_math.lib注意l代表小端模式复制到该文件夹。将必要的头文件如stm32_dsp.h,table_fft.h也复制过来。table_fft.h里面定义了FFT运算所需的旋转因子表非常重要。工程配置添加头文件路径在IDE的工程属性中将Drivers/CMSIS/DSP添加到包含路径Include Paths。添加库文件在链接器Linker设置中添加库文件路径和库名如-larm_cortexM3l_math。开启硬件FPU如果适用对于Cortex-M4F/M7确保在编译选项中开启了-mfpufpv4-sp-d16或类似选项并添加宏定义ARM_MATH_CM4或ARM_MATH_CM7。对于M3我们使用纯软件浮点或定点库。4.2 步骤二配置ADC与定时器进行双路同步采样这是实际项目中最关键的一步。我们需要用定时器触发ADC以固定频率对两路信号进行同步采样。定时器配置以TIM2为例配置为向上计数模式预分频器和自动重载值ARR根据你的采样率Fs计算。例如系统时钟72MHz要得到12.8kHz的采样率则定时器频率应为72MHz / (PSC1) / (ARR1) 12800Hz。可以设置PSC71ARR78。ADC配置使用ADC1的两个通道如通道0和通道1。工作模式选择为“双重模式”下的“同步规则模式”或“交替模式”确保两路信号采样时刻尽可能接近。更精确的做法是使用带注入组的双ADC模式。触发源选择为“定时器2的触发输出TRGO”。开启DMA。DMA应配置为循环模式从ADC数据寄存器如ADC1-DR搬运到内存中的两个缓冲区adc_buf1[NPT]和adc_buf2[NPT]。数据缓冲区在内存中定义两个int16_t类型的数组adc_buf1[NPT]和adc_buf2[NPT]。DMA填满这两个缓冲区后产生半满或全满中断在中断中设置一个标志位通知主循环可以进行FFT计算。4.3 步骤三主循环中的FFT计算流程在主循环中等待DMA采集完成标志然后进行以下操作// 1. 数据预处理将ADC值转换为DSP库需要的Q15格式 for(int i0; iNPT; i) { // 假设ADC是12位原始值0-4095。先零中心化再缩放到Q15范围附近。 // 注意防止溢出缩放因子要小于32767/4096 ≈ 8 int16_t val1 ((int16_t)adc_buf1[i] - 2048) * 8; // 缩放因子8满幅值约为±16384 int16_t val2 ((int16_t)adc_buf2[i] - 2048) * 8; // 存入复数输入数组实部为ADC值虚部为0 lBUFIN_1[i] (val1 16) 0xFFFF0000; // 高16位为实部 lBUFIN_2[i] (val2 16) 0xFFFF0000; // 低16位为0 } // 2. 执行FFT cr4_fft_256_stm32(lBUFOUT_1, lBUFIN_1, NPT); cr4_fft_256_stm32(lBUFOUT_2, lBUFIN_2, NPT); // 3. 计算幅度和相位使用修正后的函数 calc_mag_and_phase(lBUFOUT_1, mag1, phase1, NPT); calc_mag_and_phase(lBUFOUT_2, mag2, phase2, NPT); // 4. 寻找目标频率点例如50Hz int target_bin (int)(50.0 * NPT / Fs 0.5); // 计算50Hz对应的频点索引 // 确保target_bin在[1, NPT/2-1]范围内0是直流分量 // 5. 计算相位差 float phase_diff phase2[target_bin] - phase1[target_bin]; // 处理相位环绕 if(phase_diff 180.0f) phase_diff - 360.0f; if(phase_diff -180.0f) phase_diff 360.0f; // 6. 通过串口打印或使用 printf(Freq: %.1f Hz, Mag1: %.2f, Mag2: %.2f, Phase Diff: %.2f deg\r\n, (float)target_bin * Fs / NPT, mag1[target_bin], mag2[target_bin], phase_diff);4.4 步骤四优化与校准加窗处理如果采样不是整周期即频率不是频率分辨率的整数倍会发生频谱泄漏导致幅值和相位计算不准。解决方法是对采样数据加窗如汉宁窗Hamming。在数据预处理循环中将ADC值乘以窗函数系数window[i]然后再进行格式转换。幅值校准FFT计算出的幅值需要乘以一个系数才能反映真实幅值。对于矩形窗这个系数是2/NPT因为能量分散到两个频点。对于汉宁窗系数约为4/NPT因为主瓣宽度变宽。需要通过标准信号源进行校准。相位校准由于ADC通道、运放等硬件通路可能存在微小的延时差异即使输入同相信号测出的相位差也可能不为零。这需要引入一个“系统相位偏移”进行软件补偿。方法是输入两个同源同相的信号测量出一个固定的相位差offset然后在后续测量中减去这个offset。5. 常见问题、调试技巧与实测记录在实际操作中你几乎一定会遇到下面这些问题。这里是我的排查记录和解决方法。5.1 问题一FFT结果全是噪声或幅度极小现象计算出的mag数组所有值都很小没有明显的峰值。排查检查输入数据在调用FFT前通过串口打印出lBUFIN数组的前几个值。看看是否是你预期的ADC转换后的值。确认数据零中心化是否正确即没有直流偏移时数据应在0上下波动。检查数据格式确认lBUFIN中数据是否按照(实部16)的格式存放。可以用调试器查看内存。检查库链接确认DSP库是否正确链接。可以尝试调用一个简单的DSP函数如arm_add_f32测试。检查缩放因子ScalingFactor或ADC值的缩放倍数是否太大导致溢出尝试减小它。解决最常见的原因是输入数据格式错误或数值溢出。确保ADC值在转换为Q15格式前其绝对值最大值小于32767。5.2 问题二相位差结果跳动很大不稳定现象幅值测量基本准确但计算出的相位差在正负几十度之间随机跳动。排查频谱泄漏这是首要原因。如果你的信号频率不是Fs/NPT的整数倍能量会泄漏到相邻频点导致主频点的相位测量不准。观察频谱打印出mag数组看看目标频率点是不是唯一的、尖锐的峰值。如果旁边有较高的旁瓣说明泄漏严重。信噪比低如果信号本身噪声大或幅度太小噪声会严重影响相位计算。没有整周期采样确保采样点数NPT包含整数个信号周期。可以通过调整采样率Fs或使用更大的NPT来提高频率分辨率使目标频率更接近某个频点中心。相位解缠绕问题如果相位差在180度附近可能会因为计算误差在179和-181之间跳变。确保你使用了正确的atan2f和相位差范围处理规整到-180到180度。解决强制整周期采样对于已知频率的信号如工频50Hz使用锁相环PLL或高精度定时器使采样率动态跟踪信号频率确保每个采样块都是整数个周期。这是最有效的方法。加窗如果无法实现整周期采样必须加窗汉宁窗、汉明窗。加窗后相位信息在频点中心仍然是准确的但需要重新校准幅值系数。插值在最大幅值频点附近进行插值如重心法可以更精确地估计真实频率和相位。5.3 问题三测量结果存在固定偏差现象测量一个已知为0度或90度相位差的信号结果总是差一个固定值。排查实部虚部弄反再次检查power_Phase_Radians函数中lX和lY的提取是否正确以及atan2f的参数顺序是否为(Im, Re)。这是最可能的原因。硬件通道延时两路信号的模拟前端运放、滤波电路如果不完全一致会产生固定的时间差从而引入固定的相位差。用同源信号测试即可测出这个系统误差。ADC采样时刻不同步如果两路ADC不是严格同步采样例如用了交替模式而非同步模式会引入一个采样时间间隔的相位差。这个差值与频率成正比Δφ 360° * Δt * f。计算一下这个值是否与你的偏差吻合。解决修正代码中的实部虚部错误。进行系统校准输入同相信号记录测量出的相位差φ_offset在后续所有测量中减去此值φ_true φ_measured - φ_offset。5.4 调试技巧实录可视化是王道不要只盯着一个相位差数值看。把整个mag数组和phase数组通过串口发送到电脑用PythonMatplotlib或MATLAB画出来。一眼就能看出频谱是否干净峰值是否明显相位曲线是否平滑。从理想信号开始先用MygSin函数生成理想的数字正弦波进行FFT验证你的算法和库是否正确。得到正确结果后再接入真实的ADC数据。这样可以隔离硬件问题。分步验证第一步验证ADC采样值是否正确。输入一个直流电压看ADC值是否对。第二步验证定时器触发频率是否正确。用示波器看ADC的转换完成信号或DMA传输信号。第三步验证FFT输入缓冲区数据是否正确。在内存中查看lBUFIN的值是否是你采集的波形。第四步验证FFT输出和幅度计算。给一个单频正弦波看频谱是否只有一个尖峰幅度是否与输入信号电压成比例。第五步最后再验证相位。利用断点和内存查看在调用FFT函数前和后设置断点直接查看lBUFIN和lBUFOUT数组的内存内容。对比理论值能快速定位是数据准备问题还是库函数计算问题。移植这个FFT测相位功能到实际项目就像在调试一个精密的仪器。最重要的不是记住代码而是理解数据在每个环节的形态从模拟信号到ADC数字量到Q15定点数到FFT复数输出再到最终的幅度和相位。其中任何一个环节的缩放、格式或顺序出错都会导致结果面目全非。我的经验是把power_Phase_Radians函数里的实部虚部提取和atan2计算单独拿出来用一组已知的复数如10j, 01j, -10j等进行单元测试确保它能算出正确的0°, 90°, 180°角。这个基础打牢了后面的工作就顺利多了。