#include <math.h>

void fft(int nb, double *fr, double *fi, int i)
/*---------------------------------------------------------------------
 * Fast Fourier Transform and Inverse Fast Fourier Transform if the 
 * compiler does not comply with the C99 standard. 
 * If it does, try cfft instead.
 *
 * Input Parameters:
 *    nb number of data points in time domain, must be 2^n
 *    i  1 for time domain -> frequency domain (Fourier Transform)
 *       -1 for frequency domain -> time domain (inverse Four.Trans.)
 *
 * Updated parameters
 *    fr,fi real and imaginary part of complex array.
 *          in time domain:
 *             timeseries in fr, zero in fi.
 *          in frequency domain:
 *             complex Fourier coefficients in lower half of fr,fi
 *             and their conjugate complex values in inverse order 
 *             in upper half of fr,fi, i.e.:
 *             fr[0] = mean value, fi[0]=0.
 *             fr[j] = fr[nb-j] and fi[j] = -fi[nb-j] for j=1,2,...,nb/2-1
 *             fr[nb/2] = fi[nb/2] = 0.
 *
 * This corresponds to the original Cooly-Tukey algorithm.
 *
 * References:
 * "Random Vibrations and Spectral Analysis" by D.E. Newland
 *
 *---------------------------------------------------------------------
 */{
   double tmp;
   int j,jm1,k,l,lm1,lpkm1,m,n,nbd2,nbm1;
   long j_l,k_l,l_l,m_l;
   double t[2],u[2],w[2];
   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<<1))
      ;
   if((m!=nb)||(n>32))
      {
      fprintf(stderr,"fft(): 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[0]=fr[jm1];
         t[1]=fi[jm1];
         fr[jm1]=fr[lm1];
         fi[jm1]=fi[lm1];
         fr[lm1]=t[0];
         fi[lm1]=t[1];
         }
      for(k=nbd2; k<j; k=(k>>1))
         j -= k;
      j+=k;
      }

   for(m=1; m<=n; ++m)
      {
      u[0]=1.;
      u[1]=0.;
      m_l=(1L<<m);
      k_l=(m_l>>1);
      w[0]=co[m-1];
      w[1]=-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[0]=fr[lpkm1]*u[0]-fi[lpkm1]*u[1];
            t[1]=fr[lpkm1]*u[1]+fi[lpkm1]*u[0];
            fr[lpkm1]=fr[lm1]-t[0];
            fi[lpkm1]=fi[lm1]-t[1];
            fr[lm1]+=t[0];
            fi[lm1]+=t[1];
            }
         tmp=u[0];
         u[0]=u[0]*w[0]-u[1]*w[1];
         u[1]=tmp*w[1]+u[1]*w[0];
         }
      }

   if(i== 1)
      for(j=0; j<nb; ++j)
         {
         fr[j]=fr[j]/nb;
         fi[j]=fi[j]/nb;
         }

   } /* fft() */