home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource3 / 155_01 / fft.c < prev    next >
Text File  |  1980-01-01  |  4KB  |  151 lines
  1. /*    
  2. HEADER:
  3. TITLE:        Fast Fourier Transform;
  4. DATE:        05/18/1985;
  5. DESCRIPTION:    "Performs fast fourier transform using method described
  6.         by E. O. Brigham.  For details of the method, refer
  7.         to Brigham's book. THE FAST FOURIER TRANSFORM";    
  8. KEYWORDS:     Fourier, transform;
  9. FILENAME:    FFT.C;
  10. WARNINGS:    
  11.   "This program is self-contained, all that is needed is a manner of getting
  12.   the data into the array real_data (& imag_data, if applicable).  The 
  13.   transformed data will reside in these two arrays upon return with the 
  14.   original data being destroyed."
  15. AUTHORS:    Jim Pisano;
  16. COMPILERS:    DeSmet;
  17. REFERENCES:    AUTHORS:    E. O. Brigham;
  18.         TITLE:        "THE FAST FOURIER TRANSFORM";
  19.         CITATION:    "";
  20.     ENDREF
  21. */
  22.  
  23. /*    file name fft.c
  24. *    program name fft() ... Fast Fourier Transform
  25. *
  26. *    Perform fast fourier transform using method described by E. O. Brigham.
  27. *  For details of the method, refer to Brigham's book 
  28. *
  29. *    Translated to C from FORTRAN by
  30. *
  31. *        Jim Pisano
  32. *        P.O. Box 3134
  33. *        University Station
  34. *        Charlottesville, VA 22903
  35. *
  36. *  This program is in the public domain & may be used by anyone for commercial
  37. *  or non-commercial purposes.
  38. *
  39. *  real_data ... ptr. to real part of data to be transformed
  40. *  imag_data ... ptr. to imag  "   "   "   "  "      "
  41. *  inv ..... Switch to flag normal or inverse transform
  42. *  n_pts ... Number of real data points
  43. *  nu ...... logarithm in base 2 of n_pts e.g. nu = 5 if n_pts = 32. 
  44. *
  45. *  This program is self-contained, all that is needed is a manner of getting
  46. *  the data into the array real_data (& imag_data, if applicable).  The 
  47. *  transformed data will reside in these two arrays upon return with the 
  48. *  original data being destroyed.
  49. */
  50.  
  51. #include    <stdio.h>
  52. #include    <math.h>            /* declare math functions to use */
  53. #define        PI        3.1419527
  54.  
  55. fft( real_data, imag_data, n_pts, nu, inv )
  56. double *real_data, *imag_data;
  57. int n_pts, nu, inv;
  58. {
  59.     int n2, j, j1, l, i, ib, k, k1, k2;
  60.     int sgn;
  61.     double tr, ti, arg, nu1;    /* intermediate values in calcs. */
  62.     double c, s;           /* cosine & sine components of Fourier trans. */
  63.  
  64.     n2 = n_pts / 2;
  65.     nu1 = nu - 1.0;
  66.     k = 0;
  67. /* 
  68. * sign change for inverse transform 
  69. */
  70.     sgn = inv ? -1 : 1;        
  71. /*
  72. * Calculate the componets of the Fourier series of the function
  73. */
  74.     for( l = 0; l != nu; l++ )
  75.     {
  76.         do
  77.         {
  78.             for( i = 0; i != n2; i++ )
  79.             {    
  80.                 j = k / ( pow( 2.0, nu1 ) );
  81.                 ib = bit_swap( j, nu );
  82.                 arg = 2.0 * PI * ib / n_pts;
  83.                 c = cos( arg );
  84.                 s = sgn * sin( arg );
  85.                 k1 = k;
  86.                 k2 = k1 + n2;
  87.                 tr = *(real_data+k2) * c + *(imag_data+k2) * s;
  88.                 ti = *(imag_data+k2) * c - *(real_data+k2) * s;
  89.                 *(real_data+k2) = *(real_data+k1) - tr;
  90.                 *(imag_data+k2) = *(imag_data+k1) - ti;
  91.                 *(real_data+k1) = *(real_data+k1) + tr;
  92.                 *(imag_data+k1) = *(imag_data+k1) + ti;
  93.                 k++;
  94.             }
  95.             k +=  n2;
  96.         } while( k < n_pts - 1);
  97.         k = 0;
  98.         nu1 -= 1.0;
  99.         n2 /= 2;
  100.     }    
  101.     for( k = 0; k != n_pts; k++ )
  102.     {
  103.         ib = bit_swap( k, nu );
  104.         if( ib > k)
  105.         {
  106.             swap( (real_data+k), (real_data+ib) );
  107.             swap( (imag_data+k), (imag_data+ib) );
  108.         }
  109.     }
  110. /* 
  111. * If calculating the inverse transform, must divide the data by the number of
  112. * data points.
  113. */
  114.     if( inv )
  115.         for( k = 0; k != n_pts; k++)
  116.         {
  117.             *(real_data+k) /= n_pts;
  118.             *(imag_data+k) /= n_pts;
  119.         } 
  120. }
  121. /* 
  122. * Bit swaping routine in which the bit pattern of the integer i is reordered.
  123. * See Brigham's book for details
  124. */
  125. bit_swap( i, nu )
  126. int i, nu;
  127. {
  128.     int ib, i1, i2;
  129.  
  130.     ib = 0;
  131.  
  132.     for( i1 = 0; i1 != nu; i1++ )
  133.     {
  134.         i2  = i / 2;
  135.         ib = ib * 2 + i - 2 * i2;
  136.         i   = i2;
  137.     }
  138.     return( ib );
  139. }
  140. /*
  141. * Simple exchange routine where *x1 & *x2 are swapped
  142. */
  143. swap( x1, x2 )    
  144. double *x1, *x2;
  145. {
  146.     double *temp_x;
  147.  
  148.     *temp_x = *x1;
  149.     *x1 = *x2;
  150.     *x2 = *temp_x;
  151. }