dsp.h

C 语言编写的fft算法

/*/////////////////////////////////////////////////////////////////////////////////////////

--------------------------------by [e-mail]---------------------------------
							21:03 2007-10-12
/////////////////////////////////////////////////////////////////////////////////////////*/

#include <math.h>
#include <complex.h>
#include <stdlib.h>
#define	PI	3.1415926
unsigned	BinInv(unsigned target,unsigned width);
COMPLEX		*FFT(COMPLEX  *X,COMPLEX  *x, unsigned   N);
COMPLEX		*IFFT(COMPLEX  *X,COMPLEX  *x, unsigned   N);
COMPLEX		*DFT(COMPLEX  *X,COMPLEX  *x, unsigned   N);
COMPLEX		*IDFT(COMPLEX  *x,COMPLEX  *X, unsigned   N);
/*////////////////////////////////////////////////////////////////////////////////////////
/函数功能：将一个无符号整数按指定二进制位宽转化为二进制后逆序转化为十进制输出；
/入口参数：目标无符号整数；逆序位宽；
/出口参数：逆序所得无符号整数；
/注　　释：本函数用于DSP.H中FFT函数及DFT函数，并未使用其他库函数；
////////////////////////////////////////////////////////////////////////////////////////*/
unsigned	BinInv(unsigned target,unsigned width)
{
	unsigned	res=0,bit=0,i=0;

	for(i=0;i<width;i++)
	{
		bit = target % 2;		/*依次得到该整数二进制形式的低位到高位*/
		target >>= 1;		/*target /= 2*/
		res += bit * (1<<(width-1-i));	/*体现逆序所在*/
	}

	return	res;
}



/*///////////////////////////////////////////////////////////////////////////////////////
函数功能：对指定长度的采样数据进行快速傅立叶变换，并返回变换值, X = FFT(x)
入口参数：目标采样数据，保存变换值的指针，采样数据长度
出口参数：指向变换值的指针
注意事项：保存结果的空间是在函数外部申请，本函数不处理此内容
///////////////////////////////////////////////////////////////////////////////////////*/
COMPLEX	*FFT(COMPLEX  *X,COMPLEX  *x, unsigned   N)
{
	unsigned		temp=N,L=0,M=0,LE=1,LE1=0,J=0,I=0,IP=0;
	COMPLEX	U={0,0},W={0,0},T={0,0};

	while((temp>>=1)>0)	/*获得阶码M*/		
	{
		M++;
	}

	for(L=0;L<N;L++)		/*排序，得到FFT输入起始值*/
	{
		temp=BinInv(L,M);
		(X + L)->real = (x + temp)->real ;
		(X + L)->img = (x + temp)->img ;
	}
	for(L=0;L<M;L++)
	{
		LE<<=1;		/*LE=2^L,分组间隔*/
		LE1 = LE>>1;	/*LE1=LE/2，对偶跨距，实际也是一个分组中对偶点的对数*/
		U.real = 1.0;
		U.img = 1.0;
		W.real = cos(PI / LE1);
		W.img = sin(PI / LE1);	

		J = 0;
		while(J<LE1)
		{
			I = J;
			while(I < N)
			{
				IP = I + LE1;
				ComMul(X + IP,&U,&T);
				ComSub(X + I,&T,X + IP);
				ComAdd(X + I,&T,X + I);
				I += LE;
			}
			ComMul(&U,&W,&U);
			J++;
		}
	}

	return	X;
}



/*////////////////////////////////////////////////////////////////////////////////////////
函数功能：对指定长度的频域数据进行快速反傅立叶变换，并返回变换值,x = IFFT(X)
入口参数：目标频域数据，保存变换值的指针，频域数据长度
出口参数：指向变换值的指针
注意事项：保存结果的空间是在函数外部申请，本函数不处理此内容
////////////////////////////////////////////////////////////////////////////////////////*/
COMPLEX	*IFFT(COMPLEX  *x, COMPLEX  *X,unsigned   N)
{
	unsigned	temp=N,L=0,M=0,LE=1,LE1=0,J=0,I=0,IP=0;
	COMPLEX	U={0,0},W={0,0},T={0,0};

	while((temp>>=1)>0)		/*获得阶码M*/
	{
		M++;
	}

	for(L=0;L<N;L++)			/*排序，得到IFFT输入起始值*/
	{
		temp=BinInv(L,M);
		(x + L)->real = (X + temp)->real ;
		(x + L)->img = (X + temp)->img ;
	}
	for(L=0;L<M;L++)
	{
		LE<<=1;		/*LE=2^L,分组间隔*/
		LE1 = LE>>1;	/*LE1=LE/2，对偶跨距，实际也是一个分组中对偶点的对数*/
		U.real = 1.0;
		U.img = 0.0;
		W.real = cos(PI / LE1);
		W.img = sin(0-PI / LE1);	

		J = 0;
		while(J<LE1)
		{
			I = J;
			while(I < N)
			{
				IP = I + LE1;
				ComMul(x + IP,&U,&T);
				ComSub(x + I,&T,x + IP);
				ComAdd(x + I,&T,x + I);
				I += LE;
			}
			ComMul(&U,&W,&U);
			J++;
		}
	}
	for(L=0;L<N;L++)
	{
		(x + L)->real /= N;
		(x + L)->img /= N;
	}
	return	x;
}



/*///////////////////////////////////////////////////////////////////////////////////////
函数功能：
入口参数：
出口参数：
注意事项：
///////////////////////////////////////////////////////////////////////////////////////*/
COMPLEX	*DFT(COMPLEX  *X,COMPLEX  *x, unsigned   N)
{
	unsigned	n=0,k=0;
	COMPLEX	sum={0,0},temp={0,0},W={0,0};
	for(k=0;k<N;k++)
	{
		sum.real=sum.img=0;
		for(n=0;n<N;n++)
		{
			W.real=cos(2*PI*k*n/N);
			W.img=sin(-2*PI*k*n/N);
			ComMul(&W,x+n,&temp);
			ComAdd(&sum,&temp,&sum);
		}

		X[k]=sum;
	}

	return	X;
}



/*//////////////////////////////////////////////////////////////////////////////////////
函数功能：
入口参数：
出口参数：
注意事项：保存结果的空间是在函数外部申请，本函数不处理此内容
//////////////////////////////////////////////////////////////////////////////////////*/
COMPLEX	*IDFT(COMPLEX  *x,COMPLEX  *X, unsigned   N)
{
	unsigned	n=0,k=0;
	COMPLEX	sum={0,0},temp={0,0},W={0,0};
	for(n=0;n<N;n++)
	{
		sum.real=sum.img=0;
		for(k=0;k<N;k++)
		{
			W.real=cos(2*PI*k*n/N);
			W.img=sin(2*PI*k*n/N);
			ComMul(&W,X+k,&temp);
			ComAdd(&sum,&temp,&sum);
		}
		sum.real /= N;
		sum.img /= N;
		x[n]=sum;
	}

	return	x;
}