fft_code.txt

对fft变化的简单应用

/*  
    * fft.c  
    *  
    * Unix   Version   2.4   by   Steve   Sampson,   Public   Domain,   September   1988  
    *  
    * This   program   produces   a   Frequency   Domain   display   from   the   Time   Domain  
    * data   input;   using   the   Fast   Fourier   Transform.  
    *  
    * The   Real   data   is   generated   by   the   in-phase   (I)   channel,   and   the  
    * Imaginary   data   is   produced   by   the   quadrature-phase   (Q)   channel   of  
    * a   Doppler   Radar   receiver.     The   middle   filter   is   zero   Hz.     Closing  
    * targets   are   displayed   to   the   right,   and   Opening   targets   to   the   left.  
    *  
    * Note:   With   Imaginary   data   set   to   zero   the   output   is   a   mirror   image.  
    *  
    * Usage: fft   samples   input   output  
    * 1.   samples   is   a   variable   power   of   two.  
    * 2.   Input   is   (samples   *   sizeof(double))   characters.  
    * 3.   Output   is   (samples   *   sizeof(double))   characters.  
    * 4.   Standard   error   is   help   or   debugging   messages.  
    *  
    * See   also:   readme.doc,   pulse.c,   and   sine.c.  
    */  
   
  /*   Includes   */  
   
  #include   <stdio.h>  
  #ifdef   AZTEC_C  
  #include   <strings.h>  
  #else  
  #include   <malloc.h>  
  #endif  
  #include   <math.h>  
  #include   <stdlib.h>  
   
  /*   Defines   */  
   
  #define TWO_PI (2.0   *   3.14159265358979324) /*   alias   360   deg   */  
  #define Chunk (Samples   *   sizeof(double))  
   
  #define sine(x) Sine[(x)]  
  #define cosine(x) Sine[((x)   +   (Samples   >>   2))   %   Samples]  
   
  /*   Globals,   Forward   declarations   */  
   
  static   int Samples,   Power;  
  static   double *Real,   *Imag,   Maxn,   magnitude();  
  static   void scale(),   fft(),   max_amp(),   display(),   err_report();  
   
  static   int permute();  
  static   double *Sine;  
  static   void build_trig();  
   
  static   FILE *Fpi,   *Fpo;  
   
   
  /*   The   program   */  
   
  main(argc,   argv)  
  int argc;  
  char **argv;  
  {  
  if   (argc   !=   4)  
  err_report(0);  
   
  Samples   =   abs(atoi(*++argv));  
  Power   =   log10((double)Samples)   /   log10(2.0);  
   
  /*   Allocate   memory   for   the   variable   arrays   */  
   
  if   ((Real   =   (double   *)malloc(Chunk))   ==   NULL)  
  err_report(1);  
   
  if   ((Imag   =   (double   *)malloc(Chunk))   ==   NULL)  
  err_report(2);  
   
  if   ((Sine   =   (double   *)malloc(Chunk))   ==   NULL)  
  err_report(3);  
   
  /*   open   the   data   files   */  
   
  if   ((Fpi   =   fopen(*++argv,   "r"))   ==   NULL)  
  err_report(4);  
   
  if   ((Fpo   =   fopen(*++argv,   "w"))   ==   NULL)  
  err_report(5);  
   
  /*   read   in   the   data   from   the   input   file   */  
   
  fread(Real,   sizeof(double),   Samples,   Fpi);  
  fread(Imag,   sizeof(double),   Samples,   Fpi);  
  fclose(Fpi);  
   
  build_trig();  
  scale();  
  fft();  
  display();  
   
  /*   write   the   frequency   domain   data   to   the   standard   output   */  
   
  fwrite(Real,   sizeof(double),   Samples,   Fpo);  
  fwrite(Imag,   sizeof(double),   Samples,   Fpo);  
  fclose(Fpo);  
   
  /*   free   up   memory   and   let's   get   back   to   our   favorite   shell   */  
   
  free((char   *)Real);  
  free((char   *)Imag);  
  free((char   *)Sine);  
   
  exit(0);  
  }  
   
   
  static   void   err_report(n)  
  int n;  
  {  
  switch   (n)     {  
  case   0:  
  fprintf(stderr,   "Usage:   fft   samples   in_file   out_file\n");  
  fprintf(stderr,   "Where   samples   is   a   power   of   two\n");  
  break;  
  case   1:  
  fprintf(stderr,   "fft:   Out   of   memory   getting   real   space\n");  
  break;  
  case   2:  
  fprintf(stderr,   "fft:   Out   of   memory   getting   imag   space\n");  
  free((char   *)Real);  
  break;  
  case   3:  
  fprintf(stderr,   "fft:   Out   of   memory   getting   sine   space\n");  
  free((char   *)Real);  
  free((char   *)Imag);  
  break;  
  case   4:  
  fprintf(stderr,"fft:   Unable   to   open   data   input   file\n");  
  free((char   *)Real);  
  free((char   *)Imag);  
  free((char   *)Sine);  
  break;  
  case   5:  
  fprintf(stderr,"fft:   Unable   to   open   data   output   file\n");  
  fclose(Fpi);  
  free((char   *)Real);  
  free((char   *)Imag);  
  free((char   *)Sine);  
  break;  
  }  
   
  exit(1);  
  }  
   
   
  static   void   scale()  
  {  
  register   int loop;  
   
  for   (loop   =   0;   loop   <   Samples;   loop++)     {  
  Real[loop]   /=   Samples;  
  Imag[loop]   /=   Samples;  
  }  
  }  
   
   
  static   void   fft()  
  {  
  register   int loop,   loop1,   loop2;  
  unsigned i1; /*   going   to   right   shift   this   */  
  int i2,   i3,   i4,   y;  
  double a1,   a2,   b1,   b2,   z1,   z2;  
   
  i1   =   Samples   >>   1;  
  i2   =   1;  
   
  /*   perform   the   butterfly's   */  
   
  for   (loop   =   0;   loop   <   Power;   loop++)     {  
  i3   =   0;  
  i4   =   i1;  
   
  for   (loop1   =   0;   loop1   <   i2;   loop1++)     {  
  y   =   permute(i3   /   (int)i1);  
  z1   =     cosine(y);  
  z2   =   -sine(y);  
   
  for   (loop2   =   i3;   loop2   <   i4;   loop2++)     {  
   
  a1   =   Real[loop2];  
  a2   =   Imag[loop2];  
   
  b1     =   z1*Real[loop2+i1]   -   z2*Imag[loop2+i1];  
  b2     =   z2*Real[loop2+i1]   +   z1*Imag[loop2+i1];  
   
  Real[loop2]   =   a1   +   b1;  
  Imag[loop2]   =   a2   +   b2;  
   
  Real[loop2+i1]   =   a1   -   b1;  
  Imag[loop2+i1]   =   a2   -   b2;  
  }  
   
  i3   +=   (i1   <<   1);  
  i4   +=   (i1   <<   1);  
  }  
   
  i1   >>=   1;  
  i2   <<=   1;  
  }  
  }  
   
  /*  
    * Display   the   frequency   domain.  
    *  
    * The   filters   are   aranged   so   that   DC   is   in   the   middle   filter.  
    * Thus   -Doppler   is   on   the   left,   +Doppler   on   the   right.  
    */  
   
  static   void   display()  
  {  
  register   int loop,   offset;  
  int n,   x;  
   
  max_amp();  
  n   =   Samples   >>   1;  
   
  for   (loop   =   n;   loop   <   Samples;   loop++)     {  
  x   =   (int)(magnitude(loop)   *   59.0   /   Maxn);  
  printf("%d\t|",   loop   -   n);  
  offset   =   0;  
  while   (++offset   <=   x)  
  putchar('=');  
   
  putchar('\n');  
  }  
   
  for   (loop   =   0;   loop   <   n;   loop++)     {  
  x   =   (int)(magnitude(loop)   *   59.0   /   Maxn);  
  printf("%d\t|",   loop   +   n);  
  offset   =   0;  
  while   (++offset   <=   x)  
  putchar('=');  
   
  putchar('\n');  
  }  
  }  
   
  /*  
    * Find   maximum   amplitude  
    */  
   
  static   void   max_amp()  
  {  
  register   int loop;  
  double mag;  
   
  Maxn   =   0.0;  
  for   (loop   =   0;   loop   <   Samples;   loop++)     {  
  if   ((mag   =   magnitude(loop))   >   Maxn)  
  Maxn   =   mag;  
  }  
  }  
   
  /*  
    * Calculate   Power   Magnitude  
    */  
   
  static   double   magnitude(n)  
  int n;  
  {  
  n   =   permute(n);  
  return   hypot(Real[n],   Imag[n]);  
  }  
   
  /*  
    * Bit   reverse   the   number  
    *  
    * Change   11100000b   to   00000111b   or   vice-versa  
    */  
   
  static   int   permute(index)  
  int index;  
  {  
  register   int loop;  
  unsigned n1;  
  int result;  
   
  n1   =   Samples;  
  result   =   0;  
   
  for   (loop   =   0;   loop   <   Power;   loop++)     {  
  n1   >>=   1;  
  if   (index   <   n1)  
  continue;  
   
  /*   Unix   has   a   truncation   problem   with   pow()   */  
   
  result   +=   (int)(pow(2.0,   (double)loop)   +   .05);  
  index   -=   n1;  
  }  
   
  return   result;  
  }  
   
  /*  
    * Pre-compute   the   sine   and   cosine   tables  
    */  
   
  static   void   build_trig()  
  {  
  register   int loop;  
  double angle,   increment;  
   
  angle   =   0.0;  
  increment   =   TWO_PI   /   (double)Samples;  
   
  for   (loop   =   0;   loop   <   Samples;   loop++)     {  
  Sine[loop]   =   sin(angle);  
  angle   +=   increment;  
  }  
  }