电能采集128点FFT谐波分析详解
目录一、头文件 Dl645_Fft.h 完整注释二、源文件 Dl645_Fft.c 逐段详细注释三、核心逻辑分层总结便于阅读理解1. 数据链路2. FFT 核心两段逻辑3. 状态机分工无阻塞 MCU4. 工程适配重点电能采集芯片128 点基 2 DIT-FFT 谐波分析分文件注释重点标注蝶形运算、倒位序、采样转换、谐波计算逻辑。一、头文件 Dl645_Fft.h 完整注释#ifndef _DL645_Fft_H #define _DL645_Fft_H // FFT固定点数 1282^7 基2时分抽取FFT #define NUM_FFT 128 // 圆周率用于预生成旋转因子sin/cos表 #define PI 3.1415926 /** * brief FFT运算输出结果结构体 * details 存储实虚部、基波有效值、1~50次谐波含有率、总畸变THD */ typedef struct { float iRealArray[NUM_FFT]; // FFT输出频谱实部数组 float iMageArray[NUM_FFT]; // FFT输出频谱虚部数组 u16 FU[3]; // 三相基波电压有效值 UA/UB/UC u32 FI[4]; // 三相基波电流有效值 IA/IB/IC // HarmonicpercentU[x][0] THD总谐波畸变率 // HarmonicpercentU[x][n] n次电压谐波含有率(放大1e6存储小数) u32 HarmonicpercentU[3][51]; // HarmonicpercentI[x][0] THD总谐波畸变率 // HarmonicpercentI[x][n] n次电流谐波含有率 u32 HarmonicpercentI[3][51]; }sDl645FftHarmonicData_TypeDef; /** * brief 瞬时采样缓存结构体 * details 存储RN8302读取的原始24bit ADC采样、状态机标志、通道地址 */ typedef struct { u8 Channel; // 当前正在计算的通道 0~5 UA/UB/UC/IA/IB/IC u16 ReadAdress; // RN8302瞬时采样数据寄存器起始地址 u8 StarFlag; // 谐波分析状态机标志位 // 0空闲1等待采样缓存完成2执行FFT谐波计算 u8 ADSPIBusy; // SPI读取采样忙标志 u16 DataCount; // 已读取采样点计数 u16 ReadAddres; // 备用采样地址变量(冗余) u16 dwFreq; // 当前电网频率 u16 TimeOutStamp; // 采样启动时间戳用于延时等待缓存 u8 InstantaneousData[384]; // 128点×3字节384字节 24bit原始采样缓存 } sADE_Instantaneous_TypeDef; /** * brief 6路采集通道枚举定义 */ typedef enum { UAChannel 0x00, // A相电压 UBChannel 0x01, // B相电压 UCChannel 0x02, // C相电压 IAChannel 0x03, // A相电流 IBChannel 0x04, // B相电流 ICChannel 0x05 // C相电流 }eADEChannelFlag_TypeDef; // 外部函数声明 /** * brief 谐波分析状态机调度入口主循环1ms周期调用 */ extern void fnDl645Fft_Harmonic_Exec(void); /** * brief FFT模块初始化清空缓存、复位状态机 */ extern void fnDl645Fft_init(void); // 预留定点FFT接口(本工程未使用) extern void Int_FFT(s16 ReArray[], s16 ImArray[]); // 预留定点倒位序接口(本工程未使用) extern void Bit_Reverse(s16 BR_Array[]); #endif二、源文件 Dl645_Fft.c 逐段详细注释/**************************************************************************//** * file Dl645_Fft.c * brief RN8302三相电压电流128点DIT-FFT谐波分析 * 计算1~50次谐波含有率THD总谐波畸变率 * version V1.2 * date 2017/5/26 * note * Copyright (C) REnergy Micro-Electronics Co.,LTD. *--------------------------Update History------------------------------------- * V1.1 新增SPI批量读取采样缓存fnRN8302_ReadBuf * V1.2 修正总谐波畸变ftemp_P平方累加计算错误 ******************************************************************************/ #define _DL645_fft_GLOBALS #include App.h /** * brief 1~50次谐波硬件补偿系数表 * details RN8302ADC前端增益随频率衰减FFT幅值计算后乘以该系数补偿误差 * Fftcoefficient[0] 2次谐波补偿Fftcoefficient[49]50次谐波补偿 */ const float Fftcoefficient[50] { 1.0005, 1.0027, 1.0054, 1.0087, 1.0129, 1.0175, 1.0229, 1.0290, 1.0362, 1.0438, 1.0524, 1.0618, 1.0718, 1.0827, 1.0944, 1.1075, 1.1210, 1.1360, 1.1517, 1.1685, 1.1862, 1.2051, 1.2255, 1.2475, 1.2700, 1.2942, 1.3199, 1.3468, 1.3758, 1.4070, 1.4421, 1.4714, 1.5107, 1.5483, 1.5873, 1.6312, 1.6727, 1.7199, 1.7700, 1.8231, 1.8797, 1.9448, 2.0074, 2.0740, 2.1534, 2.2126, 2.2919, 2.3762, 2.4670, 2.5639, }; /** * brief SPI批量读取RN8302 16点24bit瞬时采样数据 * param wReg RN8302采样缓存寄存器起始地址 * param pBuf 数据接收缓冲区单次读取16点×3字节48字节 * note 带校验和校验校验失败丢弃本次采样数据 */ void fnRN8302_ReadBuf(u16 wReg,u8 *pBuf) { u8 i,j,temp,chksum; // 配置RN8302 SPI引脚模式 PinMode_ADCS(GPIO_MODE_OUT); // 片选输出 PinMode_ADSCLK(GPIO_MODE_OUT); // 时钟输出 PinMode_ADSDO(GPIO_MODE_OUT); // MCU→芯片输出 PinMode_ADSDI(GPIO_MODE_IN); // 芯片→MCU输入 ENABLE_RN8302(); // 拉低片选选中RN8302 // 发送寄存器低8位地址 temp (u8)(wReg 0x00ff); chksum temp; // 校验和累加起始 fnSpi2_WriteByte(temp); // 发送寄存器高4位 连续读16点指令 0x0CBRUST批量读 temp ((u8)(wReg 4)) 0xf0; temp temp | 0x0c ; chksum temp; fnSpi2_WriteByte(temp); // 循环读取16个采样点每个点3字节24bit for(j0;j 16 ; j) { // 每个采样点高字节→中字节→低字节 for(i 3; i 0;i--) { pBuf[(ij*3-1)] fnSpi2_ReadByte(); chksum pBuf[i-1]; // 累加校验和 } } // 读取芯片返回校验字节校验不匹配则放弃数据 if(fnSpi2_ReadByte()!chksum) { DISABLE_RN8302(); } // 拉高片选释放SPI总线 PinWrite_ADCS(1); PinWrite_ADCS(1); PinWrite_ADSDO(1); PinWrite_ADSCLK(1); } /** * brief FFT模块全局初始化 * details 清空采样缓存、FFT结果、复位状态机所有标志 */ void fnDl645Fft_init(void) { Harmonictemp.Channel0x00; Harmonictemp.DataCount0; Harmonictemp.StarFlag0; // 状态机置空闲 Harmonictemp.ADSPIBusy 0; // 清空瞬时采样缓存结构体 memset(Harmonictemp,0,sizeof(sADE_Instantaneous_TypeDef)); // 清空FFT谐波结果存储结构体 memset(HarmonicData,0,sizeof(sDl645FftHarmonicData_TypeDef)); } /** * brief 原始24bit补码采样数据转换为[-1,1]浮点 * param piRetValue 输出浮点数组长度NUM_FFT128 * details RN8302输出24位有符号补码范围 -2^23 ~ 2^23-1归一化浮点用于FFT运算 */ void SampleDataModifyF(float *piRetValue) { u8 i; u32 Temp; // 存放3字节拼接后的24bit无符号临时值 for(i0; iNUM_FFT; i) { piRetValue[i] 0; Temp 0; // 从384字节缓存取出当前采样点3字节拼接至Temp低24位 memcpy((u8 *)Temp , (Harmonictemp.InstantaneousData[i*3]) , 3); // 判断符号位bit231 负数补码 if(Temp 0x800000) { // 24bit补码转原码按位取反1 Temp (0xffffff - Temp) 1; // 除以2^238388608 归一化 [-1,0) piRetValue[i] -(float)Temp/8388608; } else { // 正数直接归一化 [0,1] piRetValue[i] (float)Temp/8388608; } } } /** * brief 128点基2 DIT时分抽取FFT核心函数 * param dataR 输入时域浮点数组 / 输出频谱实部 * param dataI 输入虚部(初始全0) / 输出频谱虚部 * note 两大阶段1.倒位序重排 2.7级蝶形三层循环运算原位运算仅双数组 */ void FFT(float dataR[NUM_FFT],float dataI[NUM_FFT]) { int x0,x1,x2,x3,x4,x5,x6,xx; int i,j,k,b,p,L; float TR,TI,temp; // 蝶形运算临时缓存防止原位覆盖丢失原值 // 阶段17位二进制倒位序重排 // N1282^7下标0~127使用7位二进制b6 b5 b4 b3 b2 b1 b0 // 倒序规则高低位翻转 b0 b1 b2 b3 b4 b5 b6 生成新下标xx for ( i0;iNUM_FFT;i ) { x0i0x01; // 取出bit0 x1(i/2)0x01;// 取出bit1 x2(i/4)0x01;// bit2 x3(i/8)0x01;// bit3 x4(i/16)0x01;//bit4 x5(i/32)0x01;//bit5 x6(i/64)0x01;//bit6 // 翻转后组合新下标 xx b0*64 b1*32 b2*16 b3*8 b4*4 b5*2 b6*1 xxx0*64x1*32x2*16x3*8x4*4x5*2x6; dataI[xx]dataR[i]; // 临时存入虚部数组做中转 } // 倒序完成数据拷贝回实部虚部清零(输入纯实数信号无虚分量) for ( i0;iNUM_FFT;i ) { dataR[i]dataI[i]; dataI[i]0; } // 阶段2三层循环蝶形运算 7级 // 外层循环L级数循环 L1~7 共7级蝶形 for ( L1;L7;L ) { // 计算当前级蝶形两点间距 b2^(L-1) b1; iL-1; while ( i0 ) { bb*2; i--; } // 中层循环j遍历本级所有旋转因子W_N^p共b2^(L-1)种 for ( j0;jb-1;j ) { // 计算旋转因子查表索引 p j * 2^(7-L) p1; i7-L; while ( i0 ) { pp*2; i--; } pp*j; // 内层循环k遍历所有同旋转因子的蝶形单元步长2*b for ( kj;k128;kk2*b ) { // 【关键】缓存原始AdataR[k]原位运算会覆盖必须临时保存 TRdataR[k]; TIdataI[k]; tempdataR[kb]; // 缓存B点实部原始值 // 蝶形公式1、2计算上支路输出 Y A W*B dataR[k] dataR[k] dataR[kb]*cos_tab[p] dataI[kb]*sin_tab[p]; dataI[k] dataI[k] - dataR[kb]*sin_tab[p] dataI[kb]*cos_tab[p]; // 蝶形公式3、4计算下支路输出 Y- A - W*B使用缓存TR/TI dataR[kb] TR - dataR[kb]*cos_tab[p] - dataI[kb]*sin_tab[p]; dataI[kb] TI temp*sin_tab[p] - dataI[kb]*cos_tab[p]; } } } } /** * brief 预生成旋转因子sin/cos查表数组 * details 嵌入式避免实时调用sin/cos浮点函数运算时直接查表提速 * sin_tab[p] sin(2πp/N) ; cos_tab[p] cos(2πp/N) */ void InitForFFT() { int i; for ( i0;iNUM_FFT;i ) { sin_tab[i]sin(PI*2*i/NUM_FFT); cos_tab[i]cos(PI*2*i/NUM_FFT); } } /** * brief 幅值阈值判断信号过小跳过FFT计算节省算力 * param Channel 通道号0~5 * return 1信号有效执行FFT0幅值过低跳过 */ u8 fnDl645Fft_DataCmp(u8 Channel) { u8 Btemp1; // 默认允许计算 switch(Channel) { case 0x00: // A相电压 if(Dl645Bkgrd.PriPara.WireMode 1) { if(Dl645Front.PubData.U[0]0x0070) Btemp0; } else { if(Dl645Front.PubData.U[0]0x1540) Btemp0; } break; case 0x01: // B相电压 if(Dl645Bkgrd.PriPara.WireMode 1) { if(Dl645Front.PubData.U[1]0x0070) Btemp0; } else { if(Dl645Front.PubData.U[1]0x1540) Btemp0; } break; case 0x02: // C相电压 if(Dl645Bkgrd.PriPara.WireMode 1) { if(Dl645Front.PubData.U[2]0x0070) Btemp0; } else { if(Dl645Front.PubData.U[2]0x1540) Btemp0; } break; case 0x03: // A相电流 if((Dl645Front.PubData.I[0] 0x7fffffff) 0x750) Btemp0; break; case 0x04: // B相电流 if((Dl645Front.PubData.I[1] 0x7fffffff) 0x750) Btemp0; break; case 0x05: // C相电流 if((Dl645Front.PubData.I[2] 0x7fffffff) 0x750) Btemp0; break; } return(Btemp); } /** * brief 谐波分析主状态机系统主循环周期调用 * StarFlag 三状态流转0空闲 → 1等待采样缓存 → 2FFT谐波计算 * 每15s、45s整点自动启动一轮6通道谐波检测 */ void fnDl645Fft_Harmonic_Exec(void) { float ftemp_N,ftemp_P; float fBase; // 基波(50Hz)幅值 s16 i; // 状态0空闲定时触发启动采样 if((Harmonictemp.StarFlag0) (SysClock.Second 0x15 || SysClock.Second 0x45)) { Harmonictemp.StarFlag 1; // 切换至等待采样状态 memset(HarmonicData,0,sizeof(sDl645FftHarmonicData_TypeDef)); // 清空上次谐波结果 Harmonictemp.TimeOutStampSysStamp?SysStamp:(SysStamp-1); // 记录启动时间戳 // SPI写寄存器配置RN8302启动128点瞬时采样缓存 Dl645RN8302DataComm.ucTemp8 0xe5; fnRN8302_Write(0x0180,Dl645RN8302DataComm.ucTempBuf,1); Dl645RN8302DataComm.ucTemp8 0xa0; // 开启采样写入缓存 fnRN8302_Write(0x0163,Dl645RN8302DataComm.ucTempBuf,1); Dl645RN8302DataComm.ucTemp8 0xDC; // 写保护锁存 fnRN8302_Write(0x0180,Dl645RN8302DataComm.ucTempBuf,1); Harmonictemp.Channel0x00; // 从A相电压开始计算 Harmonictemp.ReadAdress 0x0200; // RN8302采样缓存起始地址 InitForFFT(); // 预计算sin/cos旋转因子表 return; } // 状态1延时25ms等待芯片缓存128点采样 if(Harmonictemp.StarFlag 0x01) { // 未到25ms直接退出等待下一周期 if(fnStamp_Through(Harmonictemp.TimeOutStamp)25) return; Harmonictemp.StarFlag 0x02; // 切换至FFT计算状态 // 读取采样缓存状态标志缓存未就绪直接放弃本轮分析 fnRN8302_Read( 0x0163 , Dl645RN8302DataComm.ucTempBuf , 1 ); if(Dl645RN8302DataComm.ucTempBuf[0] 0x30) return; // 读取三相电压、三相电流硬件基波有效值 for(i0;i7;i) { Dl645FrontTmp.UI[i] 0 ; fnRN8302_Read( 0x0058i , (u8 *)Dl645FrontTmp.UI[i] , 4 ) ; } // 电压有效值缩放存入FU for(i0;i3;i) { HarmonicData.FU[i] ((s16)(Dl645FrontTmp.UI[i]/(10*(Dl645FirmPara.VRmsConst)))); } // 电流有效值缩放存入FI for(i0;i3;i) { HarmonicData.FI[i] (s32)(Dl645FrontTmp.UI[i3]/(Dl645FirmPara.IRmsConst)); } return; } // 状态2读取采样 → FFT运算 → 谐波含有率/THD计算 if(Harmonictemp.StarFlag 0x02) { // SPI分8次读取128点原始采样(每次16点)存入384字节缓存 for(i0;i8;i) { fnRN8302_ReadBuf(((Harmonictemp.ReadAdress)i*16),(u8 *)Harmonictemp.InstantaneousData[i*48]); } // 24bit补码采样转[-1,1]浮点时域数组 SampleDataModifyF(HarmonicData.iRealArray[0]); // FFT虚部初始清零 for(i0;iNUM_FFT;i) HarmonicData.iMageArray[i]0; // 执行128点DIT-FFT输出频谱实虚部 FFT(HarmonicData.iRealArray, HarmonicData.iMageArray); // 计算基波(50Hz)幅值k1对应1次基波 ftemp_N (float)(((float)HarmonicData.iRealArray[1] * HarmonicData.iRealArray[1]) ((float)HarmonicData.iMageArray[1] * HarmonicData.iMageArray[1])); ftemp_N sqrt(ftemp_N); fBase ftemp_N; // 基波幅值作为分母计算谐波占比 ftemp_P 0; // 累加各次谐波平方用于THD // 循环计算1~50次谐波含有率 for(i0;i50;i) { // ki2 对应i1次谐波k1基波k2二次谐波 ftemp_N (float)(((float)HarmonicData.iRealArray[i2] * HarmonicData.iRealArray[i2]) ((float)HarmonicData.iMageArray[i2] * HarmonicData.iMageArray[i2])); ftemp_N sqrt(ftemp_N); // 乘以硬件频率补偿系数放大1e6存储小数百分比 ftemp_N ( ftemp_N * 1000000 * Fftcoefficient[i]) / fBase; // 谐波平方累加THD使用平方和开根号 ftemp_P (ftemp_N * ftemp_N); // 区分电压/电流通道存储谐波百分比(BCD格式) if(Harmonictemp.Channel0x03) { HarmonicData.HarmonicpercentU[Harmonictemp.Channel][i1] fnHexToBcd_u32((u32)(ftemp_N)) ; } else { HarmonicData.HarmonicpercentI[Harmonictemp.Channel-3][i1] fnHexToBcd_u32((u32)ftemp_N) ; } } // 计算总谐波畸变率THD ftemp_P sqrt(ftemp_P); if(Harmonictemp.Channel0x03) { HarmonicData.HarmonicpercentU[Harmonictemp.Channel][0] fnHexToBcd_u32((u32)(ftemp_P)); } else { HarmonicData.HarmonicpercentI[Harmonictemp.Channel-3][0] fnHexToBcd_u32((u32)ftemp_P) ; } // 切换下一采集通道地址偏移128点 Harmonictemp.Channel; Harmonictemp.ReadAdress128; // 跳过幅值过低无信号通道 for(iHarmonictemp.Channel;i6;i) { if(fnDl645Fft_DataCmp(Harmonictemp.Channel)) break; Harmonictemp.Channel; Harmonictemp.ReadAdress128; } // 6通道全部计算完成状态机切回空闲0 if(Harmonictemp.Channel 5 ) { Harmonictemp.StarFlag 0; return; } } // 超时保护超过8000ms未完成强制复位状态机 if(fnStamp_Through(Harmonictemp.TimeOutStamp)8000) { Harmonictemp.StarFlag 0; } }三、核心逻辑分层总结便于阅读理解1. 数据链路RN8302 硬件 24bit ADC 采样 → SPI 批量读取fnRN8302_ReadBuf→ 3 字节缓存 →SampleDataModifyF补码转归一化浮点 → FFT 蝶形运算 → 各次谐波幅值 / THD 计算2. FFT 核心两段逻辑倒位序DIT 时域抽取算法强制前置步骤二进制高低位翻转重排输入数组7 级三层循环蝶形外层 L级数 1~7中层 j切换不同旋转因子 W内层 k批量蝶形加减运算TR/TI/temp缓存防止原位覆盖出错。3. 状态机分工无阻塞 MCUStarFlag0整点触发配置芯片开启采样StarFlag1延时 25ms 等待硬件缓存 128 点StarFlag2读取采样→FFT→谐波计算自动切换 6 路通道超时兜底 8000ms 强制复位防止卡死。4. 工程适配重点Fftcoefficient补偿芯片 ADC 高频衰减解决高次谐波计量误差阈值判断fnDl645Fft_DataCmp无信号跳过 FFT 节省浮点算力查表InitForFFT嵌入式 MCU 禁用实时三角函数大幅提速原位 FFT仅 2 个 float 数组极低 RAM 占用适配单片机。