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


double saw (double x)
{
	return (x + (x < 0.5 ? 0.5 : -0.5)) * 2.0 - 1.0;
}


/*
   Direct Fourier Transform
   dir =  1 forward 
   dir = -1 reverse 

*/
int DFT(int dir, int N, double *R, double *I) 
{
	int u, k, i;
	double p;
	double R2[N], I2[N];
        
	for (u = 0; u < N; u++) 
	{
		R2[u] = 0.0;
		I2[u] = 0.0;
		for (k = 0; k < N; k++) 
		{
			p = (-dir*2.0*M_PI*u*k)/N;
			
         		R2[u] += R[k]*cos(p) + I[k]*sin(p);
			I2[u] += I[k]*cos(p) - R[k]*sin(p);
    		}
		if (dir == 1) {
		  R2[u] /= N;
		  I2[u] /= N;
		}
    		
	} 

	for (i = 0; i < N; i++) {
        	R[i] = R2[i];
        	I[i] = I2[i];
        }

	return true;
}



int main(int argc, char **argv) {

	unsigned int i, n, steps = 16;
	double phase, step = 1.0 / steps;
        double phases[steps];
	int N = steps;
	
	double R[N], I[N];

	for (i = 0, phase = 0.0; i < steps; ++i, phase += step) 
	{
                
		R[i] = 0.0;
		I[i] = 0.0;	

		R[i] = sin(2.0*M_PI*phase);
		//R[i] = saw(phase);
		phases[i] = phase;
	        
	}

	for (i = 0; i < steps; ++i)
		printf ("%.4f\t%.8f\n", phases[i],R[i]);


	printf("\nCall dft-function forward %d\n", DFT(1, N, R, I));    

	for (i = 0; i < steps; ++i)
		printf ("%.8f\t%.8f\n", R[i],I[i]);	

	printf("\nCall dft-function reverse %d\n", DFT(-1, N, R, I));

	for (i = 0; i < steps; ++i)
		printf ("%.8f\t%.8f\n", R[i], I[i]);	  
	

	return EXIT_SUCCESS;
}