📄 harmonic.c
字号:
#include "includes.h"
#include "basetype.h"
#include "DrvCfg.h"
#include "cs5463.h"
#include "HRJ803.h"
#include "driver.h"
#include "Pulse.h"
#include "db.h"
#include "UartDrv.h"
#include "math.h"
//定义全局的基波复数形式数组
DWORD dwLastTicks; //为3秒钟统计一次数据，
BYTE byCurXieboTime; //当前3秒钟内统计的次数，一般要求3秒内至少6次
//BYTE bySaveDayHarm; //保存日统计数据的标志
SDateTime sLastTimeHarm; //为保存日统计数据，上次统计的时间
#define NUM_FFT 128 //fft_code_tables.h中的NUM_FFT 要也相应的修改，当NUM_FFT改变时
void Compute_Harmonic(INT16 iRealArray[], INT16 iMageArray[],BYTE byIndex ,SRestrict *psRestrict);
void SampleDataModify(DWORD* pWaveCode, DWORD dwIndex,DWORD dwFreq,INT16 *piRetValue);
typedef struct
{
WORD wUa1; //基波电压有效值 系数100 单位V
WORD wIa1; //基波电流有效值 系数1000 单位A
WORD wUb1;
WORD wIb1;
WORD wUc1;
WORD wIc1;
//*************电压电流的谐波有效值 电压系数100 单位V 电流系数1000单位A
WORD wHRUa_rms[XIEBO_NUM]; //[0]2次 [1]3次 ...
WORD wHRIa_rms[XIEBO_NUM];
WORD wHRUb_rms[XIEBO_NUM];
WORD wHRIb_rms[XIEBO_NUM];
WORD wHRUc_rms[XIEBO_NUM];
WORD wHRIc_rms[XIEBO_NUM];
//***********以下为1分钟内的 3秒钟数据 为单位 所得的95%概率大值
/* WORD wTHDua; //A相电压总畸变含有率 系数100 单位%
WORD wHRUa[XIEBO_NUM]; //A 相电压谐波含有率 系数100 单位%
WORD wTHDIa; //A 相电流总畸变电流 系数1000 单位A
WORD wHIa; //A相总畸变电流 系数1000 单位A
WORD wTHDub; //B相电压总畸变含有率
WORD wHRUb[XIEBO_NUM]; //B 相电压谐波含有率
WORD wTHDIb; //B 相电流总畸变电流
WORD wHIb; //B相总畸变电流
WORD wTHDuc; //C相电压总畸变含有率
WORD wHRUc[XIEBO_NUM]; //C 相电压谐波含有率
WORD wTHDIc; //C 相电流总畸变电流
WORD wHIc; //C相总畸变电流
*/
}SPowerQualitySum; //电能质量 三秒钟采集的谐波数据累加和,一般为6次
//static SPowerQualitySum sPowerQualitySum;
//*****************拉格朗日插值**************************
typedef struct
{
float x;
float y;
}SData;//变量x和函数值y的结构
/*特例：拉格朗日三点插值又称为抛物线插值法
y= y1*(x-x2)(x-x3)/(x1-x2)(x1-x3) + y2*(x-x1)(x-x3)/(x2-x1)(x2-x3) +y3*(x-x1)(x-x2)/(x3-x1(x3-x2)
其中，(x1,y1)、(x2,y2)、(x3,y3)为三个样点坐标,(x,y)为要获得的点
*/
//count 为插值的点数
float lagrange(float x, int count,SData d[])
{
float y=0.0;
float p=1.0;//初始化p
int k = 0;
int j =0;
for( k=0; k<count; k++)//这儿默认为count－1次插值
{
p = 1.0;
for(j=0; j<count; j++)
{//计算p的值
if(k==j)continue;//判断是否为同一个数
p=p*(x-d[j].x)/(d[k].x-d[j].x);
}
y=y+p*d[k].y;//求和
}
return y;//返回y的值
}
//*********采样数据修正****************
//这里假设AD每次采样完成所需的时间是一致的、理想的，
void SampleDataModify(DWORD* pWaveCode, DWORD dwIndex,DWORD dwFreq,INT16 *piRetValue)
{
int i;
int j;
int iCurIndex;
int ThreePointx[3];
float fADSampleInterval; //AD采样的时间间隔 若为4000Hz则为0.25ms
float fFFtNeedInterval; //当前频率dwFreq， 所对应的FFT间隔。
float fCurTime;
float ftemp;
INT16 itemp;
//例如，当前频率50HZ,NUM_FFT＝128（2个周波），AD采样200点。最好的采样率应为64*50=3200Hz,则fCurInterval＝1/3200 = 0.3125ms
//这里采用拉格朗日三点插值，把200个采样点，插值为128个点，则插值后的采样频率为3200Hz。然后再进行FFT运算，达到了频率的跟踪
SData sData[3];
INT16 iValue[MAX_WAVE_NUM];
pWaveCode +=dwIndex;
for(i=0; i<MAX_WAVE_NUM; i++) //A相视在功率 有功功率 无功功率 功率因数
{
iValue[i] = (INT16)((*pWaveCode & 0x00ffff00)>>8);
pWaveCode += sizeof(SWaveCodeItem)/sizeof(DWORD);
}//for
fADSampleInterval = (float)1.0/(float)4000.0; //0.25ms
fFFtNeedInterval = (float)1.0 / ((float)3200.0 * (float)dwFreq /(float)0x19999); //因50Hz时插值后的采样率应为3200Hz，那么若当前频率为55Hz则插值后的采样率 ＝ 3200*55/50; 而50Hz对应AD频率寄存器的0x19999。
piRetValue[0] = iValue[0]; //第一个点取 iValue的第二个点，iValue的第一个点为可能的第一次插值用，也可能用不到，依据当前的实时频率，见下面的语句
fCurTime = fFFtNeedInterval;
j=1;
for(i=1; i<NUM_FFT; i++) //fft的点数
{
iCurIndex = (int)(fCurTime / fADSampleInterval); //因第一个取的是iValue的第二个点，所以以后的要都后推
if(fCurTime < ( (float)iCurIndex*fADSampleInterval + fADSampleInterval/(float)2.0)) //本次FFT插值的时间点，位于当前iCurIndexAD采样时间点与iCurIndex＋1采样点的前半部的时间
iCurIndex--; //
sData[0].x = (float)((float)iCurIndex * fADSampleInterval);
sData[1].x = (float)(iCurIndex + 1) * (fADSampleInterval);
sData[2].x = (float)(iCurIndex + 2)* (fADSampleInterval);
sData[0].y = (float)iValue[iCurIndex];
sData[1].y = (float)iValue[iCurIndex+1];
sData[2].y = (float)iValue[iCurIndex+2];
ftemp = lagrange(fCurTime, 3,sData);
itemp =ftemp;
piRetValue[j++] = itemp;
fCurTime +=fFFtNeedInterval;
}
}
const BYTE cIndex_Coef[6] = {0,3,1,4,2,5};
//********计算基波、总的畸变率、各次谐波含有率,fft的调用
void Compute_Harmonic(INT16 iRealArray[], INT16 iMageArray[],BYTE byIndex ,SRestrict *psRestrict)
{
BYTE i = 0;
float ftemp_N;
float fBase; //基波的平方
float fSum =0;
WORD *pwHR;
//*************FFT*****************
//WindowCalc( iRealArray, 0);
Bit_Reverse(iRealArray);
Int_FFT(iRealArray, iMageArray);
//************ 计算谐波有效值 计算各次谐波含有率， 总畸变率
fBase = (float)(((float)iRealArray[2] *(float)iRealArray[2]) + ((float)iMageArray[2] * (float)iMageArray[2]));
ftemp_N = sqrt(fBase);
if( (byIndex & 0x01 ) ==0)
ftemp_N *= 25000.0*2*1.414213 ; //250*100
else ftemp_N *=6000.0 *2* 1.414213; //6000
ftemp_N /=10066329.6; //((2^24)*0.6=10066329.6)=0x999999
ftemp_N *=256; //因在进行FFT运算前，把AD的数据右移了8为，这里要左移回来，见
ftemp_N *=(float)psRestrict->iValueCoef[cIndex_Coef[byIndex]] /1000.0; //通过系数的修正
pwHR = (WORD*)&g_sPowerQuality.wUa1 + byIndex;
*pwHR = ftemp_N;
pwHR = (WORD*)&g_sPowerQuality.wHRUa[0] + (byIndex* XIEBO_NUM );
fSum = 0;
for(i=0; i<XIEBO_NUM; i++)
{
ftemp_N = ((long)iRealArray[4+i*2] *(long)iRealArray[4+i*2]) + ((long)iMageArray[4+i*2] * (long)iMageArray[4+i*2]);
fSum += ftemp_N;
ftemp_N = ftemp_N / fBase;
ftemp_N = sqrt(ftemp_N)* (float)psRestrict->iHRCoef[cIndex_Coef[byIndex]] /1000.0 *10000.0 ;
pwHR[i] = ftemp_N;
}
pwHR = (WORD*)&g_sPowerQuality.wTHDua + byIndex;
ftemp_N = 10000.0 * sqrt(fSum/fBase)*(float)psRestrict->iHRCoef[cIndex_Coef[byIndex]] /1000.0 ;
*pwHR = ftemp_N ;
}
BYTE byCurrentError_type = 0; //电流异常标志：
BYTE byUpdate_PhaseFlag = 0;
//**********计算相角 及电流开路和反极性的判断*******************
void Compute_Phase(INT16 iRealArray[],INT16 iMageArray[])
{
float Setx[6]; //基波的x轴坐标
float Sety[6]; //基波的y轴坐标
float UxOfA;
float UyOfA;
float RMS_UA1;
long ltemp;
BYTE i;
BYTE byError[3];
INT16 *ptr;
long lRealI0;
long lImageI0;
byUpdate_PhaseFlag = 1;
RMS_UA1 = (float)( (long)iRealArray[0] * (long)iRealArray[0] + (long)iMageArray[0]* (long)iMageArray[0] );
RMS_UA1 = sqrt(RMS_UA1);
if(RMS_UA1 ==0) UxOfA = UyOfA = 0; //防止 被除数为0
else
{
UxOfA = (float)iMageArray[0] / (float)RMS_UA1; //COS
UyOfA = (float)iRealArray[0] / (float)RMS_UA1; //SIN
}//
for(i=1; i<6; i++) //使得A相的基波电压位于Y轴，其他的电压电流的基波相应的旋转，
{
Setx[i] = ((float)iRealArray[i] * (float)UxOfA - (float)iMageArray[i] *(float)UyOfA);
Sety[i] = ((float)iRealArray[i] * (float)UyOfA + (float)iMageArray[i] * (float)UxOfA);
}
ptr =(INT16*)&g_sPowerQuality.iUaPhase;
*ptr =0;//以A相电压的为基准
ptr ++;
for(i=1; i<6; i++,ptr++) //计算 Ia1 Ub1 Ib1 Ic1的基波角度
{
// *ptr = 10 * atan2((double)Sety[i],(double)Setx[i])*180.0/3.1415926;
if(Setx[i]==0 && Sety[i]>=0)
{
*ptr =0;
}
else if(Setx[i]==0 && Sety[i]<0)
{
*ptr =1800; //180度
}
else if(Setx[i]>=0 && Sety[i]==0)
{
*ptr =2700;//270度
}
else if(Setx[i]<0 && Sety[i]==0)
{
*ptr =900;//90度
}
else if(Setx[i] >0 && Sety[i] >0) //位于第1象限
{
*ptr = (INT16)2700 + 10 * atan(Sety[i]/Setx[i])*180.0/3.1415926 ;
}
else if(Setx[i] <0 && Sety[i] >0) //位于第2象限
{
*ptr = (INT16)900 - (INT16)(10 * atan((double)Sety[i]/(double)-Setx[i])*180.0/3.1415926) ;
}
else if(Setx[i] <0 && Sety[i] <0) //位于第3象限
{
*ptr = (INT16)900 + (INT16)(10*atan((double)Sety[i]/(double)Setx[i])*180.0/3.1415926);
}
else if(Setx[i] >0 && Sety[i] <0) //位于第4象限
{
*ptr = (INT16)(2700 + (INT16)(10* atan((double)Sety[i]/(double)Setx[i])*180.0/3.1415926));
}
}
i =0;
//某相电压小于失流的电压条件(电压≥30%Ub)时,不进行逆相序判定,保持原状态
// Ub 相位不处于(210～270 度) 或 Uc 相位不处于(90～150 度)
if( g_sRealVal.wUa < 2200*0.3 || g_sRealVal.wUb <2200*0.3 || g_sRealVal.wUc < 2200*0.3)
{
i = 1;
}
if(i !=1) //当电压不是逆相序时才判断电流是否反极性
{
i = 0;
if(g_sPowerQuality.wIa1 < 150 ) //0.05A = 5*3%
{
byError[0] =1;
i++;
}
if(g_sPowerQuality.wIb1 < 150 )
{
byError[1] = 1;
i++;
}
if( g_sPowerQuality.wIc1 <150)
{
byError[2] =1;
i++;
}
//某相电流小于失流的起始电流条件时 或某相电压小于失流的电压条件时,则不判断是否反极性
if(i==0)//
{
lRealI0 = (long)iRealArray[1] + (long)iRealArray[3] + (long)iRealArray[5];
lImageI0 = (long)iMageArray[1] + (long)iMageArray[3] + (long)iMageArray[5];
//判断是否时任意二相电流反向
// lRealI0 = iRealArray[1] + iRealArray[3] + iRealArray[5];
// lImageI0 = iMageArray[1] + iMageArray[3] + iMageArray[5];
if( lRealI0 < 300 && lRealI0 >-300)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -