/******************************************************************************

   fourier.cp    Fast Fourier Transform for n = 2 raised to the power n2pow
   ----------    n2pow is between 1 and sizeof(long) * 8 - 1 (usually 31) 


   How does the FFT algorithm really work ?
   
   The following code is commented to help you step through the code
   of the FFT.  This code is inspired by the "Numerical Recipes" discussion
   of Chapter 12.  It tries to use a similar notation.
   
   One might think that the code takes unnecessary speed penalties at different
   point in the code (check the number of points or bit reversal for instance).
   Apart from the fact that this penalty is quite small, the idea for this
   code is to be easy to read and easy to use.
   
   Daniel Cote
   dccote@novajo.ca
   
*******************************************************************************/

#include "fourier.h"    /* module's own header */

#include <iostream>

/*
 Commented FFT routine with C++'s valarray of complex<double> 

 The types:

    typedef complex<double> complex_double;
    typedef complex<float> complex_float;

 are defined in the header for readability
 
 
*/

char FFT(valarray<complex_double>& data, long nn, long isign)
{
    long    i, j, m, istep, mmax;
    double  wtemp, theta;
    complex_double w,wp,temp;

    // Verify that nn is indeed a power of two.  A power of two is 
    // represented in binary with only one bit set to 1 and all others
    // set to zero.  Hence we look for the first bit set to 1 and check
    // if this is the only one by comparing that 2 raised to the power i
    // number to nn 

    if (nn <= 2) {
        cerr << "FFT not completed: not enough points\n";
        return 1;       // error
        }

    i = 1;
    while (nn & i) {    // exit as soon as one bit is set to 1
        i <<= 1;        // go to next bit
        }
    
    if (nn != i) {
        cerr << "FFT not completed: number of points " << nn << " is not a power of 2\n";
        return 1;       // error
        }
    
    // Swapping bit-reversed elements 
    for (long i = 0 ; i < nn ; i++) {
        // Calculating index j as the bit-reversed pattern of index i
        // using 0-based indices. For example:
        //
        // 000 swapped with 000 (i.e. not swapped)
        // 001 swapped with 100
        // 010 swapped with 010 (i.e. not swapped)
        // 011 swapped with 110
        // 100 swapped with 001 (already swapped)
        // 101 swapped with 101 (i.e. not swapped)
        // 110 swapped with 011 (already swapped)
        // 111 swapped with 111 (i.e. not swapped)
        // 
        // We check for each bit in i, starting from the last one (most significant)
        // and each time a bit is set, we set the equivalent reversed bit in j
        // Since nn is passed as the total number of values in the array, we need
        // to index from 0 to nn-1 and hence we require one less bit and we need to start
        // the bit reversed mask at nn/2, not nn.
        
        int maskingbit = nn >> 1;    // Since this is a power of two, it has only one bit set
        int reversedmaskingbit = 1;  // This is the reverse bit pattern mask
        j = 0;
        
        while (maskingbit != 0) {
            if (maskingbit & i) {
                j |= reversedmaskingbit;
            }
            maskingbit             >>= 1; // Shift one bit to the right
            reversedmaskingbit  <<= 1;    // Shift one bit to the left
        }
        // Now j is the bit reversed pattern of i.

        
        // We do not want to swap the two values twice, hence we choose
        // to swap the two values when j > i
        if (j > i) {
            swap<complex_double>(data[i],data[j]);
        }

    }

    
    // Now comes the part where the actual FFT is performed.
    // We need to multiply each element by a factor w
    // (W^k in "Numerical Recipe") where w is e^(i*theta). This has more than
    // one cycle and therefore we can calculate w only once and apply it
    // to the other cycles with an extra loop with index i.
    
    // Assuming an 8-point transform, we now have an array:
    //         t0        t4        t2        t6        t1        t5        t3        t7
    // The value of k at each point in the array is:
    //         0        k        2k        3k         4k         5k         6k        7k
    // equivalent to:
    //         0        k        2k        3k        -4k        -3k        -2k        -k
    // The multiple-points transform go as follows:
    //        
    //       t0      t4      t2      t6      t1      t5      t3      t7
    //     F^eee   F^eeo   F^eoe   F^eoo   F^oee   F^oeo   F^ooe   F^ooo 
    //       \______/         \______/        \______/        \______/
    //         F^ee             F^eo            F^oe            F^oo
    //           \______________/                 \______________/
    //                  F^e                             F^o
    // where each of these is defined at all k:
    // F^ee = F^eee + e^2*pi*i/2*k * F^eeo
    // F^eo = F^eoe + e^2*pi*i/2*k * F^eoo 
    // F^oe = F^oee + e^2*pi*i/2*k * F^oeo 
    // F^oo = F^ooe + e^2*pi*i/2*k * F^ooo 
    // F^e  = F^ee  + e^2*pi*i/4*k * F^eo
    // F^o  = F^oe  + e^2*pi*i/4*k * F^oo
    // F    = F^e   + e^2*pi*i/8*k * F^o    
    //
    //       t0       t4       t2       t6       t1        t5        t3        t7
    //     F^eee    F^eeo    F^eoe    F^eoo    F^oee     F^oeo     F^ooe     F^ooo 
    //     F^ee(0)  F^ee(k)  F^eo(2k) F^eo(3k) F^oe(-4k) F^oe(-3k) F^oo(-2k) F^oo(-k)
    //     F^e(0)   F^e(k)   F^e(2k)  F^e(3k)  F^o(-4k)  F^o(-3k)  F^o(-2k)  F^o(-k)
    //     F(0)     F(k)     F(2k)    F(3k)    F(-4k)    F(-3k)    F(-2k)    F(-k)
    //
    // what is important to notice is that F^eo(0) is identical to F^eo(2k) because
    // of the period of e^i*theta.  Hence we do not need to calculate it or store it.
    // To obtain F^e(0) = F^ee(0) + w^0*F^eo(0), we use F^e(0) = F^ee(0) + w^0*F^eo(2k)
    // Similarly, for F^e(2k) = F^ee(2k) + w^2k*F^eo(2k) 
    //                        = F^ee(2k) - w^0k*F^eo(2k).
    
    // For every tranform at given k values, the corresponding w^k factors are:
    //      1       1        1        1        1         1         1         1 
    //      1       -1       1        -1       1         -1        1         -1       
    //      1       i        -1       -i       1         i         -1        -i
    //      1     Ã2+iÃ2     i      Ã2-iÃ2    -1       -Ã2-iÃ2     -i      Ã2-iÃ2  
    //
    // Since the exponential has a period equal to half the number of points in the transform
    // we can calculate e^i*theta and invert it for the second half of the cycle.

    
    // We start at the highest level of recursion (i.e. the point of 
    // two-point transforms).  mmax represents half the period of e^i*theta (always number
    // of points divided by two)
    mmax = 1;
    while ( mmax < nn) {    // Keep going until recursion is over (nn-point transform) 
        istep = mmax << 1;                        
        theta = 6.28318530717959 / (isign * istep ); // This is the argument of w = e^i*theta.
        
        
        wtemp = sin(0.5 * theta);              // We avoid calculating the same
                                               // function twice
        wp = complex_double(-2.0 * wtemp * wtemp, sin(theta));
        w  = complex_double(1.0,0);            // w (or W^k) is initially 1
        
        for (m = 0 ; m < mmax ; m++)           // We multiply each point by (W^k)^m ...
        {
            for (i = m ; i < nn ; i+=istep)    // .. in each cycle
            {
                j = i + mmax;
                temp = w * data[j];                
                data[j] = data[i] - temp;      // Multiply by e^i*(theta+pi)
                data[i] += temp;               // Multiply by e^i*theta
            }
            
            w = w * wp + w;                    // We calculate the new value of w
        }
        mmax <<= 1;                            // We go up one level, doing a transform 
                                               // with twice as many points as before.
    }

    if (isign == -1)  // If it is an inverse transform, we normalize
        data /= nn;


    return 0;         // No error
}