#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>

#include "fft.h"

#define  TWO_PI  6.283185307179586476925287


void complex_add(complex *result, complex *x, complex *y) {
	result->re = x->re + y->re;
	result->im = x->im + y->im;
}


void complex_sub(complex *result, complex *x, complex *y) {
	result->re = x->re - y->re;
	result->im = x->im - y->im;
}


void complex_mul(complex *result, complex *x, complex *y) {
	result->re = (x->re * y->re) - (x->im * y->im);
	result->im = (x->re * y->im) + (x->im * y->re);
}


void complex_mul_scalar(complex *result, float v) {
	result->re *= v;
	result->im *= v;
}


void complex_pow(complex *result, complex *base, int exponent) {
	complex temp;
	int i;

	result->re = 1.0;
	result->im = 0.0;
	for (i = 0; i < exponent; i++) {
		temp = *result;
		complex_mul(result, base, &temp);
	}
}


float complex_abs(complex *v) {
	return (float) sqrt((v->re * v->re) + (v->im * v->im));
}


void complex_conjugate(complex *v) {
	v->im = -(v->im);
}


void do_fft(complex *x, int ln, int inverse) {
	register int i, j, k, l, n1, n2, ie, ia, n;
	register float c, s, xt, yt, p, a;

	n = 1 << ln;
	// Quick return for n=1
	if (n == 1)
		return;
	// Conjugate if backward transform
	if (inverse) {
		for (i = 0; i < n; i++)
			x[i].im = - x[i].im;
	}
	// Main loop
	p = TWO_PI / ((float) n);
	n2 = n;
	for (k = 1; k <= ln; k++) {
		n1 = n2;
		n2 /= 2;
		ie = n / n1;
		ia = 1;
		for (j = 0; j < n2; j++) {
			a = ((float) (ia - 1)) * p;
			c = (float) cos(a);
			s = (float) sin(a);
			ia += ie;
			for (i = j; i < n; i += n1) {
				l = i + n2;
				xt = x[i].re - x[l].re;
				x[i].re += x[l].re;
				yt = x[i].im - x[l].im;
				x[i].im += x[l].im;
				x[l].re = c * xt + s * yt;
				x[l].im = c * yt - s * xt;
			}
		}
	}
	// Bit reversal permutation
	j = 0;
	for (i = 0; i < n-1; i++) {
		if (i < j) {
			xt = x[j].re;
			x[j].re = x[i].re;
			x[i].re = xt;
			yt = x[j].im;
			x[j].im = x[i].im;
			x[i].im = yt;
		}
		k = n / 2;
		while (k < (j + 1)) {
			j -= k;
			k /= 2;
		}
		j += k;
	}
	// Conjugate and normalize if backward transform
	if (inverse) {
		p = ((float) n);
		for (i = 0; i < n; i++) {
			x[i].im /= -p;
			x[i].re /= p;
		}
	}
}


void fft(complex *x, int ln) {
	do_fft(x, ln, 0);
}


void ifft(complex *x, int ln) {
	do_fft(x, ln, 1);
}


int xlog2(int n) {
	return (int) rint(log(n) / log(2.0));
}


void fill_hanning_window(float *w, int n) {
	int i;
	for (i = 0; i < n; i++)
		w[i] = 0.5 + 0.5 * cos(M_PI + ((float)i * M_PI * 2 / (n - 1)));
}