/*---------------------------------------------------------------------------- 
 * 
 * NAME 
 *      fft1  - A one dimensional FFT subroutine. 
 * 
 * USAGE 
 *      fft1 (npoints, real, imag, dir) 
 *          int npoints, dir; 
 * 
 * EXPLICIT INPUTS 
 *      npoints:        the number of points in the FFT. Must be a power of 2. 
 *      real, imag:     pointers to arrays of doubles for input and output. 
 *      dir:            1 for forward transform 
 *                      -1 for inverse transform 
 * 
 * IMPLICIT INPUTS 
 *      None. 
 * 
 * EXPLICIT OUTPUTS 
 *      None. 
 * 
 * IMPLICIT OUTPUTS 
 *         The contents of the arrays are changed from the input data to the 
 *      FFT coefficients. 
 * 
 * PROCEDURE 
 *         Uses time decomposition with input bit reversal. The Cooley/Tookey 
 *      Fortran scheme for doing recursive odd/even decimation without really 
 *      doing bit reversal is used. The computation is done in place, so the 
 *      output replaces the input. 
 * 
 * HISTORY 
 * 16-Nov-82  John Schlag (jfs) at Carnegie-Mellon University 
 *      Incorporated O'Gorman's code to generalize the algorithm for 
 *      both forward and inverse use. 
 * 
 * 29-Jun-82  John Schlag (jfs) at Carnegie-Mellon University 
 *      Created back in December, 1980 for an EE631 project at the University 
 *      of Delaware. Translated from Fortran into C. Entered to the CMU system 
 * 
 *---------------------------------------------------------------------------- 
 */ 
 
#include <math.h> 
#include "image.h" 
#include "bblib.h" 
#ifdef TABLE 
#include "sincos.h" 
#endif 
 
#ifdef TABLE 
static double tsin(int angle) 
{ 
    if (angle < 0) 
    { 
  if (angle < -32) angle = - 64 - angle; 
 return -sinus[-angle]; 
    } 
    else 
    { 
 if (angle > 32) angle = 64 - angle; 
 return sinus[angle]; 
    } 
} 
 
static double tcos(int angle) 
{ 
    if (angle > 0) angle = -angle; 
 angle += 32; 
    if (angle < 0) 
 return -sinus[-angle]; 
    else 
 return sinus[angle]; 
} 
#endif 
 
int 
fft1 (int npoints, double real[], double imag[], int dir) 
{ 
    register int  i, j, k; 
    int   index, swapindex; 
    register double  tr, ti, wr, wi; 
#ifndef TABLE 
    register double  angle; 
#else 
    register int  angle; 
#endif 
 
    /* SWAP THE INPUT ELEMENTS FOR THE DECIMATION IN TIME ALGORITHM. */ 
 
    for (index = 1, swapindex = 0; index < npoints; index++) 
    { 
        k = npoints; 
        do 
        { 
     k /= 2; 
 } while ( (swapindex + k) >= npoints ); 
        swapindex = (swapindex % k) + k; 
        if (swapindex <= index)  continue; 
            tr = real[index]; 
        real[index] = real[swapindex]; 
        real[swapindex] = tr; 
        ti = imag[index]; 
        imag[index] = imag[swapindex]; 
        imag[swapindex] = ti; 
    } 
 
    /* 
     * DO THE BUTTERFLY COMPUTATIONS. 
     *      stage index:    k = 1, 2, 4, 8 . . . 
     *                      For npoints = 8, for example, there will 
     *                      be three stages of butterfly computations. 
     *      b'fly indices:  i and j will be separated by a distance 
     *                      dependent on the current stage. 
     *                      k is used as the separation constant. 
     */ 
 
    for (k = 1; k < npoints; k *= 2) 
    { 
        for (index = 0; index < k; index++) 
        { 
 
#ifndef TABLE 
            angle = M_PI * ((double) (index * -dir)) / ((double) k); 
            wi = sin(angle);  
            wr = cos(angle); 
#else 
     angle = index * -dir * 64 / k; 
     wi = tsin(angle); 
     wr = tcos(angle); 
#endif 
            for (i = index; i < npoints; i += 2*k) 
            { 
                j = i + k; 
                tr = (wr * real[j]) - (wi * imag[j]); 
                ti = (wr * imag[j]) + (wi * real[j]); 
                real[j] = real[i] - tr; 
                imag[j] = imag[i] - ti; 
                real[i] += tr; 
                imag[i] += ti; 
            } 
        } 
    } 
 
    /* For the inverse transform, scale the output by 1/N. */ 
 
    if (dir == -1) 
    { 
       for (i = 0; i < npoints; i++) 
       { 
            real[i] /= npoints; 
            imag[i] /= npoints; 
       } 
    } 
    return NO_ERROR; 
} 
 
#if 0  
 
/*  
 * This is used to compare the accuracy of a table lookup method 
 * versus the actual sinus/cosinus computation. (rdv) 
 */ 
 
main(void)
{ 
double angle, wi1, wi2, wr1, wr2; 
int iangle, i; 
 
 for (i=-63; i<64; i++) 
 { 
                angle = M_PI * (double) i / (double) 64; 
                wi1 = sin(angle); 
                wr1 = cos(angle); 
  iangle = i; 
  wi2 = tsin(iangle); 
  wr2 = tcos(iangle); 
  printf("sin: %g %g\n", wi1, wi2); 
  printf("cos: %g %g\n", wr1, wr2); 
 }  
} 
#endif 
