/*

  Fast Fourier Transform

  Probably the most basic type.  Being unable to comprehend anyone's
  description of the algorithm for years, it occured to me the other day that
  it can't be anything other than optimization of the non-fast version.  So
  I wrote out all of the operations, looked for obvious optimizations, and
  made them.  Probably the code could be optimized a lot more, but nothing
  like the drastic change between the slow fourier transform algorithm which
  is O(n*n) vs. one like this which is O(n*log(n)), and thus I don't care.
  Sadly, I don't think this is any easier to understand than anyone else's
  description of the algorigthm, but if you're trying to understand how it
  works, I left some printf()s in there you can uncomment to see a list of
  all of the operations it performs.  Probably not much help, though.

  First call fft_initalize(), first argument is the number of samples, which
  must be a power of 2, second argument is 0 for a normal transform (to convert
  sample data to frequency data) or 1 for a reverse transform (to convert
  frequency data to sample data).  Then call fft_fast(), first argument is
  a complex * for the output data, second argument is a complex * for the
  input data.

  If your input data is real data (which is probably the case if you have no
  idea what I'm talking about) then just create an array of the complex type
  (which is a standard C99 data type) and copy the values to that array, which
  will set the imaginary components to zero automatically.

  If you just want the frequency amplitudes from the return data, then just
  take cabs(output[i]) for each item in the output array, which converts the
  complex number to a real number that represents the magnitude of that
  frequency in the input data.

  Compile with "-std=gnu99" and "-lm" options.  Also, don't forget to use
  "-O3" and "--fast-math" when compiling any code you want to run fast.

  Released under the antivirial license.  Basically, you can do anything
  you want with it as long as what you want doesn't involve the GNU GPL.
  See http://www.ecstaticlyrics.com/antiviral/ for more information.

*/

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

static int samples = 0;
static int inverse = 0;
static complex *table = NULL;
static complex *data = NULL;
static complex *temp = NULL;

void fft_initialize (int s, int i) {
  if (s & (s - 1)) {
    fprintf(stderr, "FFT: Sample size is not a power of two.\n");
    exit(1);
  };
  if (s == samples && !inverse == !i) return;
  samples = s; inverse = i;
  table = realloc(table, samples * sizeof(complex));
  if (table == NULL) {
    fprintf(stderr, "FFT: Failed to allocate memory for table.\n");
    exit(1);
  };
  data = realloc(data, samples * sizeof(complex));
  if (data == NULL) {
    fprintf(stderr, "FFT: Failed to allocate memory for data table.\n");
    exit(1);
  };
  temp = realloc(temp, samples * sizeof(complex));
  if (temp == NULL) {
    fprintf(stderr, "FFT: Failed to allocate memory for temp table.\n");
    exit(1);
  };
  for (int i = 0; i < samples; i++) {
    table[i] = cexp((inverse ? +1 : -1) * 2 * M_PI * I * i / samples);
  };
};

static void recursive (complex *data, int offset, int step, int count, complex *temp) {
  if (count > 2) {
    recursive(temp, offset, step << 1, count >> 1, data);
    recursive(temp, offset + step, step << 1, count >> 1, data);
  };
  //printf("offset=%d count=%d step=%d:\n", offset, count, step);
  for (int i = 0; i < count >> 1; i++) {
    int a = offset + step * i;
    int b = offset + step * (2 * i);
    int c = offset + step * (2 * i + 1);
    int d = (step * i) & (samples - 1);
    data[a] = temp[b] + temp[c] * table[d];
    //printf("  o[%d] = i[%d] + i[%d] * a(%d)\n", a, b, c, d);
    //printf("    (%+0.2f, %+0.2f) = (%+0.2f, %+0.2f) + (%+0.2f, %+0.2f) * (%+0.2f, %+0.2f)\n", creal(data[a]), cimag(data[a]), creal(temp[b]), cimag(temp[b]), creal(temp[c]), cimag(temp[c]), creal(table[d]), cimag(table[d]));
  };
  for (int i = 0; i < count >> 1; i++) {
    int a = offset + step * (i + (count >> 1));
    int b = offset + step * (2 * i);
    int c = offset + step * (2 * i + 1);
    int d = (step * (i + (count >> 1))) & (samples - 1);
    data[a] = temp[b] + temp[c] * table[d];
    //printf("  o[%d] = i[%d] + i[%d] * a(%d)\n", a, b, c, d);
    //printf("    (%+0.2f, %+0.2f) = (%+0.2f, %+0.2f) + (%+0.2f, %+0.2f) * (%+0.2f, %+0.2f)\n", creal(data[a]), cimag(data[a]), creal(temp[b]), cimag(temp[b]), creal(temp[c]), cimag(temp[c]), creal(table[d]), cimag(table[d]));
  };
};

void fft_fast (complex *output, complex *input) {
  if (samples) {
    if ((int) trunc(log(samples) / log(2) + 0.5) & 1) {
      memmove(temp, input, samples * sizeof(complex));
    } else {
      memmove(data, input, samples * sizeof(complex));
    };
    recursive(data, 0, 1, samples, temp);
    for (int i = 0; i < samples; i++) {
      if (inverse) {
        output[i] = data[i];
      } else {
        output[i] = data[i] / samples;
      };
    };
  } else {
    fprintf(stderr, "FFT: fft_initialize() has not been called.\n");
    exit(1);
  };
};