/* Discrete Fourier Transform v0.0: naive approach (left out) 
                              v1.0: sophisticated approach: decimation
			      v1.1: introduced sign S for computing the 
			            inverse transform (if S=-1)
			      v1.2: Testing the computation of the derivative
			      v1.3: (a) The function fz(double complex *z) is 
                                        given in a separate file, fz.c
                                    (b) New functions, dfz(double complex *z) 
                                        and d2fz(double complex *z) 
                                        compute the first and second 
                                                              derivatives
                                        at the points 2.0*Pi*n/N, n=0,...,N-1
				    (c) dfk(double complex *z) computes the 
                                        the Fourier coefficients of the first 
                                        derivative
			      v1.4: Experimenting with the sine transform: OK!
			            the inverse sine transform: OK!
			      v1.5: Experimenting with the cosine transform
			            The direct is OK; the inverse is OK! 
*/

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

#define N 128 // signal length: power of two: N = 2^R
#define NHALF ((N)>>1)
#define B 32 // computer wordlength
#define R 7 // exponent of signal length
#define Pi acos(-1.0L)

#define sqr(x) (x)*(x)

double ZERO, HALF;

void bit(unsigned int a){ /* représentation binaire de a sur B bits */
  int k; 
  for(k=(B-1);k>=0;k--){
    printf("%d",(a&(1<<(k)))>>(k)  );
  }
    printf("\n");
}


void bit_l(unsigned int a, int Ba){//représentation binaire de a sur le nombre de bits de a, Ba */

  int k; 
  for(k=(Ba-1);k>=0;k--){
    printf("%d",(a&(1<<(k)))>>(k)  );
  }
    printf("\n");
}


void bit_rev(unsigned int a, int Ba){// la représentation binaire ``renversée'' de a */ 
  int k; 
  for(k=0;k<Ba;k++){
    printf("%d", (a&(1<<(k)))>>(k)  );
  }
    printf("\n");
}

//this function takes as input an integer and computes the integer whose bits are exchanged across its length Ba

int bit_rev1(unsigned int a, int Ba){
    int k, k1;
    int u;
    unsigned int a1=0;
  for(k=0;k<Ba;k++){
      k1=Ba-1-k;
      u=(a&(1<<(k)))>>(k);
      a1=a1|(u<<k1);
  }
  return(a1);
}


//this function sets the input in bit-reversed order

void bit_swap(double complex *z){
    int k;
    double complex temp;
    for(k=0;k<N;k++){
	if(k<=bit_rev1(k,R)){
	temp=z[k];
	z[k]=z[bit_rev1(k,R)];
	z[bit_rev1(k,R)]=temp;
	}
    }

}

void output(double complex *z, int TAILLE){
    int n;

    for(n=0;n<TAILLE;n++){
	printf("%d\t%g\t%g\n", n, z[n]);
    }
    printf("\n");
}



void fft(double complex *z, int S){// S = -1 for the inverse
    int n,m;
    int r;
    int l1, l2, i1, i, j;
//    double Pi=acos(-1.0L);
    double complex temp;



// Assumes the input vector z[n] is in bit-reversed order!


    l2=1;
    for(r=0;r<R;r++){ //the stages of the computation

// Automating the computation for any value of N that is a power of 2

//	printf("%d %d\n",r,1<<r);

/* The parametrization of the loops over j and i is taken from Paul Bourke, 
   http://astonomy.swin.edu.au/~pbourke */ 

	l1=l2;
	l2 <<= 1;
	for(j=0;j<l1;j++){// the loop over the sites at level r
	    for(i=j;i<N;i=i+l2){// the loop over the other sites at level r
                                // they are at distance l2 from the given point

//the butterfly graph computation: at stage r we use the 2^r root of unity 
//of orders i and i1 for the elements i and i1 of the vector z 
  
		i1=i+l1;
//		printf("%d %d %d %g %g %g %g\n",r,i,i1, cexp(2.0*Pi*I*(double)i/(double)(1<<(r+1))), cexp(2.0*Pi*I*(double)i1/(double)(1<<(r+1))) );
		temp  =  z[i] + cexp(2.0*Pi*I*(double)(S*i)/(double)(1<<(r+1)))*z[i1]; 
		z[i1] = z[i] + cexp(2.0*Pi*I*(double)(S*i1)/(double)(1<<(r+1)))*z[i1]; 
		z[i]  = temp;
	    }
	}

	


    }

}


void FFT(double complex *z){// Computes the Fourier coefficients of z[n]
    bit_swap(z);
    fft(z,1);
}

void InvFFT(double complex *z){// Computes the function values corresponding to the Fourier coefficients z[n]
    int n;
    bit_swap(z);
    fft(z,-1);
    for(n=0;n<N;n++){
	z[n]=z[n]/(double)N;}
}






#include "fz.c" // the function values in the interval 0... 2*Pi

void dfz(double complex *z){// The derivative's values 
    int n;

    
    FFT(z); // computes the Fourier coefficients of z[n]

// Computes the Fourier coefficients of the derivative
    for(n=0;n<(N/2);n++){
	z[n]=-I*n*z[n];
    }

    for(n=(N/2);n<N;n++){
	z[n]=-I*(n-N)*z[n];
    }

    InvFFT(z); // computes the values of the derivative in ``real space''
               // i.e. where the function is, originally, defined
}



void dfk(double complex *z){// The derivative's values in Fourier space 
    int n;

    
    FFT(z); // computes the Fourier coefficients of z[n]

// Computes the Fourier coefficients of the derivative
    for(n=0;n<(N/2);n++){
	z[n]=-I*n*z[n];
    }

    for(n=(N/2);n<N;n++){
	z[n]=-I*(n-N)*z[n];
    }


}

void d2fz(double complex *z){

// This computes the 2nd derivative's values in ``real space''
// i.e. where the function values are, originally, defined

    int n;

    FFT(z);
    for(n=0;n<(N/2); n++){
	z[n]=-(n*n)*z[n];
    }

    for(n=(N/2);n<N;n++){
	z[n]=-(n-N)*(n-N)*z[n];
    }

    InvFFT(z);

}


void fst(double complex *z){
// Fast Sine Transform of length N/2 from FFT of length N
    int n;
    double complex ZERO=z[0];//keep the zero'th component of the signal
//Antisymmetrize around n=N/2

    for(n=0; n<(N/2); n++){
      z[n] = z[n<<1];
    }
    z[N/2] = 0.0;
    for(n=(N/2)+1;n<N;n++){
	z[n] = -z[N-n];
    }

    FFT(z);//FFT of length N

    for(n=0;n<(N/2);n++){
	z[n] = (z[n]-ZERO)/(2.0*I);
    }
}


void FST(double complex *z){
  ZERO = z[0];   // I need tp save these...
  HALF = z[N/2]; // ...in case they're not zero
  fst(z);
  //  output(z,NHALF);

}



void Invfst(double complex *z){
// Inverse Fast Sine Transform of length N/2 from FFT of length N
    int n;
    double complex ZERO=z[0];//keep the zero'th component of the signal
//Antisymmetrize around n=N/2 (A)
    z[N/2] = 0.0;
    for(n=(N/2)+1;n<N;n++){
	z[n] = -z[N-n];
    }

    FFT(z);//FFT of length N

    for(n=0;n<(N/2);n++){
	z[n] = (z[n]-ZERO)/(2.0*I)/(double)(N/4); // I need to divide by N/4
    }

    // I've already antisymmetrized in (A)--so I just put these entries in the ight slots 
    for(n=(N/2)+1;n<N;n++){
      z[n] = -z[N-n];}// OK for odd signal; for even signal z[n] = z[N-n]
}
    
void InvFST(double complex *z){
  Invfst(z);
  z[0] = ZERO;
  z[N/2] = HALF;
  output(z,NHALF);

}

// Not yet tested
/*
void fct(double complex *z){
// Fast Cosine Transform of length N/2 from FFT of length N
    int n;
    double complex ZERO = z[0]; 
//symmetrize around n=N/2;  z[N/2] doesn't belong to my signal and I take it 0
    for(n=0;n<(N/2);n++){
      z[n] = z[n<<1];
    }
    z[N/2] = 0.0; 
    for(n=(N/2)+1;n<N;n++){
	z[n] = z[N-n];
    }
    FFT(z);
    for(n=0;n<(N/2);n++){
	z[n] = (z[n]  + ZERO)/2.0;}


}




void Invfct(double complex *z){
// Fast Inverse Cosine Transform of length N/2 from FFT of length N
    int n;
    double complex ZERO = z[0]; 
    double complex HALF = z[N/2]; // the N/2 Fourier component is NOT always 0!

    fct(z);

    for(n=0;n<(N/2);n++){
	z[n] = -ZERO + HALF*(1-2*(n%2)) + 2.0*z[n]; }

// The zero'th component doesn't have the same normalization as the others!

    z[0] = z[0]/(double)N;
    for(n=1;n<(N/2);n++){
	z[n] = 2.0*z[n]/(double)N;}

// since the signal is symmetric about n=N/2 I have 

    for(n=(N/2)+1;n<N;n++){
	z[n] = z[N-n];}


}
*/








main(){
  double complex z[N],w[N];
//    double Pi=acos(-1.0L);
    int n;

    fz(z);

    //    output(z, N);


    FST(z); 

    InvFST(z);


	
    

  
       
    //    FFT(z);
    //    output(z, N);

    //        InvFFT(z);
    //	output(z, N);


    /*  for(n=0;n<N;n++){
      printf("%d\t%g\n",n,sqr(cabs(z[n])) );
    }
    */

}