fft.txt

快速付里叶变换子程序所需 RAM 空间以输入的首地址为基址,向增加的方向扩展

72.  函数名：_FFT

功能：FFT变换

用法：void _FFT(double *fr, double *fi, int n, int flag)

参数：fr：采样点的实数表，fi：采样点的虚数表，n：采样点个数，flag：flag=0表示求Fourier变换，flag=1表示求逆Fourier变换

返回：无

void _FFT(double *fr, double *fi, int n, int flag)
{
	int mp,arg,q,cntr,p1,p2;
	int i,j,a,b,k,m;
	double sign,pr,pi,harm,x,y,t;
	double *ca,*sa;

	ca=(double *)calloc(n,sizeof(double));
	sa=(double *)calloc(n,sizeof(double));
	if(ca==NULL||sa==NULL) Error("calloc error in _FFT!");
	j=0;
	if(flag!=0)
	{
		sign=1.0;
		for(i=0;i<=n-1;i++)
		{
			fr[i]=fr[i]/n;
			fi[i]=fi[i]/n;
		}
	}
	else
		sign=-1.0;
	for(i=0;i<=n-2;i++)
	{
		if(i<j)
		{
			t=fr[i];
			fr[i]=fr[j];
			fr[j]=t;
			t=fi[i];
			fi[i]=fi[j];
			fi[j]=t;
		}
		k=n/2;
		while(k<=j)
		{
			j-=k;
			k/=2;
		}
		j+=k;
	}
	mp=0;
	i=n;
	while(i!=1)
	{
		mp+=1;
		i/=2;
	}
	harm=2*M_PI/n;
	for(i=0;i<=n-1;i++)
	{
		sa[i]=sign*sin(harm*i);
		ca[i]=cos(harm*i);
	}
	a=2;
	b=1;
	for(cntr=1;cntr<=mp;cntr++)
	{
		p1=n/a;
		p2=0;
		for(k=0;k<=b-1;k++)
		{
			i=k;
			while(i<n)
			{
				arg=i+b;
				if(k==0)
				{
					pr=fr[arg];
					pi=fi[arg];
				}
				else
				{
					pr=fr[arg]*ca[p2]-fi[arg]*sa[p2];
					pi=fr[arg]*sa[p2]+fi[arg]*ca[p2];
				}
				fr[arg]=fr[i]-pr;
				fi[arg]=fi[i]-pi;
				fr[i]+=pr;
				fi[i]+=pi;
				i+=a;
			}
			p2+=p1;
		}
		a*=2;
		b*=2;
	}
	free(ca);
	free(sa);
}