#include "fft.h"

namespace fft{
	
int bit_swap(int i, int nu) {
	int ib, i1, i2;
	ib = 0;
	for(i1=0;i1!=nu;i1++) {
		i2 = i/2;
		ib = ib*2+i-2*i2;
		i = i2;}
	return(ib); }
 
 
void swap(double x1, double x2)  {
	double temp_x;
	temp_x = x1;
	x1 = x2;
	x2 = temp_x; }
 
 
void fft(double * real_data, double * imag_data, int n_pts, int inv) {
	int nu = (int)(log10((double)n_pts)/log10((double)2));
	int n2, j, l, i, ib, k, k1, k2;
	int sgn;
	double tr, ti, arg, nu1;
	double c, s;
	n2 = n_pts / 2;
	nu1 = nu - 1.0;
	k = 0;
	sgn = inv ? -1 : 1;
	for(l=0;l!=nu;l++) {
		do {
		  
		  for (i=0;i!=n2;i++) {
              j = (int)(k/(pow(2.0,nu1)));
			  ib = bit_swap( j, nu );
			  arg = 2.0 * PI * ib / n_pts;
			  c = cos( arg );
			  s = sgn * sin( arg );
			  k1 = k;
			  k2 = k1 + n2;
			  tr = *(real_data+k2) * c + *(imag_data+k2) * s;
			  ti = *(imag_data+k2) * c - *(real_data+k2) * s;
			  *(real_data+k2) = *(real_data+k1) - tr;
			  *(imag_data+k2) = *(imag_data+k1) - ti;
			  *(real_data+k1) = *(real_data+k1) + tr;
			  *(imag_data+k1) = *(imag_data+k1) + ti;
			  k++;
		  }
		  
		  k +=  n2; 
		  
		}while(k<n_pts-1);
		
		k = 0;
		nu1 -= 1.0;
		n2 /= 2; 
	}
	
	for (k=0;k!=n_pts;k++) {
		ib = bit_swap( k, nu );
		if(ib>k) {
		   swap(*(real_data+k), *(real_data+ib));
		   swap(*(imag_data+k), *(imag_data+ib));
		}
	}
	
	if(inv)
		for (k=0;k!=n_pts;k++) {
			*(real_data+k) /= n_pts;
			*(imag_data+k) /= n_pts; 
		}
}
};//End of namespace :-)
