#include <math.h>
#include <complex.h>

void cfft(int nb, double complex *fx, int i)
/*---------------------------------------------------------------------
 * Fast Fourier Transform and Inverse Fast Fourier Transform if the
 * compiler complies with the C99 standard. Otherwise, take fft.c
 *
 * Input Parameters:
 *    int nb      = number of data points in time domain, must be 2^n
 *    int i       = 1 for time domain -> frequency domain (Fourier Transform)
 *                 -1 for frequency domain -> time domain (inverse Four.Trans.)
 *
 * Updated parameters
 *    double complex *fx In time domain:
 *                       timeseries in real part, zero in imaginary part.
 *                       In frequency domain: complex Fourier coefficients
 *                       in lower halfand their conjugate complex values
 *                       in inverse order in upper half, i.e.
 *                       creal(fx[0]) = mean value, cimag (fx[0])=0,
 *                       fx[j] = conjg(fx[nb-j]) corresponds to the
 *                       frequency j/duration, where j=1,2,...,nb/2-1.
 *                       fx[nb/2] = 0.
 *
 * This corresponds to the original Cooly-Tukey algorithm.
 *
 * References:
 * "Random Vibrations and Spectral Analysis" by D.E. Newland
 *
 *---------------------------------------------------------------------
 */{
   int j,jm1,k,l,lm1,lpkm1,m,n,nbd2,nbm1;
   long j_l,k_l,l_l,m_l;
   double complex t,u,w;

   static double co[]={
      -9.999999999999962E-001,
      -4.371139000183733E-008,
       7.071067657322372E-001,
       9.238795283293805E-001,
       9.807852793372718E-001,
       9.951847264044178E-001,
       9.987954561381469E-001,
       9.996988186794428E-001,
       9.999247018349539E-001,
       9.999811752815535E-001,
       9.999952938093143E-001,
       9.999988234516364E-001,
       9.999997058628658E-001,
       9.999999264657138E-001,
       9.999999816164282E-001,
       9.999999954041071E-001,
       9.999999988510268E-001,
       9.999999997127567E-001,
       9.999999999281892E-001,
       9.999999999820472E-001,
       9.999999999955118E-001,
       9.999999999988780E-001,
       9.999999999997194E-001,
       9.999999999999298E-001,
       9.999999999999825E-001,
       9.999999999999957E-001,
       9.999999999999989E-001,
       9.999999999999998E-001,
       9.999999999999999E-001,
            1.000000000000000,
            1.000000000000000,
            1.000000000000000
      };
   static double si[]={
      -8.742278000367459E-008,
       9.999999999999990E-001,
       7.071067966408575E-001,
       3.826834424611044E-001,
       1.950903273750643E-001,
       9.801714304836734E-002,
       4.906767569175357E-002,
       2.454122920569705E-002,
       1.227153862718945E-002,
       6.135884819898878E-003,
       3.067956848339383E-003,
       1.533980228971620E-003,
       7.669903400861504E-004,
       3.834951982431209E-004,
       1.917476026465663E-004,
       9.587380176390885E-005,
       4.793690093703264E-005,
       2.396845047540110E-005,
       1.198422523856115E-005,
       5.992112619388148E-006,
       2.996056309707521E-006,
       1.498028154855441E-006,
       7.490140774279308E-007,
       3.745070387139916E-007,
       1.872535193569991E-007,
       9.362675967849996E-008,
       4.681337983925003E-008,
       2.340668991962502E-008,
       1.170334495981251E-008,
       5.851672479906256E-009,
       2.925836239953128E-009,
       1.462918119976564E-009
      };

   /* n=ld(nb), exit if not integer */
   for(n=0,m=1; m<nb; ++n,m=m*2)
      ;
   if((m!=nb)||(n>32))
      {
      fprintf(stderr,"cfft(): number of data points not 2^n or too large.\n");
      exit(1);
      }

   /* start the Fourier transformation */
   nbd2=nb>>1;
   nbm1=nb-1;
   for(j=l=1; l<=nbm1; ++l)
      {
      if(l<j)
         {
         jm1=j-1;
         lm1=l-1;
         t=fx[jm1];
         fx[jm1]=fx[lm1];
         fx[lm1]=t;
         }
      for(k=nbd2; k<j; k=k>>1)
         j -= k;
      j+=k;
      }

   for(m=1; m<=n; ++m)
      {
      u=1.;
      m_l=1L<<m;
      k_l=m_l>>1;
      w=co[m-1] - i * I * si[m-1];
      for(j_l=1L; j_l<=k_l; ++j_l)
         {
         for(l_l=j_l; l_l<=(long)nb; l_l+=m_l)
            {
            lpkm1=l_l+k_l-1L;
            lm1=l_l-1L;
            t=fx[lpkm1]*u;
            fx[lpkm1]=fx[lm1]-t;
            fx[lm1]+=t;
            }
         u=u*w;
         }
      }
   if(i== 1)
      for(j=0; j<nb; ++j)
         fx[j]=fx[j]/nb;

   } /* cfft() */