# include "data.h"

#define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void fft (float *Data, int N, int isign) {

  int n,mmax,m,j,istep,i;
  double wtemp,wr,wpr,wpi,wi,theta;
  float tempr,tempi, *data;

  data = Data;
  n = N << 1;
  j = 0;

  for (i = 0; i < n; i+=2) {
    if (j > i) {
      FSWAP (data[j], data[i]);
      FSWAP (data[j+1], data[i+1]);
    }
    m = n >> 1;
    while (m >= 2 && j >= m) {
      j -= m;
      m >>= 1;
    }
    j += m;
  }
  mmax = 2;
  while (n > mmax) {
    istep = 2*mmax;
    theta = 6.28318530717959 / (isign*mmax);
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0;
    wi = 0.0;
    for (m = 0; m < mmax; m+=2) {
      for (i = m; i < n; i+=istep) {
	j = i + mmax;
	tempr = wr*data[j] - wi*data[j+1];
	tempi = wr*data[j+1] + wi*data[j];
	data[j] = data[i] - tempr;
	data[j+1] = data[i+1] - tempi;
	data[i] += tempr;
	data[i+1] += tempi;
      }
      wr = (wtemp = wr)*wpr - wi*wpi+wr;
      wi = wi*wpr + wtemp*wpi + wi;
    }
    mmax = istep;
  }
}

void fftold (float *Data, int N, int isign) {

  int n,mmax,m,j,istep,i;
  double wtemp,wr,wpr,wpi,wi,theta;
  float tempr,tempi, *data;

  data = Data - 1;
  n = N << 1;
  j = 1;

  for (i = 1; i < n; i+=2) {
    if (j > i) {
      FSWAP (data[j], data[i]);
      FSWAP (data[j+1], data[i+1]);
    }
    m = n >> 1;
    while (m >= 2 && j > m) {
      j -= m;
      m >>= 1;
    }
    j += m;
  }
  mmax = 2;
  while (n > mmax) {
    istep = 2*mmax;
    theta = 6.28318530717959 / (isign*mmax);
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0;
    wi = 0.0;
    for (m = 1; m < mmax; m+=2) {
      for (i = m; i <= n; i+=istep) {
	j = i + mmax;
	tempr = wr*data[j] - wi*data[j+1];
	tempi = wr*data[j+1] + wi*data[j];
	data[j] = data[i] - tempr;
	data[j+1] = data[i+1] - tempi;
	data[i] += tempr;
	data[i+1] += tempi;
      }
      wr = (wtemp = wr)*wpr - wi*wpi+wr;
      wi = wi*wpr + wtemp*wpi + wi;
    }
    mmax = istep;
  }
}

/* convert indices to zero reference */
void fftN (float *data, int *nn, int ndim, int isign) {

  int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
  int ibit,idim,k1,k2,n,nprev,nrem,ntot;
  float tempi,tempr;
  double theta,wi,wpi,wpr,wr,wtemp;

  ntot = 1;
  for (idim = 0; idim < ndim; idim++) ntot *= nn[idim];

  nprev = 1;
  for (idim = ndim - 1; idim >= 0; idim--) {
    n  =  nn[idim];
    nrem = ntot / (n*nprev);
    ip1 = nprev << 1;
    ip2 = ip1*n;
    ip3 = ip2*nrem;
    i2rev = 0;
    for (i2 = 0; i2 < ip2; i2+=ip1) {
      if (i2 < i2rev) {
	for (i1 = i2; i1 <= i2+ip1-2; i1+=2) {
	  for (i3 = i1; i3 < ip3; i3+=ip2) {	
	    i3rev = i2rev+i3-i2;
	    FSWAP(data[i3],data[i3rev]);
	    FSWAP(data[i3+1],data[i3rev+1]);
	  }
	}
      }
      ibit = ip2 >> 1;
      while (ibit >= ip1 && i2rev >= ibit) {
	i2rev -= ibit;
	ibit >>= 1;
      }
      i2rev += ibit;
    }
    ifp1 = ip1;
    while (ifp1 < ip2) {
      ifp2 = ifp1 << 1;
      theta = isign*6.28318530717959/(ifp2/ip1);
      wtemp = sin(0.5*theta);
      wpr = -2.0*wtemp*wtemp;
      wpi = sin(theta);
      wr = 1.0;
      wi = 0.0;
      for (i3 = 0; i3 < ifp1; i3+=ip1) {
	for (i1 = i3; i1 <= i3+ip1-2; i1+=2) {
	  for (i2 = i1; i2 < ip3; i2+=ifp2) {
	    k1 = i2;
	    k2 = k1+ifp1;
	    tempr = wr*data[k2]-wi*data[k2+1];
	    tempi = wr*data[k2+1]+wi*data[k2];
	    data[k2] = data[k1]-tempr;
	    data[k2+1] = data[k1+1]-tempi;
	    data[k1] += tempr;
	    data[k1+1] += tempi;
	  }
	}
	wr = (wtemp = wr)*wpr-wi*wpi+wr;
	wi = wi*wpr+wtemp*wpi+wi;
      }
      ifp1 = ifp2;
    }
    nprev *= n;
  }
}

#undef FSWAP

/* based on the PRESS routine, this fft takes an array from data[0] to data[2N-1] */
/* this function takes Data = h(t) and replaces it in situ with H(F) or vice versa.
   There are assumed to be 2*N input values with 
   Data[0,2,4,...] the real and Data[1,3,5,...] the imaginary ones.
   the output is ordered the same.  

   for h(t), values are in time sequence order.
   for H(F), values are in order F = 0, 1/N, ... 1/2 - 1/N, +/- 1/2, -1/2 + 1/N, ... -1/N 

   no normalization is performed, so a signal of amplitude A sin (w_k * t) will be 
      give an H(F) value of 0.5 * A * N, and the DC term will have H(0) = A * N.

*/ 

int IsBinary (int N) {

  int i, nbit;

  /* check if number is a binary number */
  nbit = 0;
  for (i = 0; i < 8*sizeof(N); i++) {
    nbit += (N & 0x01);
    N = (N >> 1);
  }
  if (nbit == 1) { 
    return (1); 
  } else { 
    return (0); 
  }

}  

#define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void fourn (float *data, int *nn, int ndim, int isign) {

  int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
  int ibit,idim,k1,k2,n,nprev,nrem,ntot;
  float tempi,tempr;
  double theta,wi,wpi,wpr,wr,wtemp;

  ntot = 1;
  for (idim = 1; idim <= ndim; idim++)
    ntot *= nn[idim];
  nprev = 1;
  for (idim = ndim; idim >= 1; idim--) {
    n = nn[idim];
    nrem = ntot/(n*nprev);
    ip1 = nprev << 1;
    ip2 = ip1*n;
    ip3 = ip2*nrem;
    i2rev = 1;
    for (i2 = 1; i2 <= ip2; i2+=ip1) {
      if (i2 < i2rev) {
	for (i1 = i2; i1 <= i2+ip1-2; i1+=2) {
	  for (i3 = i1; i3 <= ip3;i3+=ip2) {
	    i3rev = i2rev+i3-i2;
	    FSWAP(data[i3],data[i3rev]);
	    FSWAP(data[i3+1],data[i3rev+1]);
	  }
	}
      }
      ibit = ip2 >> 1;
      while (ibit >= ip1 && i2rev > ibit) {
	i2rev -=  ibit;
	ibit >>= 1;
      }
      i2rev += ibit;
    }
    ifp1 = ip1;
    while (ifp1 < ip2) {
      ifp2 = ifp1 << 1;
      theta = isign*6.28318530717959/(ifp2/ip1);
      wtemp = sin(0.5*theta);
      wpr  =  -2.0*wtemp*wtemp;
      wpi = sin(theta);
      wr = 1.0;
      wi = 0.0;
      for (i3 = 1;i3<=ifp1;i3+=ip1) {
	for (i1 = i3;i1<=i3+ip1-2;i1+=2) {
	  for (i2 = i1;i2<=ip3;i2+=ifp2) {
	    k1 = i2;
	    k2 = k1+ifp1;
	    tempr = wr*data[k2]-wi*data[k2+1];
	    tempi = wr*data[k2+1]+wi*data[k2];
	    data[k2] = data[k1]-tempr;
	    data[k2+1] = data[k1+1]-tempi;
	    data[k1] += tempr;
	    data[k1+1] += tempi;
	  }
	}
	wr = (wtemp = wr)*wpr-wi*wpi+wr;
	wi = wi*wpr+wtemp*wpi+wi;
      }
      ifp1 = ifp2;
    }
    nprev *= n;
  }
}

#undef FSWAP
