home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume1 / 8707 / 75 / fft.c < prev   
Encoding:
C/C++ Source or Header  |  1990-07-13  |  3.2 KB  |  150 lines
  1. /*
  2.  *    fft.c
  3.  *
  4.  *    Performs a Cooley-Tukey Fast Fourier Transform.
  5.  *
  6.  *    Original version Pascal written in "Simple Calculations with 
  7.  *    Complex Numbers by David Clark Doctor Dobbs Journal 10/84
  8.  *
  9.  *    Translated to C by R. Hellman 02/21/86
  10.  *    Released to Public Domain
  11.  *
  12.  *    Difference from original version:
  13.  *
  14.  *    As in fortran versions that I have worked with, instead of
  15.  *    specifying total array size in decimal - instead specify number
  16.  *    of array members as power of two.  Forward or reverse fft is
  17.  *    specified by the sign of this argument.
  18.  *
  19.  *    ****IMPORTANT****
  20.  *    For convenience sake members of an array are taken to begin at [1]
  21.  *    rather than [0] thus for an fft on a 64 member array of pointers
  22.  *    to complex numbers, there must be  65 members.
  23.  *
  24.  *    f is a one dimensional array of pointers to complex numbers with
  25.  *    length (2**abs(ln)).  The sign of ln determines in which direction
  26.  *    the transform will be performed.  If ln is positive, then isinverse
  27.  *    is set to -1 and a forward transform is carried out.  If ln is
  28.  *    negative then isinverse is set to 1 and a reverse transform is
  29.  *    calculated.
  30.  *
  31.  *    For numpoints = (2**abs(ln)),
  32.  *        transform[j] = sum(f[i] * w ^ ((i-1) * (j-1))),
  33.  *    where i, j are 1 to numpoints and
  34.  *        w = exp(isinverse * 2 * pi * sqrt(-1)/numpoints).
  35.  *
  36.  *    The program also computes the inverse transform, for which the
  37.  *    defining expression is:
  38.  *        inverse DFT = (1/numpoints) * sum(f[i] * w ^ ((i-1) * (j-1))).
  39.  *
  40.  *    Run time is proportional to numpoints * log2(numpoints) rather than
  41.  *    to numpoints ** 2 for the classical discrete Fourier transform.
  42.  */
  43.  
  44. #include <stdio.h>
  45. #include "complex.c"
  46.  
  47. fft(f, ln)
  48. complex    *f[];
  49. int    ln;
  50. {
  51.     int    i, isinverse, numpoints;
  52.  
  53.     if (ln > 0)  {
  54.             puts("Calculating forward Fourier transform\n");
  55.             isinverse = -1;
  56.             numpoints = pow(2.0, ln);
  57.     }
  58.         else {
  59.             puts("Calculating reverse Fourier transform\n");
  60.             isinverse = 1;
  61.             numpoints = (int) pow(2.0, -ln);
  62.         }
  63.  
  64.     bitreverse(numpoints, f);
  65.     butterfly(numpoints, isinverse, f);
  66. }
  67.  
  68. /*
  69.  *    Reverse the bits
  70.  *
  71.  */
  72.  
  73. static bitreverse(numpoints, f)
  74. int    numpoints;
  75. complex    *f[];
  76. {
  77.     int    i, j, m;
  78.     complex    *temp;
  79.  
  80.     j = 1;
  81.     for (i = 1; i <= numpoints; i++)  {
  82.         if (i < j)  {
  83.             temp = f[j];
  84.             f[j] = f[i];
  85.             f[i] = temp;
  86.         }
  87.  
  88.         m = numpoints >> 1;
  89.  
  90.         if (m < j)
  91.             while (m < j)  {
  92.                 j -= m;
  93.                 m = (m + 1) >> 1;
  94.             }
  95.  
  96.         j += m;
  97.     }
  98. }
  99.  
  100. /*
  101.  *    Calculate the butterflies for the bit reversed data
  102.  *    Normalize if a reverse fft is performed
  103.  *
  104.  */
  105.  
  106. static butterfly(numpoints, isinverse, f)
  107. int    numpoints, isinverse;
  108. complex    *f[];
  109. {
  110.     int    mmax, step, index, m, i, j;
  111.     double    theta;
  112.     complex    w, cn, temp;
  113.  
  114.     mmax = 1;
  115.  
  116.     while (numpoints > mmax)  {
  117.  
  118.         step = mmax << 1;
  119.  
  120.         for (m = 1; m <= mmax; m++)  {
  121.  
  122.             theta = M_PI * ((double)(isinverse) *
  123.                        (double)(m - 1)) / (double)(mmax);
  124.  
  125.             w.re = cos(theta);
  126.             w.im = sin(theta);
  127.  
  128.             for (i = 1; i <= ((numpoints-m) / step) + 1; i++)  {
  129.  
  130.                 index = m + ((i-1) * step);
  131.                 j = index + mmax;
  132.  
  133.                 cmult(&temp, &w, f[j]);
  134.                 csub(f[j], f[index], &temp);
  135.                 cadd(f[index], f[index], &temp);
  136.             }
  137.         }
  138.  
  139.         mmax = step;
  140.     }
  141.  
  142.     if (isinverse == 1)  {
  143.         cn.re = numpoints;
  144.         cn.im = 0.0;
  145.  
  146.         for (i = 1; i <= numpoints; i++)
  147.             cdiv(f[i], f[i], &cn);
  148.     }
  149. }
  150.