/*                               -*- Mode: C -*- 
 * Filename: dft.c
 * Copyright (C) Dan Noland 2003
 * Author: Dan Noland
 * Created: Thu Mar 11 03:19:19 2004
 *           By: Dan Noland
 * Last-Updated: Thu Mar 11 03:51:54 2004
 *     Update #: 12
 * Status: 
 */

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

#define TRUE 1
#define FALSE 0

// Take, eat, this is my code which is hacked for you...
// (In case you didn't grasp that one consider this public domain)
// nolandda Thu Mar 11 03:49:18 EST 2004

/*
  Discrete Fourier Transform
*/

int dft(long int length, double real_sample[], double imag_sample[])
{
  long int i, j;
  double arg;
  double cosarg,sinarg;
  double *temp_real=NULL,*temp_imag=NULL;

  temp_real = calloc(length, sizeof(double));
  temp_imag = calloc(length, sizeof(double));
  if (temp_real == NULL || temp_imag == NULL)
  {
    return(FALSE);
  }

  for(i=0; i<length; i+=1) 
  {
    temp_real[i] = 0;
    temp_imag[i] = 0;
    arg = -1.0 * 2.0 * 3.141592654 * (double)i / (double)length;
    for(j=0; j<length; j+=1) 
    {
      cosarg = cos(j * arg);
      sinarg = sin(j * arg);
      temp_real[i] += (real_sample[j] * cosarg - imag_sample[j] * sinarg);
      temp_imag[i] += (real_sample[j] * sinarg + imag_sample[j] * cosarg);
    }
  }

  /* Copy the data back */
  for (i=0; i<length; i+=1) 
  {
    real_sample[i] = temp_real[i];
    imag_sample[i] = temp_imag[i];
  }

  free(temp_real);
  free(temp_imag);
  return(TRUE);
}


/*
  Inverse Discrete Fourier Transform
*/

int inverse_dft(long int length, double real_sample[], double imag_sample[])
{
  long int i, j;
  double arg;
  double cosarg,sinarg;
  double *temp_real=NULL,*temp_imag=NULL;

  temp_real = calloc(length, sizeof(double));
  temp_imag = calloc(length, sizeof(double));
  if (temp_real == NULL || temp_imag == NULL)
  {
    return(FALSE);
  }

  for(i=0; i<length; i+=1) 
  {
    temp_real[i] = 0;
    temp_imag[i] = 0;
    arg = 2.0 * 3.141592654 * (double)i / (double)length;
    for(j=0; j<length; j+=1) 
    {
      cosarg = cos(j * arg);
      sinarg = sin(j * arg);
      temp_real[i] += (real_sample[j] * cosarg - imag_sample[j] * sinarg);
      temp_imag[i] += (real_sample[j] * sinarg + imag_sample[j] * cosarg);
    }
  }

  /* Copy the data back */
  for (i=0; i<length; i+=1) 
  {
    real_sample[i] = temp_real[i] / (double)length;
    imag_sample[i] = temp_imag[i] / (double)length;
  }

  free(temp_real);
  free(temp_imag);
  return(TRUE);
}

int main()
{
  double* dr = calloc(100, sizeof(double));
  double* di = calloc(100, sizeof(double));
  dft(100, dr, di);
  inverse_dft(100, dr, di);
  return(1);
}