home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume1 / 8707 / 48 / fftn.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-07-13  |  3.3 KB  |  130 lines
  1. #include <math.h>
  2. #include "complex.h"
  3.  
  4. /*    Use 4-multiply complex multiply for MSC    */
  5. #ifdef DOS    /*    IBM C    */
  6. #define MULT4
  7. #endif
  8. #ifdef MSDOS    /*    Microsoft C    */
  9. #define MULT4
  10. #endif
  11.  
  12. #define PI 3.1415926535979
  13.  
  14. fftn(data, nn, ndim, isign)
  15.  
  16. #ifdef SPEED
  17. struct fcomplex data[];
  18. #else
  19. struct complex data[];
  20. #endif
  21.  
  22. unsigned nn[];
  23. int ndim, isign;
  24. {
  25.     int idim;
  26.     unsigned i1, i2rev, i3rev, ibit;
  27.     unsigned ip2, ifp1, ifp2, k2, n;
  28.     unsigned nprev = 1, nrem, ntot = 1;
  29.     register unsigned i2, i3;
  30.     double theta;
  31.     struct complex w, wp;
  32.  
  33. #ifdef SPEED
  34.     float wtemp;
  35.     struct fcomplex temp, wt;
  36. #else
  37.     double wtemp;
  38.     struct complex temp, wt;
  39. #endif
  40.  
  41. #ifndef MULT4    /*    Temporary variables needed for 3-multiply complex mult.    */
  42.     double t1, t2;
  43. #endif
  44.  
  45.     /*      Compute total number of complex values  */
  46.     for (idim = 0; idim < ndim; ++idim)
  47.         ntot *= nn[idim];
  48.  
  49.     for (idim = ndim - 1; idim >= 0; --idim) {
  50.         n = nn[idim];
  51.  
  52.         nrem = ntot / (n * nprev);
  53.         ip2 = nprev * n;        /*      Unit step for next dimension */
  54.         i2rev = 0;              /*      Bit reversed i2 */
  55.  
  56.         /*      This is the bit reversal section of the routine */
  57.         /*      Loop over current dimension     */
  58.         for (i2 = 0; i2 < ip2; i2 += nprev) {
  59.             if (i2 < i2rev)
  60.                 /*      Loop over lower dimensions      */
  61.                 for (i1 = i2; i1 < i2 + nprev; ++i1)
  62.                     /*      Loop over higher dimensions  */
  63.                     for (i3 = i1; i3 < ntot; i3 += ip2) {
  64.                         i3rev = i3 + i2rev - i2;
  65.                         temp = data[i3];
  66.                         data[i3] = data[i3rev];
  67.                         data[i3rev] = temp;
  68.                     }
  69.  
  70.             ibit = ip2;
  71.             /*      Increment from high end of i2rev to low */
  72.             do {
  73.                 ibit >>= 1;
  74.                 i2rev ^= ibit;
  75.             } while (ibit >= nprev && !(ibit & i2rev));
  76.         }
  77.  
  78.         /*      Here begins the Danielson-Lanczos section of the routine */
  79.         /*      Loop over step sizes    */
  80.         for (ifp1 = nprev; ifp1 < ip2; ifp1 <<= 1) {
  81.             ifp2 = ifp1 << 1;
  82.             /*  Initialize for the trig. recurrence */
  83.             theta = isign * 2.0 * PI / (ifp2 / nprev);
  84.             wp.x = sin(0.5 * theta);
  85.             wp.x *= -2.0 * wp.x;
  86.             wp.y = sin(theta);
  87.             w.x = 1.0;
  88.             w.y = 0.0;
  89.  
  90.             /*  Loop by unit step in current dimension  */
  91.             for (i3 = 0; i3 < ifp1; i3 += nprev) {
  92.                 /*      Loop over lower dimensions      */
  93.                 for (i1 = i3; i1 < i3 + nprev; ++i1)
  94.                     /*  Loop over higher dimensions */
  95.                     for (i2 = i1; i2 < ntot; i2 += ifp2) {
  96.                         /*      Danielson-Lanczos formula */
  97.                         k2 = i2 + ifp1;
  98.                         wt = data[k2];
  99.  
  100. #ifdef MULT4
  101. /*    Complex multiply using 4 real multiplies.  Faster in MSC    */
  102.                         data[k2].x = data[i2].x - (temp.x = w.x * wt.x - w.y * wt.y);
  103.                         data[k2].y = data[i2].y - (temp.y = w.x * wt.y + w.y * wt.x);
  104. #else
  105. /*    Complex multiply using 3 real multiplies.  Should usually be faster.    */
  106.                         data[k2].x = data[i2].x - (temp.x =
  107.                             (t1 = w.x * wt.x) - (t2 = w.y * wt.y));
  108.                         data[k2].y = data[i2].y - (temp.y =
  109.                             (w.x + w.y) * (wt.x + wt.y) - t1 - t2);
  110. #endif
  111.                         data[i2].x += temp.x;
  112.                         data[i2].y += temp.y;
  113.                     }
  114.                 /*      Trigonometric recurrence        */
  115.                 wtemp = w.x;
  116. #ifdef MULT4
  117. /*    Complex multiply using 4 real multiplies.    */
  118.                 w.x += w.x * wp.x - w.y * wp.y;
  119.                 w.y += wtemp * wp.y + w.y * wp.x;
  120. #else
  121. /*    Complex multiply using 3 real multiplies.    */
  122.                 w.x += (t1 = w.x * wp.x) - (t2 = w.y * wp.y);
  123.                 w.y += (wtemp + w.y) * (wp.x + wp.y) - t1 - t2;
  124. #endif
  125.             }
  126.         }
  127.     nprev *= n;
  128.     }
  129. }
  130.