dxdb.c

在CCS开发工具下交流采样并实现FFT算法

#include <math.h>
#define N  1024
#pragma DATA_SECTION(Uads,"Uad_Sample");
float Uads[8*N];
#pragma DATA_SECTION(Iads,"Iad_Sample");
float Iads[8*N];
#pragma DATA_SECTION(Har,"Har_analyse");
float Har[N];
#pragma DATA_SECTION(W,"W_xishu");
float W[N];
#define pi 3.141592653589793238
float Urms=0;
float Irms=0;
float P=0,Q=0,S=0;
float PF=0,PH=0,F=0;

void fftr2_dit(float * x,float * w,short n);
void bit_rev(float * x,short n);

void main()
{
  short i=0;
  float e=2.0*pi/N;
  float Usum=0;
  float Isum=0;
  float Psum=0;
  float Qsum=0;
  for(i=0;i<4*N;i++)
  {
     
     Uads[2*i]=1.4142135623730950488*cosf(2*pi*i/N);
     Iads[2*i]=1.4142135623730950488*cosf(2*pi*i/N-pi/3);
     Uads[2*i+1]=0.0;
     Iads[2*i+1]=0.0;
  }
  for(i=0;i<N/2;i++)
  {
     W[2*i]=cosf(i*e);
     W[2*i+1]=sinf(i*e);
  }
  //电压有效值
  for (i=0;i<N;i++)
  {
     Usum =Usum+Uads[2*i]*Uads[2*i];
  }
  Urms=sqrt(Usum/N);
  //电流有效值
  for (i=0;i<N;i++)
  {
    Isum =Isum+Iads[2*i]*Iads[2*i];
  }
  Irms=sqrt(Isum/N);
  //有功功率
  for (i=0;i<N;i++)
  {
    Psum=Psum+Uads[2*i]*Iads[2*i];
  }
  P=Psum/N;
  //无功功率
  for (i=0;i<N;i++)
  {
    Qsum=Qsum+Uads[2*i]*Iads[2*i+2*N/4];
  }
  Q=Qsum/N;
  //视在功功率
  S=sqrtf(P*P+Q*Q);
  //功率因数
  PF=P/S;
  //相位差
  PH=acosf(PF);
  //谐波分析
 bit_rev(W,(N/2));
 fftr2_dit(Uads,W,N);
 bit_rev(Uads,N);
 for(i=0;i<N;i++)
 {
    Har[i]=sqrtf(Uads[2*i]*Uads[2*i]+Uads[2*i+1]*Uads[2*i+1])/512;
 }
}

//位反转
void bit_rev(float * x,short n)
{
    short i,j,k;
    float rtemp,itemp;
    j=0;
    for(i=1;i<(n-1);i++)
    {
      k = n>>1;
      while(k<=j)
      {
        j-=k;
        k>>=1;
      }
      j+=k;
      if(i<j)
      {
        rtemp=x[j*2];
        x[j*2]=x[i*2];
        x[i*2]=rtemp;
        itemp=x[j*2+1];
        x[j*2+1]=x[i*2+1];
        x[i*2+1]=itemp;
      }
    }
}
//时分基2 FFT
void fftr2_dit(float * x,float * w,short n)
{
    short n2,ie,ia,i,j,k,m;
    float rtemp=0,itemp=0,c=0,s=0;
    n2=n;
    ie=1;
    for(k=n;k>1;k>>=1)
    {
       n2>>=1;
       ia =0;
       for(j=0;j<ie;j++)
       {
         c = w[2*j];
         s = w[2*j+1];
         for(i=0;i<n2;i++)
         {
           m = ia + n2;
           rtemp= c*x[2*m]+s*x[2*m+1];
           itemp= c*x[2*m+1]-s*x[2*m];
           x[2*m]=x[2*ia]-rtemp;
           x[2*m+1]=x[2*ia+1]-itemp;
           x[2*ia]=x[2*ia]+rtemp;
           x[2*ia+1]=x[2*ia+1]+itemp;
           ia++;
         }
         ia +=n2;        
       }
       ie<<=1;
    }
}