www.pudn.com > howtofft_src.zip > Fourier.cpp

/********************************************************/ 
/* WARNING:												*/ 
/* This code cannot be used in any aplication			*/ 
/* without permition of the author						*/ 
/* for more information please read the license in the	*/ 
/* Numerical Recipies in C book or go to www.nr.com		*/ 
/* this is mearly an example of how to use it			*/ 
/********************************************************/ 
 
 
#include "StdAfx.h" 
#include  
#include ".\fourier.h" 
 
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr 
 
CFourier::CFourier(void) 
{ 
	pi=4*atan((double)1);vector=NULL; 
} 
 
CFourier::~CFourier(void) 
{if(vector!=NULL) 
		delete [] vector; 
} 
 
// FFT 1D 
void CFourier::ComplexFFT(float data[], unsigned long number_of_samples, unsigned int sample_rate, int sign) 
{ 
 
	//variables for the fft 
	unsigned long n,mmax,m,j,istep,i; 
	double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi; 
 
	//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 float [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 i) { 
			SWAP(vector[j],vector[i]); 
			SWAP(vector[j+1],vector[i+1]); 
			if((j/2)<(n/4)){ 
				SWAP(vector[(n-(i+2))],vector[(n-(j+2))]); 
				SWAP(vector[(n-(i+2))+1],vector[(n-(j+2))+1]); 
			} 
		} 
		m=n >> 1; 
		while (m >= 2 && j >= m) { 
			j -= m; 
			m >>= 1; 
		} 
		j += m; 
	} 
	//end of the bit-reversed order algorithm 
 
	//Danielson-Lanzcos routine 
	mmax=2; 
	while (n > mmax) { 
		istep=mmax << 1; 
		theta=sign*(2*pi/mmax); 
		wtemp=sin(0.5*theta); 
		wpr = -2.0*wtemp*wtemp; 
		wpi=sin(theta); 
		wr=1.0; 
		wi=0.0; 
		for (m=1;m(pow(vector[fundamental_frequency],2)+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=(long)floor((float)fundamental_frequency/2); 
 
}