home *** CD-ROM | disk | FTP | other *** search
/ The AGA Experience 2 / agavol2.iso / software / sound / utils / gfft-2.03 / source / gfft-2.03-source.lha / rfft.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-01-02  |  10.9 KB  |  275 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:        rfft.c
  13.  * Purpose:     fast Fourier analysis function with real input array
  14.  * Author:      Charles Peterson (CPP)
  15.  * History:     12-July-1993 CPP; Created.
  16.  *              8-Feb-95 CPP (1.31); Progress Requester
  17.  * Comment:
  18.  *   Inspired by the algorithm expressed in realft, page 417-418 in
  19.  *   Numerical Recipes in C  (William H. Press, Brian P. Flannery, Saul A.
  20.  *   Teukolsky, William T. Vetterling.  Cambridge University Press, 1988.)
  21.  *   This, however, is an original implementation without any copied code.
  22.  *
  23.  *   To fully appreciate this routine, read cfft.c, and the book, first.
  24.  *   cfft.c is where the transforming is done; here, we just use some
  25.  *   tricks to transform twice as many datapoints in (almost) the same
  26.  *   amount of time, because we are starting with real (not complex) data.
  27.  *   We start by pretending that the odd data points are the imaginary
  28.  *   values for the even data points (0, 2, ...).  Then we use a
  29.  *   symmetry modified version of the Danielson-Lanczos equation to
  30.  *   unfold the data into the positive half of the full transform.
  31.  *
  32.  *   F[k] = 1/2 * (H[k] + conj (H[N/2-k])) - 
  33.  *          i/2 * (H[k] + conj (H[N/2-k])) * e ** (2 * PI * i * k / N)
  34.  *   e ** (i * theta) = cos (theta) + i * sin (theta)
  35.  *
  36.  *   The negative half of the transform is simply related to the positive
  37.  *   half by the equation F[n] = conj (F[-n]) (because h(t) is real).
  38.  *
  39.  *   LEGEND
  40.  *   = here denotes equality
  41.  *   * denotes multiply, ** denotes power (familiar & easy to read; sorry)
  42.  *   cconj is the complex conjugate: conj( {x, y} ) = {x, -y}
  43.  *   i = sqr (-1)
  44.  *   e and PI are the named transcendental numbers
  45.  *   cos and sin are the named trigonometric functions
  46.  *   N is the total number of samples input; k is a number 0..N-1
  47.  *   F is the output complex transform; H is the pseudo-transform returned
  48.  *       by sending cfft the packed samples
  49.  *
  50.  *   Single precision floating point, except for trigonometric values
  51.  *   which are computed in double precision.
  52.  *
  53.  *   The values F[0] and F[N/2] are real and independent.  In order to
  54.  *   pack all data into the original space, Press, et. al., chose to
  55.  *   put the value for F[N/2] into the space for the imaginary part of
  56.  *   F[0].  I didn't like this idea at first, but as the higher level
  57.  *   procedures are written to be general, creating and maintaining an
  58.  *   array 1 unit larger for the 'real' case gets to be quite complicated.
  59.  *   And, it's so nice when input and output arrays are exactly the same
  60.  *   size always.  (I'm breathing easier already.)
  61.  *
  62.  *   (-: CRACKPOT OPINION AND RATIONALIZATION MODE ON
  63.  *
  64.  *   Heck, I shouldn't be programming stuff like this in the psuedo-
  65.  *   assembler known as C.  I want a nice high level OO language with
  66.  *   dynamic arrays and such.  On the other hand, I hear this little voice
  67.  *   saying I'm a fool not to use an assembler--particularly one tailored
  68.  *   to a particular DSP.  I argue critical parts can always be recoded,
  69.  *   and right now I want to get the algorithms and structure correct, and
  70.  *   make this a shining example of readable, understandable, and
  71.  *   learnable--yes, literate--code so that even if people don't use it,
  72.  *   some may learn from it.  Unfortunately, we've not yet been delivered
  73.  *   to the Promised Land of perfect languages, with high and low levels
  74.  *   seemlessly integrated and unified by intuitive concepts and syntax
  75.  *   expanding to infinity.  Hmmn, maybe God didn't want that to happen
  76.  *   yet.  Anyway, right now C is the closest thing to UNIversal, and a
  77.  *   fair compromise language too.  And compromises one has to keep on
  78.  *   making in C as well--they go with the turf.  Dammit.
  79.  *
  80.  *   CRACKPOT OPINION AND RATIONALIZATION MODE OFF :-)
  81.  *
  82.  *   Another catch...the input data is real, and the output is (mostly)
  83.  *   complex.  And, since we don't want to allocate another array and copy
  84.  *   data around, what do we do?  Well, Press, et. al., simply always use
  85.  *   float arrays and use a 'convention' as to how the complex data is
  86.  *   stored therein.  (They did that for four1--the equivalent of my
  87.  *   cfft--as well, which I do consider Blasphemy, since all data there is
  88.  *   complex.)  I use Complex_float arrays there and here as well, so at
  89.  *   least I am correct 3 out of 4 times.  Not only that, but it makes the
  90.  *   code much more traceable to the algorithms, I believe.  (Check it out
  91.  *   for yourself.  I went to a lot of trouble writing cfft.c and rfft.c,
  92.  *   so I hope somebody reads them, even if they don't bother to read the
  93.  *   comments.)  But how do you call this function, then?  Well, you could
  94.  *   allocate a Complex_float data array to begin with, and pack the real
  95.  *   input data into it initially.  Or, you could try declaring rfft()
  96.  *   differently in the caller--using a float data array for i/o.  (This
  97.  *   passed my current SAS 6.0 'ANSI' compiler and linker--without even a
  98.  *   warning--until I put a prototype for realft in complex.h.)  Or, you
  99.  *   could cast the data array/pointer to Complex_float at the point of
  100.  *   call in the caller.  That worked for SAS 6.0 as well, with an
  101.  *   appropriate warning.  I like that approach; I hope they don't change
  102.  *   the compiler to make it not work, and I hope GCC and other compilers
  103.  *   are equally forgiving.  Of course, I take no responsibility for any of
  104.  *   this, or anything.  I know nothing, nuuthing.
  105.  *
  106.  */
  107.  
  108. #include <math.h>
  109.  
  110. #ifdef _AMIGA
  111. #ifdef _M68881
  112. #include <m68881.h>
  113. #endif
  114. #endif
  115.  
  116. #include "complex.h"
  117.  
  118. void rfft (Complex_float datac[], long n, int isign)
  119. {
  120. /* ARGUMENTS:
  121.  *
  122.  * datac       input/output data array
  123.  * -----       transformed according to isign (see below)
  124.  *
  125.  *             * (For isign = 1) Input array is a real array packed (?) into
  126.  *             * the complex array.  For i = 0  to  n, a hypothetical input
  127.  *             * array dataf[0..n] is packed as follows:
  128.  *             * datac[i/2].real = dataf[i]
  129.  *             * datac[i/2].imaginary = dataf[i+1]
  130.  *             * Or, if you have a reasonable compiler, you might simply
  131.  *             * be able to cast the real (float) array as a the required
  132.  *             * Complex_float array and be done with it.  (Please check
  133.  *             * that it is working properly if you do that.)
  134.  *
  135.  *             * Output array is the positive frequency half of the
  136.  *             * full transform, i.e. F[0]..F[N/2].  It is complex, so
  137.  *             * the declared type Complex_float applies.
  138.  *
  139.  *             *** NOTE: The value for F[0] is real (not complex) and is
  140.  *             *** returned in datac[0].real.  The value for F[N/2] is
  141.  *             *** real (not complex) and is returned in datac[0].imaginary.
  142.  *
  143.  *             * (For isign = -1, the roles of input and output above are
  144.  *             * reversed.)
  145.  *
  146.  *  n          number of real data elements packed in input
  147.  * ---         (2 less than half the number of complex numbers to output)
  148.  *             *** n MUST BE A POWER OF 2  ***
  149.  *             *** THIS IS NOT CHECKED FOR ***
  150.  *
  151.  * isign    if 1, datac is replaced with its discrete Fourier transform
  152.  * -----        if -1, datac is replaced with n times its INVERSE
  153.  *              discrete Fourier transform
  154.  *
  155.  */
  156.  
  157.     Complex_double w;
  158.     Complex_double wk;
  159.     Complex_double wktemp;
  160.     Complex_float wkfloat;
  161.     double theta;
  162.     long k;
  163.  
  164.     long halfn = n / 2;
  165.     long halfhalfn = halfn / 2;
  166.     float shalf = 0.5;
  167.  
  168.     if (n < 2)
  169.     {
  170.     return;  /* Nothing like satisfying the boundary condition.    */
  171.              /* This is the correct result too (as long as you     */
  172.              /* i/o only one value: datac[0].real), and avoids a   */
  173.              /* nasty divide-by-zero and use of uninitialized data */
  174.     }
  175.  
  176.     theta = PI / (double) (isign * halfn);
  177.  
  178.     if (isign == 1)
  179.     {
  180.     cfft (datac, halfn, isign);
  181.     }
  182.     else
  183.     {
  184.     /*
  185.      * Inverse transforming is done near end of routine.
  186.      * For now, just reset this one factor which differs between
  187.      * normal and inverse cases
  188.      * (set explicitly to encourage multiply-by-half optimizations).
  189.      */
  190.     shalf = (-0.5);
  191.     }
  192. /*
  193.  * We're skipping k=0 and so we get right down to trigonometric values
  194.  */
  195.     wkfloat.real = (float) (wk.real = w.real = cos (theta));
  196.  
  197.     wkfloat.imaginary = (float) (wk.imaginary = w.imaginary = sin (theta));
  198.  
  199.     for (k = 1; k < halfhalfn; k++)
  200.     {
  201.     int n2k = halfn - k;
  202.     Complex_float hk;
  203.     Complex_float hn2k;
  204.     /*
  205.      * First, we peel the two transforms (for k and N/2-k) apart
  206.      */
  207.         hk.real = 0.5 * (datac[k].real + datac[n2k].real);
  208.     hk.imaginary = 0.5 * (datac[k].imaginary - datac[n2k].imaginary);
  209.     hn2k.real = shalf * (datac[k].imaginary + datac[n2k].imaginary);
  210.     hn2k.imaginary = (-shalf) * (datac[k].real - datac[n2k].real);
  211.     /*
  212.      * Then, we recombine them using the modified D/L equation
  213.      */
  214.         datac[k].real = hk.real + (hn2k.real * wkfloat.real) -
  215.       (hn2k.imaginary * wkfloat.imaginary);
  216.     datac[k].imaginary = hk.imaginary + (hn2k.imaginary * 
  217.                          wkfloat.real) +
  218.       (hn2k.real * wkfloat.imaginary);
  219.  
  220.     datac[n2k].real = hk.real - (hn2k.real * wkfloat.real) +
  221.       (hn2k.imaginary * wkfloat.imaginary);
  222.     datac[n2k].imaginary = (-hk.imaginary) +
  223.       (hn2k.imaginary * wkfloat.real) + (hn2k.real * wkfloat.imaginary);
  224.     /*
  225.      * Now, compute the next power of e
  226.      */
  227.         C_MULT (wk, w, wktemp);  /* i/o must not overlap */
  228.     wkfloat.real = (float) (wk.real = wktemp.real);
  229.     wkfloat.imaginary = (float) (wk.imaginary = wktemp.imaginary);
  230.     }
  231. /*
  232.  * Now we compute and store the values for 0 and N/2
  233.  */
  234.     if (isign == 1)
  235.     {
  236.     float tempr = datac[0].real;
  237.     datac[0].real += datac[0].imaginary;
  238.     datac[0].imaginary = tempr - datac[0].imaginary;
  239.     }
  240.     else
  241.     {
  242.     /*
  243.      * Now, here's the completion of the inverse part
  244.      * You'll have to take this part on faith (I am)
  245.      */
  246.     float tempr = datac[0].real;
  247.     datac[0].real = 0.5 * (tempr + datac[0].imaginary);
  248.     datac[0].imaginary = 0.5 * (tempr - datac[0].imaginary);
  249.  
  250.     cfft (datac, halfn, isign);
  251.     }
  252. }
  253.  
  254. /*
  255.  * The number of inner loops would be N log2 N where N is # of bins
  256.  * However, thanks to the above function, the number is (N/2) log N
  257.  * That's why fft_inner_loops is defined here
  258.  * 
  259.  */
  260.  
  261. long fft_inner_loops (long number_bins)
  262. {
  263.     long log2;
  264.     long alog;
  265.     long loopcount;
  266.  
  267.     for (alog = 1, log2 = 0; 
  268.      alog < number_bins; 
  269.      log2++, alog *= 2);
  270.  
  271.     loopcount = (number_bins / 2) * log2;
  272.  
  273.     return loopcount;
  274. }      
  275.