Fourier transform

ftw.c

#include <fftw3.h>
#include <math.h>

#define N 256

// This thing does a fourier transform on a generated array of data
// and makes a CSV out of it so I could make pretty graphs in an
// outside tool.
//
// in the output:
//
// i = the index of the input array, you can call it t or whatever
// |inR + inI(i)| = inC
// |outR + outR(i)| = outC
//
// Compile with:
// gcc ftw.c -lfftw3 -lm
//
// You need fftw3-dev
//
// Change N above for higher resolution and slower data

void
showstats(fftw_complex val)
{
  double R, I;
  double A2, B2, C2, C;

  R = val[0];
  I = val[1];

  A2 = R * R;
  B2 = I * I;
  C2 = A2 + B2;
  C = sqrt(C2);
  printf(",%f,%f,%f", R, I, C);
}

int
main (int argc, char* argv)
{
  fftw_complex *in, *out;
  fftw_plan p;
  int i;

  in  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  p = fftw_plan_dft_1d(N, in, out, FFTW_BACKWARD, FFTW_MEASURE);

  for(i = 0; i < N; i++) {
    in[i][0] = cos(M_PI * (double)i / 5) + 
               cos(M_PI * (double)i / 17) / 11 +
               cos(M_PI * (double)i / 9) / 7;
    in[i][1] = 0; //sin(M_PI * (double)i / 4);
  }

  fftw_execute(p);

  printf("\"i\",\"inR\",\"inI\",\"inC\",\"outR\",\"outI\",\"outC\"\n");
  for(i = 0; i < N; i++) {
    printf("%d", i);
    showstats(in[i]);
    showstats(out[i]);
    printf("\n");
  }

  fftw_destroy_plan(p);

  fftw_free(out);
  fftw_free(in);

  return 0;
}