fft.h

fft算法

/* ===================fft.h========================= 

   Header file for the FFT routines 

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
   
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
   
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA

   Any questions, please contact me by email:[e-mail]

================================================== */

#define uchar unsigned char
#define uint unsigned int
#define PI 3.141592


///****************related to the complex************************

//================================================
//define the struct of complex
//================================================
struct complex
{
	double real,image;
};

//================================================
//fuction:	printc(complex a)
//use:		the display fuction for complexs
//return:	0
//================================================
void printc(complex a)
{
	printf("%f%s%f",a.real," + j",a.image);
}

//================================================
//fuction:	c_add(complex a,complex b)
//use:		the add fuction for complexs
//return:	complex c
//================================================
complex c_add(complex a,complex b)
{
	complex c;
	c.real = a.real + b.real;
	c.image = a.image + b.image;
	return c;
}

//================================================
//fuction:	c_minus(complex a,complex b)
//use:		the minus fuction for complexs
//return:	complex c
//================================================
complex c_minus(complex a,complex b)
{
	complex c;
	c.real = a.real - b.real;
	c.image = a.image - b.image;
	return c;
}

//================================================
//fuction:	c_multiply(complex a,complex b)
//use:		the multiply fuction for complexs
//return:	complex c
//================================================
complex c_multiply(complex a,complex b)
{
	complex c;
	c.real = a.real*b.real - a.image*b.image;
	c.image = a.real*b.image + a.image*b.real;
	return c;
}

//================================================
//fuction:	c_mode(complex a)
//use:		get the mode value of complex
//return:	double
//================================================
double c_mode(complex a)
{
	return sqrt(pow(a.real,2) + pow(a.image,2));
}
///***************************************************


///*****************the FFT algorithms****************

/*==================the W facotr=================
there are two ways of getting the W factor:
	1.by using the formula, we can computer them once the points of FFT decided, 
	by the function named W_factor(int points),and stored into a matrix;
	2.by using a w table matrix, in which the values are fixed by the points of FFT,
	then we can get it by searching the matrix.
	here we choose the first one
==================================================*/
//================================================
//function:	w_factor(int N)
//use:		the first way of getting the w factor
//return:	the w factor matrix
//================================================
complex *w_factor(int N)
{
	complex *w_temp = new complex[N/2];
	int k;
	for(k=0;k<N/2;k++)
	{
		w_temp[k].real = cos(2*PI*k/N);
		w_temp[k].image = - sin(2*PI*k/N);
	}
	return w_temp;
}

//================================================
//function:	invert(int n,int BITS)
//use:		invert the n's binary order
//refer:	BITS is the bits of binary n
//return:	the washed n
//================================================
int invert(int n,int BITS)
{
	int temp = 0,i;
	temp = temp|(n&1);	//give temp N's lowest bit
	for(i=1;i<BITS;i++)
	{
		temp = temp<<1;
		n = n>>1;
		temp = temp|(n&1);
	}
	return temp;
}

//================================================
//function:	wash(double *time,int N,int BITS)
//use:		wash the order of the time sequence
//refer:	here N must be 2,4,8,16...
//return:	the washed time sequence
//================================================
void wash(double *time,int N,int BITS)
{
	int i = N,j;
	double temp;
	for(i=0;i<N/2;i++)		//BE CAREFUL, HERE IT MUST BE N/2, NOT N
	{
		j = invert(i,BITS);
		temp  = time[i];
		time[i] = time[j];
		time[j] = temp;
	}
}

//================================================
//function:	r2_fft(double *time,int N)
//use:		a general radix-2 FFT algorithm, 
//			the inpute must be real number
//refer:	here N must be 2,4,8,16...
//return:	the frequent matrix frequent[]
//================================================
complex *r2_fft(double *time,int N)
{
	int BITS=0,i=N,j,k,power_2;
	complex *frequent = new complex[N];
	complex *W = new complex[N/2];
	complex c_temp,c_temp_o,c_temp_e;
	//get the BITS of binary N
	while(i!=1)
	{
		i = i>>1;
		BITS ++;
	}
	W = w_factor(N);		//generate the W factors
	wash(time,N,BITS);		//wash the order of the time sequence
	//1-point DFT
	for(j=0;j<N;j++)
	{
		frequent[j].real = time[j];
		frequent[j].image = 0;
	}
	//the FFT body
	for(i=0;i<BITS;i++)
	{
		power_2=1<<i;
		for(j=0;j<N;j=j+power_2*2)
		{
			for(k=0;k<power_2;k++)
			{
				c_temp=c_multiply(frequent[k+j+power_2],W[N*k/power_2/2]);
				c_temp_o=c_add(frequent[k+j],c_temp);
				c_temp_e=c_minus(frequent[k+j],c_temp);
				frequent[k+j]=c_temp_o;
				frequent[k+j+power_2]=c_temp_e;
			}
		}
	}
	return frequent;
}
///***************************************************