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

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

int iDivUp(int a, int b){
    return ((a % b) != 0) ? (a / b + 1) : (a / b);
}
#include "DFT_GPU_Kernel.h"

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;

	//Device Variables
	float *signal_real_dev;
	float *signal_imag_dev;
	float *DFT_real_dev;
	float *DFT_imag_dev;

	dim3 threads;
	dim3 grid;


	//Create host memory
	int N_DFT=(int)pow(2.0f,(int)(log((float)N)/log(2.0f) + 0.5f));
	size_t Total_mem=N_DFT*sizeof(float);
	signal_real=(float *)malloc(Total_mem);
	signal_imag=(float *)malloc(Total_mem);
	DFT_real=(float *)malloc(Total_mem);
	DFT_imag=(float *)malloc(Total_mem);

	//Create device memory
	cudaMalloc((void **)&signal_real_dev,Total_mem);
	cudaMalloc((void **)&signal_imag_dev,Total_mem);
	cudaMalloc((void **)&DFT_real_dev,Total_mem);
	cudaMalloc((void **)&DFT_imag_dev,Total_mem);
	

	//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 in CPU
/*	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)));///CLOCKS_PER_SEC
	*/

	//Applying the DFT in GPU
	time_ini = clock();
	//Copy Host memory to Device Memory
	cutilSafeCall(cudaMemcpy(signal_real_dev,signal_real,Total_mem,cudaMemcpyHostToDevice));
	cutilSafeCall(cudaMemcpy(signal_imag_dev,signal_imag,Total_mem,cudaMemcpyHostToDevice));
	//Define the blocks and grid
	threads.x= /*HERE*/;
	threads.y= /*HERE*/;
	threads.z= /*HERE*/;
	grid.x= /*HERE*/;
	grid.y= /*HERE*/;
	grid.z= /*HERE*/;


	DFT_GPU(DFT_real_dev,DFT_imag_dev,signal_real_dev,signal_imag_dev,N_DFT,1,grid,threads);
	//cudaDeviceSynchronize();

	//Copy device memory to host Memory
	cutilSafeCall(cudaMemcpy(DFT_real,DFT_real_dev,Total_mem,cudaMemcpyDeviceToHost));
	cutilSafeCall(cudaMemcpy(DFT_imag,DFT_imag_dev,Total_mem,cudaMemcpyDeviceToHost));
	
	cudaDeviceSynchronize();
	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){
	if((fp=fopen("../Signal&Magnitude_DFT_GPU.txt","wb"))==NULL){
		printf("\nNo se puede abrir el archivo...");
		system("pause");
	}
	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);						
		
		//fprintf(fp,"%f %f\n",signal_real[i],signal_imag[i]);						
	}

	fclose(fp);

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

	free(signal_real);
	free(signal_imag);
	free(DFT_real);
	free(DFT_imag);

	cudaFree(signal_real_dev);
	cudaFree(signal_imag_dev);
	cudaFree(DFT_real_dev);
	cudaFree(DFT_imag_dev);
	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;
      }
   } 
}


