rfft.c

来自「2048点实序列fft 采用优化地方法节省存储空间和运算时间。包括2048点实序」· C语言 代码 · 共 117 行

//文件名:RFFT.c
#include "config.h"
//实数列基2蝶形快速FFT 

/*实型一维数组，长度为n，开始存放要交换数据的实部，最后存放变换结果的
前n/2+1个值，其存储顺序为[re(0),re(1),...,re(n/2),im(n/2-1),...,im(1)] 实序列im(0)=0*/
//data,长度为N的一维数组，是入口，也是出口
void rFFT(float *data,int n) 
{ 
	int i,j,k,m,i1,i2,i3,i4,n1,n2,n4; 
	float a,e,cc,ss,xt,t1,t2; 
	for(j=1,i=1;i<16;i++)	//the max number of transformed data is 2^16 
	{ 
		m=i; 
		j*=2; 
		if (j==n) 
		  	break; 
	} 
	n1=n-1; 
	for(j=0,i=0;i<n1;i++) 
	{ 
		if(i<j) 
		{ 
			xt=data[i]; 
			data[i]=data[j]; 
			data[j]=xt; 
		} 
		k=n/2; 
		while(k<(j+1)) 
		{ 
			j-=k; 
			k/=2; 
		} 
		j+=k; 
	} 
	for(i=0;i<n;i+=2) 
	{ 
		xt=data[i]; 
		data[i]=xt+data[i+1]; 
		data[i+1]=xt-data[i+1]; 
	} 
	n2=1; 
	for(k=2;k<=m;k++) 
	{ 
		n4=n2; 
		n2=2*n4; 
		n1=2*n2; 
		e=4*asin(1)/n1; 
		for(i=0;i<n;i+=n1) 
		{ 
			xt=data[i]; 
			data[i]=xt+data[i+n2]; 
			data[i+n2]=xt-data[i+n2]; 
			data[i+n2+n4]=(-1)*data[i+n2+n4]; 
			a=e; 
			for(j=1;j<=(n4-1);j++) 
			{ 
				i1=i+j; 
				i2=i-j+n2; 
				i3=i+j+n2; 
				i4=i-j+n1; 
				cc=cos(a); 
				ss=sin(a); 
				a=a+e; 
				t1=cc*data[i3]+ss*data[i4]; 
				t2=ss*data[i3]-cc*data[i4]; 
				data[i4]=data[i2]-t2; 
				data[i3]=(-1)*data[i2]-t2; 
				data[i2]=data[i1]-t1; 
				data[i1]=data[i1]+t1; 
			} 
		} 
	} 
}


/****************************************************
排序算法，对tab_mod[]中的数据进行从大到小的排序
***************************************************/
void sort(float* mod,float* count_in)
{
	float temp;
	uint i,j;
	for(i=0;i<=Num/2;i++)
	  	count_in[i]=i;
	for(i=0;i<=Num/2;i++)
	{
	  	for(j=i;j<=Num/2;j++)
		{
			if(mod[i]<mod[j])
			{
				temp=mod[i];
				mod[i]=mod[j];
				mod[j]=temp;
				temp=count_in[i];
				count_in[i]=count_in[j];
				count_in[j]=(unsigned char)(temp);
			}
		}
	}
}


void rAmplitude(void)
{
  	uint i;
	for(i=1;i<Num/2;i++)
	  	rRe[i]=(sqrt(rRe[i]*rRe[i]+rRe[Num-i]*rRe[Num-i]));
	if(rRe[0]<0)
	  	rRe[0]=-rRe[0]/2;
	else  rRe[0]=rRe[0]/2; 	
	if(rRe[Num/2]<0)
	  	rRe[Num/2]=-rRe[Num/2];
	/*for(i=1;i<Num/2;i++)
	 	rRe[i]=(uint)(rRe[i]/rRe[0]*2500);		//范围校正
	rRe[0]=2500;*/	  
}