#include <stdio.h>
#include <math.h>
#include "fft1d.h"

void scale(int n, Complex *x, double c)
{
  int i;
  for(i=0;i<n;i++){
    x[i].re *= c;
    x[i].im *= c;      
  }
}

/* 
   Initial permutation
   a[x] <-> a[y], where bit_i (x) = bit_n-i (y)
   ex: a[001] <-> a[100] for n = 8
*/

void initPerm (int n, Complex *x)
{
  int i,j,m;
  Complex t;

  j = n/2;
  for (i=1; i<n; i++) { 
    if (j > i) {
      // printf("data[%d] <-> data[%d]\n", i, j);
      t = x[i];
      x[i] = x[j];
      x[j] = t;
      }
    m = n/2;
    while (m >= 1 && j >= m) {
      j -= m;
      m = m/2;
    }
    j += m;
  }
}

/*
  w(k) = cos(k 2pi/n) + i sin(k 2pi/n)
  Alternative impl : computing w on the fly propagates errors.
*/

void computeW (int n, Complex *w)
{
  int i;
  for (i=0; i<n/2; i++) {
    w[i].re = cos(2*M_PI*i/n);
    w[i].im = sin(2*M_PI*i/n);
  }
}

/*
  butterfly (a, b, w)
  {
    newa = a+b*w;
    newb = a-b*w;
    a = newa;
    b = newb;
  }  
*/

void butterfly (Complex &a, Complex &b, Complex &w)
{
  double tr, ti;

  tr = b.re * w.re - b.im * w.im;
  ti = b.re * w.im + b.im * w.re;

  b.re = a.re - tr; 
  b.im = a.im - ti;
  a.re += tr; 
  a.im += ti;
}

void fftCore (int n, Complex *x, Complex *w)
{
  int i,j, step, iw, diw;

  // diw = n / 2*step
  // n log n levels
  for (step=1, diw=n/2; step < n; step*=2, diw/=2) {

    for (i=0; i<n; i+=2*step) {
      iw = 0;

      for (j=0; j<step; j++) {
	butterfly (x[i+j], x[i+j+step], w[iw]);
	iw += diw;
      }
    }
  }
}

void fft (int n, Complex *x, Complex *w)
{
  initPerm (n, x);
  fftCore (n, x, w);
  scale (n, x, 1.0/sqrt(n));
}

void ifft (int n, Complex *x, Complex *w)
{
  int i;
  for (i=0; i<n; i++)
    x[i].im = -x[i].im;

  fft(n, x, w);
}