/* simple Fourier transformation for regularly spaced data points */
/* input data is a pair of vectors reF[n], imF[n], representing   */
/* F(n*dx)=F(x)=reF[n]+i*imF[n] ; dx=1/N; interval is 0<=x<=1     */
/* output is a pair of vectors reG[m], imG[m] representing the    */
/* transformation. For a real function, the Fourier cosine        */
/* transform will be                                              */
/* F(n*dx)=1/(sqrt(N))*sum_m reG[m]*cos(2*pi*m*n*x/N)             */

#include<math.h>
#include<stdio.h>
#define PI 3.1415926
#define N 200         /* increase or decrease as needs be */
void FT( double *rF, double *iF, double *rg, double *ig);

double reF[N], imF[N], reG[N], imG[N];
int i,j;
double dx, x;

main(){
double *reF_ptr,*imF_ptr,*reG_ptr,*imG_ptr;

/* build a sample function; F(x)=2.0/(1.0+(x-0.5)*(x-0.5)) */
dx=1.0/(double)N;
 for(i=0;i<N;i++){
x=(double)i*dx;
reF[i]=0.1/(0.01+(x-0.5)*(x-0.5));
imF[i]=0.0;
 reG[i]=0.0; /* initialize transforms too */
 imG[i]=0.0;}

 /* point the pointers */
reF_ptr=&reF[0];
reG_ptr=&reG[0];
imF_ptr=&imF[0];
imG_ptr=&imG[0];

/* perform the transform */
FT(reF_ptr, imF_ptr, reG_ptr, imG_ptr);


/* inverse transform */
for(i=0;i<N;i++){
*(imG_ptr+i)=-(*(imG_ptr+i));      /* changes i to -i in exponential */
}
FT(reG_ptr, imG_ptr, reF_ptr, imF_ptr);

/* output the inverse Fourier transform */
/* is original function computed as a Fourier sum */
 for(i=0;i<N;i++){
printf("%f\t%f\n", (double)i/(double)N,*(reF_ptr+i));
 }

return;
}


void FT( double *rF, double *iF, double *rg, double *ig)
{
int n,m;
double arg, factor, norm,C,S;
factor=2.0*PI/(double)N;
norm=sqrt((double)N);

 for(m=0;m<N;m++){
/* zero the array sums */
*(rg+m)=0.0;
*(ig+m)=0.0;
/* now build the Fourier sums */
 for(n=0;n<N;n++){
arg=factor*((double)m)*((double)n);
C=cos(arg);
S=sin(arg);
*(rg+m)=*(rg+m)+(*(rF+n))*C+(*(iF+n))*S;
 *(ig+m)=*(ig+m)+(*(iF+n))*C-(*(rF+n))*S;}
 /* normalize the sums */
*(rg+m)=*(rg+m)/norm;
 *(ig+m)=*(ig+m)/norm;}

return ;
}