#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#define logb(v, b)	(log(v)/log(b))
#define SAMPLES 8
#define PI 3.1415926535897932
 
// complex datatype
typedef struct complex_t {
	double re;		// Real part
	double im;		// Imag part
} complex;
 
// function headers
complex c_from_polar(double r, double radians);
complex c_add(complex left, complex right);
complex c_sub(complex left, complex right);
complex c_mult(complex left, complex right);
double c_magnitude(complex c);
 
 
int* twiddleIndexGen(int* twindex, int log2N);
complex* twiddleGen(complex* twiddles, int* twindex, int log2N);
complex* DFT_naive(complex* x, int N);
 
// definitions
 
// convert polar value to complex
complex c_from_polar(double r, double radians){
	complex result;
	result.re = r * cos(radians);
	result.im = r * sin(radians);
	return result;
}
 
// add two complex numbers
complex c_add(complex left, complex right){
	complex result;
	result.re = left.re + right.re;
	result.im = left.re + right.re;
	return result;
}	
 
// subtract two complex numbers
complex c_sub(complex left, complex right){
	complex result;
	result.re = left.re - right.re;
	result.im = left.re - right.re;
	return result;
}
 
// multiply two complex numbers
complex c_mult(complex left, complex right){
	complex result;
	result.re = left.re * right.re - left.im * right.im;
	result.im = left.re * right.im - left.im * right.re;
	return result;
}	
 
// compute magnitude of complex number
double c_magnitude(complex c){
	return sqrt(c.re * c.re + c.im * c.im);
}
 
// compute DFT for input x
complex* DFT_naive(complex* x, int N){
	complex* X = (complex*) malloc(sizeof(struct complex_t) * N);
	complex* w = (complex*) malloc(sizeof(struct complex_t) * N);
	int k, n;
	for (k = 0; k < N; k++){
		w[k] = c_from_polar(1, -2*PI*k/N);
	}
	for (k = 0; k < N; k++){
		X[k].re = X[k].im = 0.0;
		for (n = 0; n < N; n++){
			//X[k] = c_add(X[k[, c_mult(x[n], w[(n*k)%N]));
		}
	}
	free(w);
	return X;
}
 
/*		twiddleIndexGen: 
 *	Precalculates the final index order
 *	 of the twiddle factors based on 
 *	 the number of samples.
 *	Uses the bit reversal method
 */
int* twiddleIndexGen(int* twindex, int log2N){
	int i, j;
	for (j = 0; j < (1<<log2N); j++){		// Index value
		int this_num = j;
		for (i = 0; i < log2N; i++){		// Bit		
			twindex[j] <<= 1;				
			twindex[j] += this_num & 1;		// is this bit 1?
			this_num >>= 1;					
		}
		printf("%d ", twindex[j]);
	}
	printf("\n");
	return twindex;
}
 
complex* twiddleGen(complex* twiddles, int* twindex, int log2N){
	twindex = twiddleIndexGen(twindex, log2N);
	int i, j;
	int N = (1<<log2N);
	for (i = 0; i < N; i++){
		j = twindex[i];
		twiddles[j] = c_from_polar(1, -2*PI*j/N);
		printf("%f, %f\n", twiddles[j].re, twiddles[j].im);
	}
	return twiddles;
}
 
/*		RFFT: 
 *	Recursive algorithm. Traverses through the array
 *	 of sample values in a predefined order before
 *	 traversing back and computing FFTs using twiddle
 *   factors.
 */
complex* RFFT(complex* top, complex* out, complex** left, complex** right, int l){
	int i;
	printf("\nL= %i\n", l);
	for (i = 0; i < SAMPLES/(1<<l); i++){
		left[l][i] = top[2*i];				// Fill LHS
		printf("%.0f ", left[l][i].re);
		right[l][i] = top[2*i+1];				// Fill RHS
		printf("%.0f\n", right[l][i].re);
	}
	if (1<<l>=SAMPLES){
 
	} else {
		RFFT(left[l], out, left, right, l+1);		// Traverse left
		RFFT(right[l], out, left, right, l+1);	// Traverse right
	}
	return top;
}
 
int main(int argc, char **argv) {
	int i;
	int ls = logb(SAMPLES, 2);	// # levels
	complex *main = (complex *)malloc(sizeof(complex)*SAMPLES);		// main array
	complex *out  = (complex *)malloc(sizeof(complex)*SAMPLES);		// output array
	complex **left = (complex **)malloc(sizeof(complex)*SAMPLES);	// array of left traverse arrays
	complex **right = (complex **)malloc(sizeof(complex)*SAMPLES);	// array of right traverse arrays 
	/*  (only need one for each per level since only working with a single pair at one time)
	 *  each level down the tree requires only half the number of values as the previous
	 *  so array size is halved. Total array size is roughly 4x sample size
	 */
 
	// Populate MAIN array
	for (i = 0; i < SAMPLES; i++) {
		main[i] = (complex){i,0};
		//printf("%.0f,%.0f\n\n", main[i].re, main[i].im);
	}
 
	// Generate twiddle indexes and twiddle factors
	int *twindex = (int*)malloc(sizeof(int)*SAMPLES);
	complex *twiddles = (complex*)malloc(sizeof(complex)*SAMPLES);
	twiddles = twiddleGen(twiddles, twindex, ls);
 
	// Define LEFT/RIGHT array sizes, each is 1/2 
	// as long as the previous array
	for (i = ls; i >= 0; i--){
		left[ls-i+1] = (complex*)malloc(sizeof(complex)*(1<<i));
		right[ls-i+1] = (complex*)malloc(sizeof(complex)*(1<<i));
	}
 
	//Begin Recursion
	int level = 1;
	RFFT(main, out, left, right, level);								
 
	// Cleanup
	for (i = ls; i >= 0; i--){							
		free(left[ls-i+1]);
		free(right[ls-i+1]);
	}
	free(left);
	free(right);
	free(twiddles);
	free(twindex);
	free(out);
	free(main);
	return 0;
}