#include "TrueBASIC.h"
int main();
void FFT(double g[][2], int p);

void FFT(double g[][2], int p)
{
    int L, N, N2, del_i, i, ip, j, jmax, k, Itmp1_;
    double Wp[2], angle, factor[2], temp[2];

    // fast Fourier transform of complex input data set g
    // transform returned in g
    N = (1 << p);                           // number of data points
    N2 = N/2;
    // rearrange input data according to bit reversal
    j = 0;
    for(i = 0; i < N-1; ++i)
    {
        // g(i,1) is real part of f((i-1)*del_t)
        // g(i,2) is imaginary part of f((i-1)*del_t)
        // set g(i,2) = 0 if data real
        if(i < j)
        {
            temp[0] = g[j][0];
            temp[1] = g[j][1];
            g[j][0] = g[i][0];
            g[j][1] = g[i][1];
            g[i][0] = temp[0];
            g[i][1] = temp[1];
        }
        k = N2;
        while(k <= j)
        {
            j -= k;
            k /= 2;
        }
        j += k;
    }
    // begin Danielson-Lanczos algorithm
    jmax = 1;
    for(L = 1; L <= p; ++L)
    {
        del_i = 2*jmax;
        Wp[0] = 1;                          //  Wp initialized at W^0
        Wp[1] = 0;
        angle = pi/jmax;
        factor[0] = cos(angle);             // ratio of new to old W^p    
        factor[1] = -sin(angle);
        for(j = 0; j < jmax; ++j)
        {
            for(i = j; i < N; i += del_i)
            {
                // calculate transforms of length 2^L
                ip = i + jmax;
                temp[0] = g[ip][0]*Wp[0] - g[ip][1]*Wp[1];
                temp[1] = g[ip][0]*Wp[1] + g[ip][1]*Wp[0];
                g[ip][0] = g[i][0] - temp[0];
                g[ip][1] = g[i][1] - temp[1];
                g[i][0] = g[i][0] + temp[0];
                g[i][1] = g[i][1] + temp[1];
            }
            // find new W^p
            temp[0] = Wp[0]*factor[0] - Wp[1]*factor[1];
            temp[1] = Wp[0]*factor[1] + Wp[1]*factor[0];
            for(Itmp1_ = 0; Itmp1_ < 2; ++Itmp1_)
                Wp[Itmp1_] = temp[Itmp1_];
        }
        jmax = del_i;
    }
}