fft2.c

基2 fft算法的C程序；C程序实现。程序内部附说明

///进行DFT，基于时间抽取的
void  CMathAnalysis::Fft2_DIT(float *pr,float *pi,int n,int nflag)
{
	//float *pf_imagine,*pf_real;
	int m,j,n1,k,n2,p,q;
	float fxt,fyt,fee,fsin,fcos,fa,fxx,fyy,ftype;

	m=CalcLog2(n);
	/*pf_imagine=new float[n];
	pf_real=new float[n];
	memcpy(pf_real,pwave,n*sizeof(float));
	memset(pf_imagine,0,n*sizeof(float));*/
	//inverse bits
	j=0;
	n1=n-1;
	n2=n/2;
	for(int i=0;i<n1;++i)
	{
		if(i<j)
		{
			fxt=pr[j];
			pr[j]=pr[i];
			pr[i]=fxt;

			fxt=pi[j];
			pi[j]=pi[i];
			pi[i]=fxt;

		}
		k=n2;
		while(k<=j)
		{
			j=j-k;
			k=k/2;
		}
		j=j+k;
	}
	//calculate FFT	
	if(nflag==0)
		ftype=-2*PI;
	else
		ftype=2*PI;
	for(k=1;k<=m;++k)
	{
		n2=pow(2,k-1);
		n1=pow(2,k);
		fee=ftype/n1;
		for(j=0;j<n2;++j)
		{
			fa=j*fee;
			fcos=cos(fa);
			fsin=sin(fa);
			for(p=j;p<n;p+=n1)
			{
				q=p+n2;

				fxt=pr[p];
				fyt=pi[p];

				fxx=pr[q]*fcos-pi[q]*fsin;
				fyy=pr[q]*fsin+pi[q]*fcos;

				pr[p]=fxt+fxx;
				pi[p]=fyt+fyy;
				pr[q]=fxt-fxx;
				pi[q]=fyt-fyy;
			}
		}
	}
	if(nflag==1)
	{
		for(int loopn=0;loopn<n;++loopn)
		{	
			pr[loopn]=pr[loopn]/n;
			pi[loopn]=pi[loopn]/n;
		}

	}
	//计算模，并进行标准化
	/*k=n/2;
	fxx=2.0/n;
	for(j=0;j<k;++j)
	{
	pout[j]=fxx*sqrt(pf_real[j]*pf_real[j]+pf_imagine[j]*pf_imagine[j]);
	}

	delete pf_imagine;
	delete pf_real;*/
}