Spectral analysis of accelerometer data

In case anybody interested.

Not killing CPU, sampling at 16Hz therefore both - number of samples and sample_rate to be set to 16, sign not used here.

Hopefully didn't make any mistake :)

	function FFT(data,number_of_samples,sample_rate,sign) {
		//variables for the fft
		var n;
		var mmax;
		var m;
		var j;
		var istep;
		var i;//unsigned long
		
		//double 
		var wtemp;
		var wr;
		var wpr;
		var wpi;
		var wi;
		var theta;
		var tempr;
		var tempi;

		var vector = null;
		var swap;
		
	//the complex array is real+complex so the array 
    //as a size n = 2* number of complex samples
    //real part is the data[index] and 
    //the complex part is the data[index+1]

	//new complex array of size n=2*sample_rate
	//if(vector!=NULL)
    //    delete [] vector;

		vector = new [2*sample_rate];

	//put the real array in a complex array
	//the complex part is filled with 0's
	//the remaining vector with no data is filled with 0's
		for(n=0; n<sample_rate;n++)
		{
			if(n<number_of_samples) {
				vector[2*n]=data[n];
			} else {
				vector[2*n]=0;
			}
			vector[2*n+1]=0;
		}

	//binary inversion (note that the indexes 
    //start from 0 witch means that the
    //real part of the complex is on the even-indexes 
    //and the complex part is on the odd-indexes)
		n=sample_rate << 1;
		j=0;
		for (i=0;i<n/2;i+=2) {
			if (j > i) {
				//SWAP(vector[j],vector[i]);
				swap=vector[j];
				vector[j] = vector[i];
				vector[i] = swap;
				
				//SWAP(vector[j+1],vector[i+1]);
				swap=vector[j+1];
				vector[j+1] = vector[i+1];
				vector[i+1] = swap;
				
				
				if((j/2)<(n/4)){
					//SWAP(vector[(n-(i+2))],vector[(n-(j+2))]);
					swap=vector[(n-(i+2))];
					vector[(n-(i+2))] = vector[(n-(j+2))];
					vector[(n-(j+2))] = swap;
					
					
					//SWAP(vector[(n-(i+2))+1],vector[(n-(j+2))+1]);
					swap=vector[(n-(i+2))+1];
					vector[(n-(i+2))+1] = vector[(n-(j+2))+1];
					vector[(n-(j+2))+1] = swap;
				}
			}
			m=n >> 1;
			while (m >= 2 && j >= m) {
				j -= m;
				m >>= 1;
			}
			j += m;
		}
	//end of the bit-reversed order algorithm

	//f routine
		mmax=2;
		while (n > mmax) {
			istep=mmax << 1;
			theta=sign*(2.0f*Math.PI/((mmax)*1.0f));
			wtemp=Math.sin(0.5f*theta);
			wpr = -2.0f*wtemp*wtemp;
			wpi=Math.sin(theta);
			wr=1.0;
			wi=0.0;
			for (m=1;m<mmax;m+=2) {
				for (i=m;i<=n;i+=istep) {
					j=i+mmax;
					tempr=wr*vector[j-1]-wi*vector[j];
					tempi=wr*vector[j]+wi*vector[j-1];
					vector[j-1]=vector[i-1]-tempr;
					vector[j]=vector[i]-tempi;
					vector[i-1] += tempr;
					vector[i] += tempi;
				}
				//wr=(wtemp=wr)*wpr-wi*wpi+wr;
				wtemp=wr;
				wr=wtemp*wpr-wi*wpi+wr;
				wi=wi*wpr+wtemp*wpi+wi;
			}
			mmax=istep;
		}
	//end of the algorithm
	
	//determine the fundamental frequency
	//look for the maximum absolute value in the complex array
    
    //System.println("Out:" + vector);
		var fundamental_frequency=0;
		for(i=2; i<=sample_rate; i+=2)
		{
			if((Math.pow(vector[i],2)+Math.pow(vector[i+1],2))>(Math.pow(vector[fundamental_frequency],2)+Math.pow(vector[fundamental_frequency+1],2))){
			
				fundamental_frequency=i;
			
			}
		}
	//since the array of complex has the format [real][complex]=>[absolute value]
	//the maximum absolute value must be ajusted to half
		fundamental_frequency= (fundamental_frequency)/2;
		return fundamental_frequency;
	}

//end FFT