// DFT_CPU.cpp : Defines the entry point for the console application.
//

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

#define PI 3.1416

void DFT(float *Dst_real,float *Dst_imag,float *signal_real,float *signal_imag,int N_DFT,int dir);

int main( int argc, char** argv ){
	int N=8000;
	float A=0.5f; //Amplitude
	float f=10.0f;//Frequency

	int i;
	float *signal_real;
	float *signal_imag;

	float *DFT_real;
	float *DFT_imag;

	float MagDFT;
	int time_ini;


	//Create host memory
	int N_DFT=(int)pow(2.0f,(int)(log((float)N)/log(2.0f) + 0.5f));

	signal_real=(float *)malloc(N_DFT*sizeof(float));
	signal_imag=(float *)malloc(N_DFT*sizeof(float));

	DFT_real=(float *)malloc(N_DFT*sizeof(float));
	DFT_imag=(float *)malloc(N_DFT*sizeof(float));


	//Creating the signal
	for(i=0;i<N_DFT;i++){
		signal_real[i]=A*sinf(2*PI*f*i/((float)N_DFT));
		signal_imag[i]=0.0f;
	}

	//Applying the DFT
	time_ini = clock();
	DFT(DFT_real,DFT_imag,signal_real,signal_imag,N_DFT,1);//1-->Forward, -1-->Inverse
	printf ("Time clock: (%f mseconds).\n",(float)(clock()-time_ini)/1000.0f);


	//Applying the inverse DFT
	//DFT(signal_real,signal_imag,DFT_real,DFT_imag,N_DFT,-1);//1-->Forward, -1-->Inverse

	//Saving the signal and the Magnitude of DFT
	FILE *fp;
	if((fp=fopen("../Signal&Magnitude_DFT_CPU.txt","wb"))==NULL){
		printf("\nNo se puede abrir el archivo...");
	}
	fprintf(fp,"%d\n",N_DFT);
	for(i=0;i<N_DFT;i++){
		MagDFT=sqrtf(DFT_real[i]*DFT_real[i] + DFT_imag[i]*DFT_imag[i]);
		fprintf(fp,"%f %f\n",signal_real[i],MagDFT);												
	}

	fclose(fp);

	printf("\n\nReady...\n\n");

	free(signal_real);
	free(signal_imag);
	free(DFT_real);
	free(DFT_imag);
	return(0);
}

/*
Discrete Fourier Transform
dir=1 forward transform
dir=-1 inverse transform
*/
void DFT(float *Dst_real,float *Dst_imag,float *signal_real,float *signal_imag,int N_DFT,int dir){
   int i,k;
   float theta;
   float cos_theta,sin_theta;
   
   for (i=0;i<N_DFT;i++){
      Dst_real[i] = 0.0f;
      Dst_imag[i] = 0.0f;
      theta = dir * 2.0f * PI * i / (float)N_DFT;
      for (k=0;k<N_DFT;k++){
         cos_theta = cosf(k * theta);
         sin_theta = sinf(k * theta);
         Dst_real[i] += (signal_real[k] * cos_theta + signal_imag[k] * sin_theta);
         Dst_imag[i] += (signal_imag[k] * cos_theta - signal_real[k] * sin_theta);
      }
   }
   if (dir == 1) {
      for (i=0;i<N_DFT;i++) {
         Dst_real[i]/=(float)N_DFT;
         Dst_imag[i]/=(float)N_DFT;
      }
   } 
}
