dft的51程序.c

DFT算法的Keil C51实现

#include <reg52.H>
#include <MATH.H>
#include <absacc.h>


float xdata tab_cos[128];
float xdata tab_sin[128];
float xdata tab_wave[128];
float xdata tab_dftR[128];
float xdata tab_dftI[128];
float xdata tab_adftR[128];
float xdata tab_adftI[128];
float xdata tab_mod[128];
float xdata tab_phase[128];

void fft(float dataR[],float dataI[]);
void afft(float dataR[],float dataI[]);
void get_mod_and_pha(float dataR_mod[],float dataI_pha[]);

main()
{  unsigned char i;
   float xdata *ptrd;
   float xdata *ptrs;
   float xdata tmp,tmp1,tmp2;
   //tmp=pow(3,4);
   tmp=6.28319/256;
   tmp=tmp*tmp;
   ptrd=0x0fff;
   ptrs=&tmp2;
   *ptrd=tmp;
   *ptrs=*ptrd;
   tmp=6.28319/128;
   for(i=0;i<128;i++)
   {
      //tab_wave[i]=(float)i/128;
      tab_wave[i]=cos(tmp*i);  //设置用于测试的波
      tab_cos[i]=cos(tmp*i);   //预先设置余弦表
      tab_sin[i]=sin(tmp*i);   //设置正弦表

      tab_dftR[i]=tab_wave[i];//读入波形
   }

   fft(tab_dftR,tab_dftI);
   for(i=0;i<128;i++)
   {
      tab_adftR[i]=tab_dftR[i];
      tab_adftI[i]=tab_dftI[i];
   }
   afft(tab_adftR,tab_adftI);
   for(i=0;i<128;i++)
   {
      tab_mod[i]=tab_dftR[i];
      tab_phase[i]=tab_dftI[i];
   }
   get_mod_and_pha(tab_mod,tab_phase);
   while(1){}
}

void fft(float dataR[],float dataI[])         //注：外部必须设置正确的正、余弦表，返回是
{   unsigned char x0,x1,x2,x3,x4,x5,x6,xx;
    unsigned char i,L,j,k,b,p;
    float TR,TI,temp;
    /********** following code invert sequence ************/
    for(i=0;i<128;i++) //蝶形算法前需要对数列顺序进行重排
    {   x0=x1=x2=x3=x4=x5=x6=0;
        x0=i&0x01;
        x1=(i>>1)&0x01;
        x2=(i>>2)&0x01;
        x3=(i>>3)&0x01;
        x4=(i>>4)&0x01;
        x5=(i>>5)&0x01;
        x6=(i>>6)&0x01;
        xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;
        dataI[xx]=dataR[i];
    }
    for(i=0;i<128;i++)
    { dataR[i]=dataI[i]; dataI[i]=0; }
    /************** following code FFT *******************/
    for(L=1;L<=7;L++)
    { /* for(1) */
        b=1; i=L-1;
        while(i>0)
        {b=b*2;i--;}               /* b= 2^(L-1) */
        for(j=0;j<=b-1;j++)        /* for (2) */
        {   p=1; i=7-L;
            while(i>0)             /* p=pow(2,7-L)*j; */
            {p=p*2;i--;}
            p=p*j;
            for(k=j;k<128;k=k+2*b) /* for (3) 蝶形算法*/
            {    TR=dataR[k]; TI=dataI[k]; temp=dataR[k+b];
                 dataR[k]=TR+dataR[k+b]*tab_cos[p]+dataI[k+b]*tab_sin[p];
                 dataI[k]=TI-dataR[k+b]*tab_sin[p]+dataI[k+b]*tab_cos[p];
                 dataR[k+b]=TR-dataR[k+b]*tab_cos[p]-dataI[k+b]*tab_sin[p];
                 dataI[k+b]=TI+temp*tab_sin[p]-dataI[k+b]*tab_cos[p];
            } /* END for (3) */
        }     /* END for (2) */
    }         /* END for (1) */
}             /* END FFT */

void afft(float dataR[],float dataI[])
{   unsigned char x0,x1,x2,x3,x4,x5,x6,xx;
    unsigned char i,L,j,k,b,p;
    float TR,TI;
    for(L=7;L>0;L--)
    { /* for(1) */
        b=1; i=L-1;
        while(i>0)
        {b=b*2;i--;}               /* b= 2^(L-1) */
        for(j=0;j<=b-1;j++)        /* for (2) */
        {   p=1; i=7-L;
            while(i>0)             /* p=pow(2,7-L)*j; */
            {p=p*2;i--;}
            p=p*j;
            for(k=j;k<128;k=k+2*b) /* for (3) 反蝶形算法**/
            {    TR=dataR[k]-dataR[k+b]; TI=dataI[k]-dataI[k+b];
                 dataR[k]=dataR[k]+dataR[k+b];
                 dataI[k]=dataI[k]+dataI[k+b];
                 dataR[k+b]=TR*tab_cos[p]-TI*tab_sin[p];
                 dataI[k+b]=TR*tab_sin[p]+TI*tab_cos[p];
            } /* END for (3) */
        }     /* END for (2) */
    }         /* END for (1) */
    for(i=0;i<128;i++)
    {
        dataR[i]=dataR[i]/128;
    }
    for(i=0;i<128;i++)   //数列反重排
    {   x0=x1=x2=x3=x4=x5=x6=0;
        x0=i&0x01;
        x1=(i>>1)&0x01;
        x2=(i>>2)&0x01;
        x3=(i>>3)&0x01;
        x4=(i>>4)&0x01;
        x5=(i>>5)&0x01;
        x6=(i>>6)&0x01;
        xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;
        dataI[xx]=dataR[i];
    }
    for(i=0;i<128;i++)
    { dataR[i]=dataI[i]; dataI[i]=0; }
}
void get_mod_and_pha(float dataR_mod[],float dataI_pha[])//fft算法后，对频谱取模和相位角
{  unsigned char i;
   float tmp;
   for(i=0;i<128;i++)
   {  tmp=dataR_mod[i];
      dataR_mod[i]=sqrt(tmp*tmp+dataI_pha[i]*dataI_pha[i]);    //取模值
      dataI_pha[i]=atan2(dataI_pha[i],tmp);                    //取相位
   }
}