资讯动态

STM32嵌入式电能质量监测:FFT+THD+实时频谱图三合一

发布时间:2026/10/18 14:55:22 来源:尧图企业网站定制
简介本资源是一套基于STM32F103ZET6正点原子精英版的嵌入式信号分析完整工程面向单片机开发初学者与数字信号处理进阶学习者解决嵌入式端实时FFT频谱分析、THD失真度计算及LCD图形化显示等核心问题。压缩包共229个文件含49个头文件h、48个C源码c构成主程序框架与算法模块36个目标文件o和依赖文件d/dep体现Keil MDK编译结构另有LCD驱动、USMART调试组件、定时器与Flash操作等关键外设代码整体5.28MB结构规范适合作为STM32嵌入式信号处理项目范例。已有3969人学习下载提供可直接编译运行的完整Keil工程含uvprojx、uvoptx、sct等配置涵盖从ADC采样、ARM CMSIS-DSP库调用FFT、谐波能量提取、THD百分比计算到ILI9341液晶频谱图动态刷新的全链路实现代码注释清晰模块划分明确是掌握嵌入式DSP实战能力的优质参考。1. 一个STM32上跑通的“电能质量快照”FFT频谱总谐波失真THD实时波形图三件套全在一块最小系统板上跑稳了你手头有一块带ADC和FSMC的STM32F407开发板接了个电压互感器采样市电想立刻知道当前电网有没有严重谐波、THD是否超标、5次/7次谐波有多大——别翻HAL库手册、别配CubeMX时钟树、更别搭MATLAB串口助手。这个项目就是为这种“现场拍板”场景设计的从ADC采样开始到屏幕右下角弹出红色THD8.7%警告整个链路不到300行核心C代码全部固化在单片机里上电即跑无需PC端任何辅助软件。它不是教学Demo而是某高校电力电子实验室真实用在电机驱动老化测试台上的嵌入式监测模块——采样率10kHz、1024点实数FFT、THD按IEC 61000-4-7标准计算、频谱图用ILI9341驱动的2.4寸屏滚动刷新。新手照着README改两行宏就能点亮熟手拆开fft_core.c会发现它绕过了CMSIS-DSP的复数FFT黑匣子用查表法预存旋转因子原位运算把407主频压到72MHz也能扛住10ms一帧的吞吐。如果你正被“嵌入式FFT太慢”“THD算不准”“频谱图卡顿”三个问题反复折磨这包源码就是你的后悔药。2. 工程结构与技术栈从裸机启动到GUI渲染的六层流水线这个项目没用RTOS也没上FatFS或LwIP所有模块都扎根在STM32标准外设库非HAL的土壤里。它用最朴素的分层思想把信号处理流程切成了六个可验证环节ADC采样触发 → 环形缓冲区管理 → FFT计算引擎 → THD数值提取 → 频谱图数据映射 → LCD逐行刷新。每一层都提供明确的输入输出契约比如fft_calc()函数只认int16_t *input和uint16_t len两个参数返回float *mag_spectrum绝不碰GPIO或DMA寄存器。这种设计让调试像剥洋葱当频谱图出现横纹你只需在lcd_draw_spectrum()入口打个断点看传进来的数组值是否合理若THD恒为0则直接跳进thd_calculate()检查基波幅值是否被误判成噪声底。下面拆解这六层如何咬合。2.1 启动与硬件初始化避开SysTick陷阱的精准定时项目放弃SysTick做采样定时改用TIM2的更新中断触发ADC规则组转换。原因很现实SysTick被systick_handler占着一旦开启FreeRTOS或加调试printf采样间隔立刻漂移。而TIM2配置为向上计数模式ARR7199对应10kHz采样率CC1通道捕获用于同步LCD刷新避免频谱图撕裂。// system_init.c void TIM2_Config(void) { TIM_TimeBaseInitTypeDef TIM_TimeBaseStructure; RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE); TIM_TimeBaseStructure.TIM_Period 7199; // 72MHz / (71991) 10kHz TIM_TimeBaseStructure.TIM_Prescaler 0; // 不分频 TIM_TimeBaseStructure.TIM_ClockDivision 0; TIM_TimeBaseStructure.TIM_CounterMode TIM_CounterMode_Up; TIM_TimeBaseInit(TIM2, TIM_TimeBaseStructure); TIM_ITConfig(TIM2, TIM_IT_Update, ENABLE); // 开启更新中断 TIM_Cmd(TIM2, ENABLE); }提示TIM_Period必须是7199而非7200因为计数器从0开始实际周期ARR1。曾有A同学填错导致采样率变成9998.6Hz后续FFT频点偏移0.14Hz查了三天才发现是这里少减了1。2.2 ADC采集与环形缓冲双缓冲机制防数据覆盖ADC配置为连续扫描模式规则组含1个通道PA0每次转换结束触发DMA传输。关键在adc_buffer.c里实现的双缓冲buffer_a[1024]和buffer_b[1024]交替使用。DMA完成一半中断HTI时切换到buffer_a完成中断TCI时切换到buffer_b。这样CPU处理buffer_a时DMA正往buffer_b灌新数据彻底消除丢点。// adc_buffer.c __IO uint16_t buffer_a[1024], buffer_b[1024]; __IO uint16_t *current_buffer buffer_a; __IO uint8_t buffer_flag 0; // 0: a is ready, 1: b is ready void DMA1_Channel1_IRQHandler(void) { if (DMA_GetITStatus(DMA1_IT_HT1)) { // Half Transfer DMA_ClearITPendingBit(DMA1_IT_HT1); current_buffer buffer_b; buffer_flag 1; } if (DMA_GetITStatus(DMA1_IT_TC1)) { // Transfer Complete DMA_ClearITPendingBit(DMA1_IT_TC1); current_buffer buffer_a; buffer_flag 0; } }逻辑说明buffer_flag是CPU和DMA共享的状态标志必须声明为__IOvolatile。CPU主循环中通过while(!buffer_flag)等待buffer_a就绪拿到指针后立即调用fft_calc(current_buffer)全程不加锁——因为DMA只改指针CPU只读指针无竞态。2.3 FFT计算引擎实数FFT优化版省掉一半计算量项目没用CMSIS-DSP的arm_rfft_fast_f32()而是采用专为实数序列优化的混合基FFT。核心思想是N点实数序列的DFT其输出满足共轭对称性只需计算前N/21个点。代码里fft_real.c的real_fft_1024()函数将1024点分解为32×32先做32组32点行变换再做32组32点列变换旋转因子全部预存在twiddle_table[512]数组中避免浮点sin/cos实时计算。// fft_real.c const float twiddle_table[512] { /* 生成脚本见tools/twiddle_gen.py */ }; void real_fft_1024(int16_t *input, float *output) { int i, j, k; float temp_r, temp_i; // Step 1: Bit-reversal permutation on input bit_reverse_1024(input); // Step 2: 32-point DFT on each row (32 rows) for (i 0; i 32; i) { dft_32(input[i*32], output[i*32]); } // Step 3: Multiply by twiddle factors for (i 0; i 32; i) { for (j 0; j 32; j) { k i * 32 j; temp_r output[k] * twiddle_table[j] - output[k1] * twiddle_table[j1]; temp_i output[k] * twiddle_table[j1] output[k1] * twiddle_table[j]; output[k] temp_r; output[k1] temp_i; } } // Step 4: 32-point DFT on each column (32 columns) for (j 0; j 32; j) { dft_32_col(output[j], 32, 32); } }参数说明input是1024点int16_t原始采样值已减去直流偏置output是1024点float复数结果实部虚部交错存储。dft_32()是手工展开的32点蝶形运算比通用FFT快3倍。注意twiddle_table索引步长是2因为每项存实部虚部一对值。2.4 THD计算模块严格遵循IEC 61000-4-7的基波识别逻辑THD不是简单算(谐波有效值/基波有效值)IEC标准要求基波频率必须在49.5~50.5Hz工频或59.5~60.5Hz美标窗口内锁定谐波只计入第2~40次对应100Hz~2kHz且基波幅值取该频带内最大谱线而非理论50Hz点。thd_calculate.c里thd_calculate()函数先扫描mag_spectrum[1..40]跳过DC分量找到最大值索引fund_idx再确认其对应频率freq fund_idx * 10因采样率10kHz/1024≈9.766Hz/点是否在49.5~50.5Hz内。若不在自动扩展搜索窗到±2点仍失败则报错BASE_FREQ_LOCK_FAIL。// thd_calculate.c float thd_calculate(float *mag_spectrum) { uint16_t fund_idx 0; float fund_mag 0.0f; float harmonic_sum 0.0f; // Step 1: Find fundamental peak in 49.5~50.5Hz window (idx 5~6) for (uint16_t i 5; i 6; i) { if (mag_spectrum[i] fund_mag) { fund_mag mag_spectrum[i]; fund_idx i; } } // Step 2: If not found, expand search to idx 4~7 if (fund_mag 1e-3f) { for (uint16_t i 4; i 7; i) { if (mag_spectrum[i] fund_mag) { fund_mag mag_spectrum[i]; fund_idx i; } } } // Step 3: Calculate THD sqrt(sum(harmonics^2)) / fundamental for (uint16_t i 2; i 40; i) { if (i ! fund_idx) { harmonic_sum mag_spectrum[i] * mag_spectrum[i]; } } return sqrtf(harmonic_sum) / (fund_mag 1e-6f); // 防除零 }注意mag_spectrum[i]是FFT后取模得到的幅值谱已做sqrt(real^2 imag^2)。代码末尾的1e-6f是工程上常用的防除零技巧比if(fund_mag0)更鲁棒——浮点比较永远别用。2.5 频谱图映射对数压缩颜色分级让微弱谐波也可见原始幅值谱动态范围太大强基波vs弱5次谐波可能差60dB直接画线会丢失细节。spectrum_map.c采用分段对数映射0~10dB线性10~60dB对数60dB饱和。输出数组spectrum_mapped[320]对应LCD宽320像素每像素代表1024/320≈3.2个频点用滑动平均合并。// spectrum_map.c void spectrum_map_to_lcd(float *mag_spectrum, uint8_t *spectrum_mapped) { float max_val 0.0f; for (int i 1; i 40; i) { // 只映射前40次谐波 if (mag_spectrum[i] max_val) max_val mag_spectrum[i]; } if (max_val 1e-4f) return; for (int x 0; x 320; x) { uint16_t start_idx (x * 1024) / 320; uint16_t end_idx ((x 1) * 1024) / 320; float avg_mag 0.0f; for (uint16_t i start_idx; i end_idx i 40; i) { avg_mag mag_spectrum[i]; } avg_mag / (end_idx - start_idx); // Log mapping: 0~10dB linear, 10~60dB log, 60dB clamp float db_val 20.0f * log10f(avg_mag / max_val 1e-6f); if (db_val 0.0f) db_val 0.0f; if (db_val 10.0f) { spectrum_mapped[x] (uint8_t)(db_val * 2.55f); // 0~255 } else if (db_val 60.0f) { spectrum_mapped[x] (uint8_t)(10.0f (db_val - 10.0f) * 4.25f); } else { spectrum_mapped[x] 255; } } }逻辑说明spectrum_mapped[x]值域0~255后续lcd_draw_spectrum()将其转为RGB565颜色如0黑255红。log10f()调用前加1e-6f防log(0)这是嵌入式日志压缩的血泪经验——某次传感器断线导致全零输入没加这行直接触发FPU异常。2.6 LCD驱动与刷新FSMCDMA双加速频谱图丝滑滚动ILI9341屏通过FSMC总线挂载数据线D0~D15接FSMC_D0~D15。关键优化在lcd_refresh.c频谱图不用逐像素写而是把320×240的映射数组frame_buffer[320*240]预先生成好再用DMA2 Channel5将整块内存刷到LCDGRAM。lcd_draw_spectrum()只负责更新frame_buffer的底部一行240行高每次刷新底部10行其余区域保持静止降低带宽压力。// lcd_refresh.c uint16_t frame_buffer[320*240]; // 全局帧缓冲 void lcd_draw_spectrum(uint8_t *spectrum_mapped) { // Update bottom 10 lines of frame_buffer for (int y 230; y 240; y) { for (int x 0; x 320; x) { uint8_t val spectrum_mapped[x]; // Map 0~255 to RGB565: 00x0000(black), 2550xF800(red) uint16_t color (val 3) 11; // R5 frame_buffer[y * 320 x] color; } } // Trigger DMA transfer to LCDGRAM LCD_SetCursor(0, 230, 319, 239); // Set window to bottom 10 lines LCD_WriteRAM_Prepare(); DMA2_Channel5-CMAR (uint32_t)frame_buffer 230 * 320 * 2; // Start addr DMA2_Channel5-CNDTR 320 * 10; // Transfer 3200 words DMA_Cmd(DMA2_Channel5, ENABLE); }参数说明LCD_SetCursor()设置LCDGRAM写入窗口LCD_WriteRAM_Prepare()发送WR信号DMA2将frame_buffer的指定区域搬运过去。CNDTR是16位计数器单位是uint16_t所以320×103200。实测此方案比CPU逐字写快12倍频谱图滚动延迟8ms。3. 编译与烧录Keil MDK-ARM v5.37下的零配置落地项目工程文件为.uvprojx格式直接用Keil uVision5打开即可编译。但有几个关键配置点必须手动核对否则会链接失败或FFT结果全零——这些不是文档里写的“建议”而是某开发者踩坑后写死在工程注释里的硬约束。3.1 启动文件与堆栈设置禁用microlib堆栈大小精确到字节工程使用startup_stm32f407xx.s启动文件但默认的__initial_sp指向0x2001FFFF128KB SRAM顶而项目实际只用SRAM1112KB需手动修改; startup_stm32f407xx.s Stack_Size EQU 0x00001000 ; 原为0x00002000过大导致malloc失败 AREA STACK, NOINIT, READWRITE, ALIGN3 Stack_Mem SPACE Stack_Size __initial_sp EQU 0x2001C000 ; 原为0x2001FFFF改为SRAM1末尾提示Stack_Size0x10004KB是实测最小安全值。若开更多调试printf需同步增大但超过0x2000会导致heap区被挤占——heap起始地址在0x2001C000 - 0x1000 0x2001B000大小固定为0x2000。3.2 C/C选项浮点单元与优化等级的生死线在Options for Target → C/C页必须勾选Use MicroLIB取消勾选microlib不支持log10f等高级数学函数Use FPU选择FPv4-SP-D16F407的FPU型号Optimization设为Level 3-O3但必须添加--fpmodefastKeil特有启用快速浮点Additional C/C flags: --fpmodefast -DARM_MATH_CM4 -D__FPU_PRESENT1原因--fpmodefast允许编译器用VSQRT.F32指令替代软件开方arm_sqrt_f32()执行时间从840周期降到12周期。某次未加此flagTHD计算耗时从1.2ms暴涨到9.7ms直接导致10ms帧率崩盘。3.3 Linker设置自定义分散加载确保FFT表进FLASHtwiddle_table[512]必须放在FLASH区ROM否则上电后值为随机数。在Options for Target → Linker页勾选Use Memory Layout from Target Dialog然后在Scatter File中添加LR_IROM1 0x08000000 0x00080000 { ; load region size_region ER_IROM1 0x08000000 0x00080000 { ; load address execution address *.o (RO) .twiddle_table_section RO ; 关键强制twiddle_table进ROM } RW_IRAM1 0x20000000 0x00020000 { ; RW data *.o (RW ZI) } }并在fft_real.c顶部添加段声明// fft_real.c #pragma push #pragma location.twiddle_table_section const float twiddle_table[512] { /* ... */ }; #pragma pop3.4 调试配置SWD速度与半主机禁用Options for Target → Debug页选择ST-Link DebuggerSettings → Trace页SWD Clock设为4000 kHz过高易丢包过低下载慢Enable SWO取消勾选SWO占用PB3与LCD背光冲突Use Semihosting取消勾选半主机依赖ARM标准库与microlib冲突注意若勾选Semihostingprintf(THD%.2f%%\r\n, thd*100)会卡死在__sys_write因为项目没实现_sys_open等底层函数。3.5 烧录验证三步确认法5分钟定位90%问题烧录后若屏幕黑或THD恒为0按顺序执行查串口1PA9/PA10用USB-TTL接电脑波特率115200看是否输出[INFO] ADC init OK。若无输出检查USART1_Config()中RCC_APB2PeriphClockCmd(RCC_APB2Periph_USART1, ENABLE)是否漏写。测PA0电压万用表测PA0对地电压应为1.65V±0.1VADC参考电压VREFINT1.2V经分压。若为0V检查互感器供电是否接反。抓TIM2波形示波器接PA1TIM2_CH1应看到10kHz方波。若无检查TIM2_Config()中RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE)是否漏写。4. 常见问题排查那些让你怀疑人生的5个瞬间项目虽小但嵌入式信号处理的坑又深又窄。以下是某开发者在实验室白板上记下的5条血泪记录每一条都对应一次凌晨三点的万用表测量。4.1 现象频谱图左右颠倒基波峰出现在最右侧原因bit_reverse_1024()函数中位反转逻辑错误将0b00000000001反转成0b10000000000512但实际应为0b100000000001024点下0b00000000001→0b100000000001024错1024点索引0~10230b000000000011反转10位得0b1000000000512。解决重写bit_reverse_1024()用查表法替代实时计算const uint16_t bit_rev_table[1024] {0,512,256,768,...}; // 预生成1024项 for (int i 0; i 1024; i) { uint16_t rev bit_rev_table[i]; if (i rev) { int16_t temp input[i]; input[i] input[rev]; input[rev] temp; } }4.2 现象THD值随采样率变化10kHz时为5%5kHz时突变为12%原因thd_calculate()中基波搜索窗硬编码为i5~6对应48.8~58.6Hz但5kHz采样率下频点间隔5000/1024≈4.88Hz50Hz对应索引50/4.88≈10.2应在i10~11。解决动态计算搜索窗uint16_t fund_start (49.5f * 1024.0f) / sample_rate; // 49.5Hz对应索引 uint16_t fund_end (50.5f * 1024.0f) / sample_rate; for (uint16_t i fund_start; i fund_end; i) { ... }4.3 现象LCD频谱图出现规律性竖条纹间隔32像素原因spectrum_map_to_lcd()中start_idx和end_idx计算用整数除法1024/3203导致每32像素合并的频点数不均前319像素各3点最后1像素剩1点avg_mag计算失真。解决改用浮点累加按比例分配权重float weight 1024.0f / 320.0f; for (int x 0; x 320; x) { float start_f x * weight; float end_f (x 1) * weight; int start_i (int)floorf(start_f); int end_i (int)ceilf(end_f); // 然后对[start_i, end_i)内每个i按重叠长度加权 }4.4 现象烧录后屏幕全白但串口有打印原因ILI9341的RESET引脚PC0在LCD_Init()中被拉低后未释放或PC0配置为推挽输出但未初始化为高电平。解决在LCD_GPIO_Config()中强制置位GPIO_ResetBits(GPIOC, GPIO_Pin_0); // 先拉低 Delay_ms(10); GPIO_SetBits(GPIOC, GPIO_Pin_0); // 再拉高 Delay_ms(120); // 等待LCD稳定4.5 现象FFT结果全为0但ADC采样值正常原因real_fft_1024()函数中bit_reverse_1024()操作的是int16_t *input但dft_32()期望float *类型不匹配导致内存越界覆写。解决统一数据类型在fft_real.c顶部加转换void real_fft_1024(int16_t *input_i16, float *output) { static float input_f32[1024]; for (int i 0; i 1024; i) { input_f32[i] (float)input_i16[i] / 32768.0f; // 归一化到[-1,1] } // 后续所有计算基于input_f32 }5. 二次开发实战从“能跑”到“能改”的三条路径拿到这个项目新手的目标是让它在自己的板子上亮起来而熟手的价值在于把它变成自己项目的子模块。我经历过三次改造第一次给电机驱动加谐波监测第二次集成到光伏逆变器并网质检仪第三次移植到GD32E507上跑通。每一次都印证了一个原则不要动FFT核心只动数据管道。下面给出三条经过验证的改造路径每条都附可立即粘贴的代码片段。5.1 路径一接入新传感器——替换ADC通道5分钟完成你的板子用PB1接电流霍尔传感器而非原项目的PA0电压。改动仅三处修改ADC1_Config()中通道ADC_RegularChannelConfig(ADC1, ADC_Channel_9, 1, ADC_SampleTime_15Cycles);PB1对应ADC_Channel_9修改dma_config()中DMA地址DMA_SetCurrDataCounter(DMA2_Channel1, 1024); DMA_SetAddress(DMA2_Channel1, (uint32_t)ADC1-DR, (uint32_t)buffer_a, 1024);在main()循环中fft_calc()前加直流偏置校准// 每10秒采1024点算均值作为offset static uint32_t offset_cal_counter 0; if (offset_cal_counter 100) { uint32_t sum 0; for (int i 0; i 1024; i) sum buffer_a[i]; adc_offset sum / 1024; offset_cal_counter 0; } // 在fft_calc前减偏置 for (int i 0; i 1024; i) { buffer_a[i] - adc_offset; }5.2 路径二扩展谐波分析——从40次到50次改两行宏IEC 61000-4-7:2002要求分析至50次谐波2.5kHz。只需改thd_calculate.c中两个宏#define MAX_HARMONIC_ORDER 50 // 原为40 #define SPECTRUM_LENGTH 512 // 原为40改为512因50次谐波需索引50但FFT输出0~512然后在thd_calculate()循环中for (uint16_t i 2; i MAX_HARMONIC_ORDER; i)。注意mag_spectrum数组需扩大到float mag_spectrum[512]否则越界写入破坏栈。5.3 路径三移植到GD32E507——时钟与外设重映射表GD32E507与STM32F407引脚兼容但寄存器略有差异。关键差异如下表其他代码0修改模块STM32F407寄存器GD32E507寄存器修改说明ADC采样时间ADC_SampleTime_15CyclesADC_SAMPLETIME_13POINT5GD32用13.5周期替代15周期FSMC地址建立FSMC_Bank1_NORSRAMInitStruct.FSMC_AddressSetupTime 0x02FSMC_Bank1_NORSRAMInitStruct.FSMC_AddressSetupTime 0x03GD32建立时间多1周期FPU使能SCB-CPACR 0x00F00000SCB-CPACR时钟使能RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOA, ENABLE)rcu_periph_clock_enable(RCU_GPIOA)GD32用rcu_前缀移植后system_clock_config()中将RCC_PLLMul_9改为RCU_PLL_MUL9其余不变。实测GD32E507在120MHz主频下FFT耗时比F407快18%因GD32的乘法器是单周期。从那以后我每次接到新MCU平台都强制走一遍“寄存器映射表核对”流程——哪怕文档说100%兼容也要拿示波器抓TIM2波形、用逻辑分析仪看FSMC时序。因为信号处理模块的脆弱性不在算法而在时序毛刺引发的采样抖动而抖动不会报错只会让THD值每天漂移0.3%直到客户投诉“你们的设备不准”。希望帮到你。本文还有配套的精品资源点击获取

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑