harmonic.c

FFT的实现在LPC2378上。首先采用了拉个朗日插值算法对采集的电网波形频率跟踪

		{   
			if(PtrCur[i]  > PtrMax[i] || PtrMax[i]==0)
			{
				PtrMax[i] = PtrCur[i];
				PtrMaxTime[i] = sTime.byMonth;
				PtrMaxTime[i] <<= 8;
				PtrMaxTime[i] |= sTime.byDay;
				PtrMaxTime[i] <<= 8;
				PtrMaxTime[i] |= sTime.byHour;
				PtrMaxTime[i] <<= 8;
				PtrMaxTime[i] |= sTime.byMinute;
			}
		}
        if(j==0)
			PtrCur = (WORD*)&g_sPowerQuality.wHRUb[0];
		else
			PtrCur = (WORD*)&g_sPowerQuality.wHRUc[0];
	}
	OSTimeDly(2);
	//到了一分钟,需要统计次分钟的数据是否越限，根据当前的统计次数 与 95%概率来计算
	if( (OSTimeGet() - dwLastTicks) >= OS_TICKS_PER_SEC * 60)
	{
        dwLastTicks = OSTimeGet(); //
		//谐波越限日统计数据
		 //*********A相谐波谐波总畸变率限值判断
		if(g_sPowerQuality.wTHDua  >psRestrict->wTHDLevelU *10)
		{
			g_sPowerQuality.wTHDua_OverT++;
		}
		//*************A相谐波各次电压含有率限值判断
		for(i=0; i<XIEBO_NUM; i++)
		{
			if(i & 0x01)  //奇次谐波
			{
				if(g_sPowerQuality.wHRUa[i]  >psRestrict->wHDLevelUo *10)
				{
					g_sPowerQuality.wHRUa_OverT[i] ++;
				}
			}
		    else //偶次谐波
		    {
			   if(g_sPowerQuality.wHRUa[i] >psRestrict->wHDLevelUe *10)
			   {
				  g_sPowerQuality.wHRUa_OverT[i] ++;
			   }
		    }
		}
		if((long)g_sPowerQuality.wTHDIa * (long)g_sPowerQuality.wIa1 > (long)100000*(long)psRestrict->iTotalHarmonicContentLevel )
		{
			g_sPowerQuality.wTHDIa_OverT++;
		}
		for(i=0; i<XIEBO_NUM; i++)
		{
            if((long)g_sPowerQuality.wHRIa[i] * (long)g_sPowerQuality.wIa1 > (long)100000*(long)psRestrict->iHarmonicContentLevel[i] )
            g_sPowerQuality.wHIa_OverT[i]++;
		}
	   //***************B相谐波谐波总畸变率限值判断
		if(g_sPowerQuality.wTHDub >psRestrict->wTHDLevelU*10)
		{
			g_sPowerQuality.wTHDub_OverT++;
		}
		//B相谐波各次电压含有率限值判断
		for(i=0; i<XIEBO_NUM; i++)
		{
		   if(i & 0x01)  //奇次谐波
		   {
			  if(g_sPowerQuality.wHRUb[i] >psRestrict->wHDLevelUo *10)
			  {
				g_sPowerQuality.wHRUb_OverT[i] ++;
			  }
		   }
		  else //偶次谐波
		  {
			 if(g_sPowerQuality.wHRUb[i] >psRestrict->wHDLevelUe *10)
			 {
				g_sPowerQuality.wHRUb_OverT[i] ++;
			 }
		  }
		}
		OSTimeDly(2);
		if((long)g_sPowerQuality.wTHDIb * (long)g_sPowerQuality.wIb1 > (long)100000*(long)psRestrict->iTotalHarmonicContentLevel )
		{
			g_sPowerQuality.wTHDIb_OverT++;
		}
		for(i=0; i<XIEBO_NUM; i++)
		{
            if((long)g_sPowerQuality.wHRIb[i] * (long)g_sPowerQuality.wIb1 > (long)100000*(long)psRestrict->iHarmonicContentLevel[i] )
            g_sPowerQuality.wHIb_OverT[i]++;
		}
		//***************C相谐波谐波总畸变率限值判断
		if(g_sPowerQuality.wTHDuc >psRestrict->wTHDLevelU *10)
		{
			g_sPowerQuality.wTHDuc_OverT++;
		}
		//C相谐波各次电压含有率限值判断
		for(i=0; i<XIEBO_NUM; i++)
		{
		   if(i & 0x01)  //奇次谐波
		   {
			 if(g_sPowerQuality.wHRUc[i] >psRestrict->wHDLevelUo *10)
			 {
				g_sPowerQuality.wHRUc_OverT[i] ++;
			 }
		   }
		   else //偶次谐波
		   {
			 if(g_sPowerQuality.wHRUc[i] >psRestrict->wHDLevelUe *10)
			 {
				g_sPowerQuality.wHRUc_OverT[i] ++;
			 }
		   }
		}
		OSTimeDly(1);
		if((long)g_sPowerQuality.wTHDIc * (long)g_sPowerQuality.wIc1 > (long)100000*(long)psRestrict->iTotalHarmonicContentLevel )
		{
			g_sPowerQuality.wTHDIa_OverT++;
		}
		for(i=0; i<XIEBO_NUM; i++)
		{
            if((long)g_sPowerQuality.wHRIc[i] * (long)g_sPowerQuality.wIc1 > (long)100000*(long)psRestrict->iHarmonicContentLevel[i] )
            g_sPowerQuality.wHIa_OverT[i]++;
		}
    }
    GetTime(&sTime);
    sTime4.byMonth = sTime.byMonth;
	sTime4.byDay = sTime.byDay;
    sTime4.byHour = sTime.byHour;
    sTime4.byMinute = sTime.byMinute;
  //  lCurSeconds = GetSecondSinceEpoch(&sTime);
    if(sTime.byDay !=  sLastTimeHarm.byDay)	//判断是否要保存天数据
	{
        SavePQRec(&sTime4, &g_sPowerQuality);
        memcpy(&sLastTimeHarm, &sTime, sizeof(SDateTime));
		memset(&g_sPowerQuality,0,sizeof(SPowerQuality));//所有数据清零为下次准备
	}
}

void InitHarmonic(void)
{
   dwLastTicks = OSTimeGet();
   byCurXieboTime = 0;
   GetTime(&sLastTimeHarm);
  
}
//*************谐波分析计算 统计************
void ProHarmonic(SRestrict *psRestrict)
{
     WORD i;
	 WORD j = 0;
	 BYTE byIndex;
	 long lTemp;
	 float fTemp;
     INT16 iRealValue[NUM_FFT];
     INT16 iMageValue[NUM_FFT];    
	 INT16 iComplex_Real[6]; //基波电压电流的实部，  Ua Ia Ub Ib Uc Ic
	 INT16 iComplex_Image[6];//基波电压电流的虚部  Ua Ia Ub Ib Uc Ic
	  // 以下为线电压的计算
	 INT16 iSampleUa[NUM_FFT/4];
	 INT16 iSampleUb[NUM_FFT/4];
	 INT16 iSampleUc[NUM_FFT/4];
	 INT16 iSampleUab[NUM_FFT/4];
	 INT16 iSampleUbc[NUM_FFT/4];
	 INT16 iSampleUca[NUM_FFT/4];
	 INT16 *ptr;
	 byIndex = 0;
     for(i=0; i<6; i++)
     {
        SampleDataModify((DWORD*)&g_sWaveCode.sWaveCodeItem[0], i, g_sWaveCode.dwFreq, &iRealValue[0]);//采样数据修正
        if(i==0 || i==2 || i==4)  //为计算线电压取得UA UB UC的 波形 
        {
            ptr = &iSampleUa[0] +( i*NUM_FFT/4);
        	for(j=0; j<NUM_FFT/4; j++)
        	{
          	   ptr[j] = iRealValue[j*4] ;
      	 	}
      	}
      	OSTimeDly(1);
        Compute_Harmonic(iRealValue,iMageValue,i,psRestrict);//计算基波、总的畸变率、各次谐波含有率、
		iComplex_Real[byIndex] = iRealValue[2];
		iComplex_Image[byIndex] = iMageValue[2];
		byIndex ++; 
		OSTimeDly(1);
     }
     
	 for(j=0; j<NUM_FFT/4; j++)
	 {
	    iSampleUab[j] = iSampleUa[j] - iSampleUb[j];
	 }
	 for(j=0; j<NUM_FFT/4; j++)
	 {
	    iSampleUbc[j] = iSampleUb[j] - iSampleUc[j];
	 }
	 for(j=0; j<NUM_FFT/4; j++)
	 {
	    iSampleUca[j] = iSampleUc[j] - iSampleUa[j];
	 }
     OSTimeDly(1);
	 //***********计算AB线电压的有效值*********
	 lTemp  = 0;
	 for(j=0; j<NUM_FFT/4; j++)
	 {
         lTemp += (long)iSampleUab[j] * (long)iSampleUab[j];	   
	 }
	 lTemp /= (NUM_FFT/4);
	 fTemp = lTemp;
	 fTemp = sqrt(fTemp) * 2500.0 / 10066329.6;
     g_sRealVal.wUAB = fTemp;

	 //*********计算BC线电压的有效值**********
     lTemp  = 0;
	 for(j=0; j<NUM_FFT/4; j++)
	 {
         lTemp += (long)iSampleUbc[j] * (long)iSampleUbc[j];	   
	 }
	 lTemp /= (NUM_FFT/4);
	 fTemp = lTemp;
	 fTemp = sqrt(fTemp) * 2500.0 / 10066329.6;
     g_sRealVal.wUBC = fTemp;
     OSTimeDly(1);
    //**********计算CA线电压的有效值**********
	 lTemp  = 0;
	 for(j=0; j<NUM_FFT/4; j++)
	 {
         lTemp += (long)iSampleUca[j] * (long)iSampleUca[j];	   
	 }
	 lTemp /= (NUM_FFT/4);
	 fTemp = lTemp;
	 fTemp = sqrt(fTemp) * 2500.0 / 10066329.6;
     g_sRealVal.wUCA = fTemp;

     Compute_Phase(iComplex_Real,iComplex_Image);//计算相角
     Compute_Unbalance(iComplex_Real,iComplex_Image);//计算不平衡度
     
     WriteCom(0,(BYTE*)&g_sPowerQuality.wU1, sizeof(WORD));
			   WriteCom(0,(BYTE*)&g_sPowerQuality.wU2, sizeof(WORD));
			   WriteCom(0,(BYTE*)&g_sPowerQuality.wU0, sizeof(WORD));
			    WriteCom(0,(BYTE*)&g_sPowerQuality.wI1, sizeof(WORD));
			   WriteCom(0,(BYTE*)&g_sPowerQuality.wI2, sizeof(WORD));
			   WriteCom(0,(BYTE*)&g_sPowerQuality.wI0, sizeof(WORD));
	/*		   Print("777");
			    dwTemp = g_sPowerQuality.wIa1;
               WriteCom(0,(BYTE*)&dwTemp,2);
                dwTemp = g_sPowerQuality.wIb1;
               WriteCom(0,(BYTE*)&dwTemp,2);
                dwTemp = g_sPowerQuality.wIc1;
               WriteCom(0,(BYTE*)&dwTemp,2);
			   dwTemp = g_sPowerQuality.wHRIa[0];
               WriteCom(0,(BYTE*)&dwTemp,2);
               dwTemp = g_sPowerQuality.wHRIa[1];
               WriteCom(0,(BYTE*)&dwTemp,2);
               
                Print("444");
                dwTemp = g_sPowerQuality.wUa1;
               WriteCom(0,(BYTE*)&dwTemp,2);
               
               dwTemp = g_sPowerQuality.wUb1;
               WriteCom(0,(BYTE*)&dwTemp,2);
               dwTemp = g_sPowerQuality.wUc1;
               WriteCom(0,(BYTE*)&dwTemp,2);
			   dwTemp = g_sPowerQuality.wHRUa[0];
               WriteCom(0,(BYTE*)&dwTemp,2);
               dwTemp = g_sPowerQuality.wHRUa[1];
               WriteCom(0,(BYTE*)&dwTemp,2);*/
  /*   i = g_sPowerQuality.iUaPhase;
     WriteCom(0,(BYTE*)&i,4);
     i = g_sPowerQuality.iUbPhase;
     WriteCom(0,(BYTE*)&i,4);
     i = g_sPowerQuality.iUcPhase;
     WriteCom(0,(BYTE*)&i,4);
     i = g_sPowerQuality.wU1;
     WriteCom(0,(BYTE*)&i,4);
     i = g_sPowerQuality.wU2;
     WriteCom(0,(BYTE*)&i,4);
     i = g_sPowerQuality.wU0;
     WriteCom(0,(BYTE*)&i,4);*/
}