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