MSP430 FFT算法实现

2万 · 2012-11-6 22:31

*****************main programe********************/
#include  <msp430x14x.h>

#include    <math.h>
#include    <stdio.h>
#include    <stdlib.h>
#include    "typedef.h "


float    result[8];
struct    compx s[8];
int    Num=8;
const float pp=3.141592653589793;

extern void FFT(struct compx *s,int num);

double mpr[8],mpi[8],mfr[8],mfi[8];
int n,k,l,il;
extern void kbfft();
extern void FFT_dingdian();
main()
{

// ????
int i;
FILE *fin,*fout;
if((fin=fopen( "indata.txt ", "r "))==NULL)
{
      printf( "can 't open infile ");
      exit(0);
};
for(i=0;i <Num;i++)
{
      //fscanf(fin, "%lf ",&s.real);
      s.real=cos(2*3.1415926*(30)*(i)/80.0);//
      s.imag=(double)0;
      mpr=s.real;
      mpi=0;

};
fclose(fin);
//////////////

//FFT(s,Num);

//n=8;k=3;l=0;il=0;
//kbfft(mpr,mpi,n,k,mfr,mfi,l,il);

FFT_dingdian();

//????
if((fout=fopen( "fftresult.txt ", "w "))==NULL)
{
      printf( "can 't open outfile ");
      exit(0);
};

for(i=0;i <Num;i++)
{
      // result=sqrt(pow(s.real,2)+pow(s.imag,2));
      fprintf(fout, "%lf+(%lf)i ",s.real,s.imag);
};

fclose(fout);

}

//---定点---
#include  <msp430x14x.h>
#include "math.h"
#define NN 8
//************************余弦表**********************************
const  float  COS_tab[91]=
               {1.0000,0.9998,0.9993,0.9986,0.9975,0.9962,0.9945,0.9925,0.9903,
               0.9877,0.9848,0.9816,0.9781,0.9744,0.9703,0.9659,0.9614,0.9563,
               0.9511,0.9455,0.9397,0.9336,0.9272,0.9205,0.9135,0.9063,0.8988,
               0.8910,0.8829,0.8746,0.8660,0.8572,0.8480,0.8387,0.8290,0.8192,
               0.8090,0.7986,0.7880,0.7771,0.7760,0.7547,0.7431,0.7314,0.7193,
               0.7071,0.6947,0.6820,0.6691,0.6561,0.6428,0.6293,0.6157,0.6018,
               0.5878,0.5736,0.5592,0.5446,0.5299,0.5150,0.5000,0.4848,0.4695,
               0.4540,0.4384,0.4226,0.4067,0.3907,0.3746,0.3584,0.3420,0.3256,
               0.3090,0.2924,0.2757,0.2588,0.2419,0.2250,0.2079,0.1908,0.1736,
               0.1564,0.1392,0.1219,0.1045,0.0872,0.0698,0.0523,0.0349,0.0174,
               0.0000};

char Re(char M,char k)
{
  char m=1;
  int n;
  while(M>0){m=m*2;M--;}//m==pow(2,M)
  n=360 * k /m;
  return n;

}

void FFT_dingdian( void )
{
signed char flag1=1,flag2=1;
char  i,j,k,R,I;
char p,p1;
//float i5,i6,i7,i8;
char M,L;
int a,b,x=1,y=1,z=1;
int N1;

struct comp
{
   float re;
   float im;
}aa[NN]={{1,0},{0,0},{-1,0},{0,0},{1,0},{0,0},{-1,0},{0,0}};//,{9,0},{10,0},
         //{11,0},{12,0},{13,0},{14,0},{15,0},{16,0}};//

struct comp bb,bb1,bb2;

L=0;
N1=NN;
while(N1>0)
{N1=N1/2;L++;}// L=log(NN)/log(2);
L=L-1;

//输入变址
for(i=0,j=0;i<NN;i++)
{
   if(i>=j){k=NN/2;}
   else
   {
      k=NN/2;
      bb=aa;
      aa=aa[j];
      aa[j]=bb;

   }
   while(k<=j){j=j-k;k=k/2;
               if(k==j)break;}

   j=j+k;
}
  //变址结束

  //由下开始FFT计算

for(M=1;M<=L;M++)
{
      x=M-1;
      a=1;
      while(x>0){a=a*2;x--;}//a==2的M-1次方
   for(k=0;k<a;k++)//a==2的M-1次方,即每一级有2的M-1次方个不同的旋转因子，步进值是1
                                    //---同时也是蝶形结两结点间距
      {
         y=M;
         b=1;
         while(y>0){b=b*2;y--;}//b==2的M次方
         //p=0;
      for(p=k;p<NN-1;p=p+b)//b==2的M次方----相同旋转因子的蝶形结间距
            {


            p1=p+a;//a==2的M-1次方

            R=Re(M,k);//----------------------=360*k/(2^M)--k最大为2^M-1
            I=Re(M,k);//Re(M,k)==Im(M,k)

            if(R>=0&&R<=90){flag1=1;}
            else if(R>90&&R<=180){ R=180-R;flag1=-1;}
            else if(R>180&&R<=270){ R=R-180;flag1=-1;}
            else if(R>270&&R<=360){ R=360-R;flag1=1;}


            if(I>=0&&I<=90){I=90-I;flag2=1;}
            else if(I>90&&I<=180){I=I-90;flag2=1;}
            else if(I>180&&I<=270){I=270-I;flag2=-1;}
            else if(I>270&&I<=360){I=I-270;flag2=-1;}

            bb.re=aa[p].re;
            bb.im=aa[p].im;


            bb1.re=aa[p1].re*(flag1)*COS_tab[R];
            bb1.im=aa[p1].im*(flag2)*COS_tab[I];
            bb2.re=aa[p1].re*(flag2)*COS_tab[I];
            bb2.im=aa[p1].im*(flag1)*COS_tab[R];
         aa[p].re=bb.re+bb1.re  + bb1.im;
         aa[p].im=bb.im+bb2.im - bb2.re ;
         aa[p1].re=bb.re-bb1.re- bb1.im;
         aa[p1].im=bb.im-bb2.im + bb2.re ;
         }
      }
   }

  //FFT计算结束
return;
}