home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 2 / goldfish_vol2_cd1.bin / files / misc / sci / gfft / source / cfft.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-05-25  |  12.4 KB  |  362 lines
  1. /***************************************************************************
  2.  *          Copyright (C) 1994  Charles P. Peterson                  *
  3.  *         4007 Enchanted Sun, San Antonio, Texas 78244-1254             *
  4.  *              Email: Charles_P_Peterson@fcircus.sat.tx.us                *
  5.  *                                                                         *
  6.  *          This is free software with NO WARRANTY.                  *
  7.  *          See gfft.c, or run program itself, for details.              *
  8.  *              Support is available for a fee.                      *
  9.  ***************************************************************************
  10.  *
  11.  * Program:     gfft--General FFT analysis
  12.  * File:        cfft.c
  13.  * Purpose:     fast Fourier analysis function with complex i/o array
  14.  * Author:      Charles Peterson (CPP)
  15.  * History:     4-June-1993 CPP; Created.
  16.  *              9-March-1994 CPP; added SAVE_TRIG_OPT sections.
  17.  * Comment:
  18.  *   Inspired by the algorithm expressed in four1.c, page 411-412 in
  19.  *   Numerical Recipes in C  (William H. Press, Brian P. Flannery, Saul A.
  20.  *   Teukolsky, William T. Vetterling.  Cambridge University Press, 1988.)
  21.  *   In turn, inspiration for their work came from N. Brenner of Lincoln
  22.  *   Laboratories, Danielson, Lanczos, J. W. Cooley, J. W. Tukey, and many
  23.  *   others, including Fourier, for whom the underlying analysis is named.
  24.  *   This, however, is an original implementation without any copied code.
  25.  *   I hope it is much easier to read and understand than its predecessors.
  26.  *   If nothing else, it is very verbose, and doesn't require an intuitive
  27.  *   understanding of all trigonometric identities.  Things got a little
  28.  *   worse when I added an option to save the previously used trig values,
  29.  *   but most people would consider that (in principle) worthwhile.
  30.  *
  31.  *   Single precision floating point, except for trigonometric values
  32.  *   which are computed in double precision.
  33.  */
  34.  
  35. #define SAVE_TRIG_OPT  /* Comment this out for easiest use outside GFFT */
  36. /*
  37.  * Comment the above line out for easiest use outside GFFT.
  38.  * You will lose the optional saving of trig values, but (see note below)
  39.  * that really doesn't make a big difference in performance.
  40.  *
  41.  * Otherwise, when using this code outside GFFT you will have to duplicate
  42.  * the GFFT memory allocation function, or modify this file accordingly
  43.  * (e.g. with an error return if the memory demanded cannot be allocated).
  44.  *
  45.  * Note: When SAVE_TRIG_OPT is defined, as in GFFT, there is a BOOLEAN
  46.  * global variable SaveMemory which disables the saving of trig values
  47.  * when set.
  48.  */
  49.  
  50. #include <math.h>
  51.  
  52. #ifdef _AMIGA
  53. #ifdef _M68881
  54. #include <m68881.h>
  55. #endif
  56. #endif
  57.  
  58. #include "complex.h"
  59.  
  60. #ifdef SAVE_TRIG_OPT
  61. #include "gfft.h"       /* gmalloc safe allocation (longjmps on error) */
  62. #include "settings.h"   /* SaveMemory option */
  63. #endif
  64.  
  65. void cfft (Complex_float datac[], long n, int isign)
  66. {
  67. /* ARGUMENTS:
  68.  *
  69.  * datac       input/output Complex_float data array
  70.  *             transformed according to isign (see below)
  71.  *
  72.  * n            number of complex data elements
  73.  *              *** n MUST BE A POWER OF 2 ***
  74.  *              *** THIS IS NOT CHECKED FOR ***
  75.  *
  76.  * isign    if 1, datac is replaced with its discrete Fourier transform
  77.  *              if -1, datac is replaced with n times its INVERSE
  78.  *                discrete Fourier transform
  79.  */
  80.  
  81. #ifdef SAVE_TRIG_OPT
  82.     static Complex_float *wkvector = NULL;
  83.     static last_n = 0;
  84. #endif
  85.  
  86. /* 
  87.  * This function is divided into 2 or 3 main blocks to minimize the scope
  88.  * of variables involved and hopefully to encourage better register use.
  89.  */
  90.  
  91. /* In the first block we sort the input array of complex numbers into
  92.  * bit reversal order (e.g. the value at binary index ...001 is swapped
  93.  * with the value at binary index 100...)  
  94.  *
  95.  * While i counts forwards from 1, j counts "inversely" from binary 100...
  96.  *
  97.  * Skip i==0 and i==n-1 (They never get swapped anyway)
  98.  */
  99.  
  100.     {
  101.     long const halfn = n / 2;
  102.     long const nminus1 = n - 1;
  103.     long register j = halfn;
  104.     long register i = 1;
  105.  
  106.     for (; i < nminus1; i++)
  107.     {
  108.         if (j > i)    /* If move required, and not done before */
  109.         {
  110.         CF_SWAP (datac[i], datac[j]);
  111.         }
  112.     /* 
  113.          * The following sub-block is for inverse counting.
  114.          * To count inversely, we clear highest set bits until reaching
  115.          * the first clear bit, then set it.
  116.          */
  117.         {
  118.         long register scanbit = halfn;
  119.         while (j >= scanbit && scanbit > 0)
  120.         {
  121.             j -= scanbit;
  122.             scanbit >>= 1;
  123.         }
  124.             j += scanbit;
  125.  
  126.         }   /* End inverse counting sub-block */
  127.  
  128.     }   /* End for (i = 1, 2, ..., n-1) */
  129.  
  130.     }   /* End sorting block */
  131. /*
  132.  * The following block(s) are the Danielson-Lanczos section of the routine.
  133.  * Here we combine the transforms, initially of one value each, by two
  134.  * at a time giving us transforms twice as large, until everything has
  135.  * been combined into one big transform.
  136.  *
  137.  * Thanks to clever sorting, each transform of 'even' elements is
  138.  * followed by a transform of 'odd' elements during each combination.
  139.  *
  140.  * mathematical summary (in which = denotes equality or definition):
  141.  * (A transform of length one simply copies input to output)
  142.  * F[k] = F[k][even] + (w ** k) * F[k][odd]
  143.  * (w ** k) = - (w ** (k + N/2))
  144.  * w = e ** (2 * PI * i / N)
  145.  * e ** (i * theta) = cos (theta) + i * sin (theta)
  146.  * i = sqr (-1)
  147.  * e and PI are the named transcendental numbers
  148.  * cos and sin are the named trigonometric functions
  149.  * N is the number of complex samples; k is a number 0..N-1
  150.  */
  151.  
  152. #ifdef SAVE_TRIG_OPT
  153.     if (SaveMemory)
  154. #endif
  155.  
  156.    /*
  157.     * This is the unaccelerated version
  158.     * The trig constants are computed as they are needed
  159.     * each time cfft is called.
  160.     */
  161.  
  162.     {
  163.     long osize;    /* Size of old transforms being combined */
  164.     long nsize;    /* Size of new transforms being produced */
  165.     Complex_double w;
  166.     Complex_double wtemp;
  167.     Complex_double wk;        /* (w ** k) */
  168.     Complex_float wk_times_f_k_odd;
  169.     Complex_float wkfloat;    /* Floating copy for innermost loop */
  170.  
  171.     double theta;
  172.     long ieven, iodd, k;
  173.  
  174.     for (osize = 1,nsize = 2;  osize < n;  osize = nsize,nsize *= 2)
  175.     {
  176.         theta = 2.0 * PI / (isign * nsize);  /* PI is from math.h */
  177.         w.real = cos (theta);
  178.         w.imaginary = sin (theta);
  179.         wk.real = 1.0;
  180.         wk.imaginary = 0.0;
  181.         wkfloat.real = 1.0;
  182.         wkfloat.imaginary = 0.0;
  183.  
  184.         for (k = 0; k < osize; k++)
  185.         {
  186.         for (ieven = k; ieven < n; ieven += nsize)
  187.         {
  188.             iodd = ieven + osize;
  189.         /*
  190.                  * Note that the the upper and lower parts of each new
  191.                  * transform being produced use the same components--
  192.                  * taken through a 'second cycle' since they are periodic
  193.          * in their original N (nsize/2).  Only W**k differs,
  194.          * and, fortuitously, W**k = -W**(k+nsize/2).  Also
  195.          * fortuitously, the input slots for F_k_even and F_k_odd
  196.          * are separated by nsize/2, as are the output slots for 
  197.                  * F_k and F_(k+nsize/2), so we can use and write over the 
  198.                  * same slots during each pass in the inner loop, using the
  199.                  * same value of W**k (we subtract instead of negating).
  200.          */
  201.             C_MULT (wkfloat, datac[iodd], wk_times_f_k_odd);
  202.             C_SUB (datac[ieven], wk_times_f_k_odd, datac[iodd]);
  203.             C_ADD (datac[ieven], wk_times_f_k_odd, datac[ieven]);
  204.  
  205.         } /* end for (sums using same wk factor) */
  206.  
  207.         /*  C_MULT (wk, w, wk) but must use temp or i/o coincide */
  208.  
  209.         C_MULT (wk, w, wtemp);
  210.         wkfloat.real = (float) (wk.real = wtemp.real);
  211.         wkfloat.imaginary = (float) (wk.imaginary = 
  212.                          wtemp.imaginary);
  213.  
  214.         } /* end for (k = 1, 2, ..., osize) */
  215.  
  216.     } /* end for (osize, nsize...) */
  217.  
  218.     } /* end of Danielson-Lanczos section (SaveMemory slow version) */
  219.  
  220. #ifdef SAVE_TRIG_OPT
  221.     else /* if (!SaveMemory) */
  222.     {
  223.  
  224. /* To speed things up under most cases (where the same N is used for
  225.  * more than one transform, we are going to first create the vector of
  226.  * trigonometrically derived wk constants.  If N is unchanged, and the
  227.  * the vector has already been created, we simply re-use the previous
  228.  * vector.  If it weren't for this section (which uses GFFT allocation
  229.  * routines) there would be no memory allocated in cfft, and no need for
  230.  * the gfft.h header file).  Note that this also increases the dynamic
  231.  * memory requirements of gfft significantly because the vector necessary
  232.  * must be large enough for 2*n complex floating point numbers, making it
  233.  * one of the biggest single dynamic memory uses--particularly for large N
  234.  * which some people might be interested in.  So, I have decided
  235.  * to make this feature optional, if included at all.  As currently
  236.  * implemented, it also requires the 'safe' gmalloc, which longjmps
  237.  * on error (so as not to confuse the interface here and elsewhere).
  238.  *
  239.  * 11-March-94 CPP; I'm somewhat disappointed with the 5-7% overall
  240.  * performance improvement attributable to this modification, in spite of
  241.  * the fact that lstat shows cfft as using 67% (58% accelerated) of the
  242.  * total time for a 30 sec batch spectrum analysis (run on a vanilla 68000
  243.  * with no FFP even).  But, I should have anticipated this minor change;
  244.  * the trigs become a relatively small part of the overall calculations.
  245.  * Don't let anyone tell you any different.  (Meanwhile, the dreaded
  246.  * ok_spectrum and ok_read, which look so horribly inefficient, were only
  247.  * taking around 2% of the total time each.)  The number one line is the
  248.  * complex subtraction following the main fft complex multiplication; it
  249.  * probably waits there for the complex multiplication to finish.  That,
  250.  * not surprisingly, it the chief bottleneck.  As it should be.
  251.  *
  252.  */
  253.  
  254.     if (n > 1 && (!wkvector || n != last_n))
  255.     {
  256.         long osize;    /* Size of old transforms being combined */
  257.         long nsize;    /* Size of new transforms being produced */
  258.         Complex_double w;
  259.         Complex_double wk;          /* (w ** k) */
  260.         Complex_double wtemp;
  261.  
  262.         double theta;
  263.         long k;
  264.         long vsize;  /* vector size */
  265.         long ntemp;
  266.         long wki = 0;
  267.     /*
  268.      * The vector size is N + log2(N) - 1 where log2 is 'log base 2'
  269.      * we calculate this using integer arithmetic
  270.      */
  271.         vsize = -1;
  272.         for (ntemp = n; ntemp > 1; ntemp /= 2) vsize++; /* log base 2 */
  273.         vsize += n;
  274.  
  275.         if (wkvector)
  276.         {
  277.         gfree (wkvector);
  278.         }
  279.         wkvector = gmalloc (vsize * sizeof(Complex_float), 
  280.                 NOTHING_SPECIAL);
  281.  
  282.         for (osize = 1,nsize = 2;  osize < n;  osize = nsize,nsize *= 2)
  283.         {
  284.         theta = 2.0 * PI / (isign * nsize);  /* PI is from math.h */
  285.         w.real = cos (theta);
  286.         w.imaginary = sin (theta);
  287.  
  288.         wk.real = 1.0;
  289.         wk.imaginary = 0.0;
  290.  
  291.         wkvector[wki].real = 1.0;
  292.         wkvector[wki++].imaginary = 0.0;
  293.  
  294.         for (k = 0; k < osize; k++)
  295.         {
  296.         /*  C_MULT (wk, w, wk) but must use temp or i/o coincide */
  297.  
  298.             C_MULT (wk, w, wtemp);
  299.             wkvector[wki].real = (float) (wk.real = wtemp.real);
  300.             wkvector[wki++].imaginary = (float) (wk.imaginary = 
  301.                              wtemp.imaginary);
  302.  
  303.         } /* end for (k = 1, 2, ..., osize) */
  304.  
  305.         } /* end for (osize, nsize...) */
  306.  
  307.         last_n = n;
  308.  
  309.     } /* end of trig constants computation block */
  310.  
  311.     /*
  312.      * Now, we compute this fft
  313.      */
  314.         {
  315.         long osize;    /* Size of old transforms being combined */
  316.         long nsize;    /* Size of new transforms being produced */
  317.         Complex_float wk_times_f_k_odd;
  318.  
  319.         long ieven, iodd, k;
  320.         register Complex_float *wkp = wkvector;
  321.  
  322.         for (osize = 1,nsize = 2; osize < n;  osize = nsize,nsize *= 2)
  323.         {
  324.         for (k = 0; k < osize; k++)
  325.         {
  326.             for (ieven = k; ieven < n; ieven += nsize)
  327.             {
  328.             iodd = ieven + osize;
  329.         /*
  330.          * Note that the the upper and lower parts of each new
  331.          * transform being produced use the same components--
  332.          * taken through a 'second cycle' since they are periodic
  333.          * in their original N (nsize/2).  Only W**k differs,
  334.          * and, fortuitously, W**k = -W**(k+nsize/2).  Also
  335.          * fortuitously, the input slots for F_k_even and F_k_odd
  336.          * are separated by nsize/2, as are the output slots for 
  337.          * F_k and F_(k+nsize/2), so we can use and write over the 
  338.          * same slots during each pass in the inner loop, using the
  339.          * same value of W**k (we subtract instead of negating).
  340.          */
  341.             C_MULT (*wkp, datac[iodd], wk_times_f_k_odd);
  342.             C_SUB (datac[ieven],wk_times_f_k_odd, datac[iodd]);
  343.             C_ADD (datac[ieven],wk_times_f_k_odd, datac[ieven]);
  344.  
  345.             } /* end for (sums using same wk factor) */
  346.  
  347.             wkp++;
  348.  
  349.         } /* end for (k = 1, 2, ..., osize) */
  350.  
  351.         wkp++;
  352.  
  353.         } /* end for (osize, nsize...) */
  354.  
  355.     } /* end of transform computation */
  356.  
  357.     } /* end of !SaveMemory (fast memory hog) version */
  358.  
  359. #endif /* whether SAVE_TRIG_OPT defined or not */
  360.  
  361. } /* end of cfft */
  362.