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

#define _PI 3.14159

#define _N 100

struct ret
{
	double *val;
	unsigned int valc;
};

static double absolute(double a, double b);

struct ret dft(unsigned long N, double *val);

int main(void)
{
  double val[_N];
  for(int i = 0; i < _N ; i++)
  {
    val[i] = 0;
  }

  for(int i = 0; i < _N; i+=4)
  {
    val[i] = 1;
  }

  struct ret a = dft(_N, (double*)&val);

  for(int i = 0; i < a.valc; i++)
  {
    printf("%f,",a.val[i]);
  }

  printf("\n");

  return 0;
}

//discrete fourier transform
// N: Number of values
// *val: Array of values
struct ret dft(unsigned long N, double *val)
{
	//Dynamic memory allocation
	struct ret _ret;
	double kmax = N / 2;
	double *ival = malloc(sizeof(double) * kmax); //Array of imaginary values
	double *rval = malloc(sizeof(double) * kmax); //Array of real values

	_ret.val = malloc(sizeof(double) * kmax);
	_ret.valc = kmax;

	//Calculate Fourier-Transform for every k
	for (unsigned int k = 0; k < kmax; k++)
	{
		rval[k] = 0;
		ival[k] = 0;
		//Actual discrete Fourier-Transform
		for (unsigned int n = 0; n < N; n++)
		{
			double angle = (2 * _PI * k * n) / N;
			double rx =   val[n] * cos(angle);//real part
			double ix = - val[n] * sin(angle);//imaginary part

			rval[k] += rx;
			ival[k] += ix;
		}
		_ret.val[k] = absolute(rval[k], ival[k]);
	}

	return _ret;
}

static double absolute(double a, double b)
{
	return sqrt(pow(a, 2) + pow(b, 2));
}