fft.c

自己编的数值计算方法的c语言实现源程序。

#include"math.h"

void fft(x,y,n,sign) 
/* x,y 长度为n的双精度实型一维数组。x: 开始放变换数据实部，最后放变换结果实部；
y： 开始放变换数据虚部，最后方变换结果虚部；
n： 数据长度，必须为2^m； sign＝1时fft()为DFT；sign＝－1时fft()为IDFT; */
int n,sign;
double x[],y[];
{	int i,j,k,l,m,n1,n2,MAX;
	double c,cl,e,s,sl,t,tr,ti;
	
	/*检查n是否为2^m，得到m的数值;  MAX值可根据具体数据设置*/
	MAX=16;
	for(j=1,i=1;i<MAX;i++)
	{	m=i;
		j=2*j;
		if(j==n) break;
	}
	
	/*对输入数据x[i]进行倒序*/
	n1=n-1;
	for(j=0,i=0;i<n1;i++)
	{	if(i<j)
		{	tr=x[j];
			ti=y[j];
			x[j]=x[i];
			y[j]=y[i];
			x[i]=tr;
			y[i]=ti;
		}
		k=n/2;
		while(k<(j+1))
		{	j=j-k;
			k=k/2;
		}
		j=j+k;
	}
	
	/*计算各级蝶形*/
	n1=1;
	for(l=1;l<=m;l++)  /*l==L(1~M), m==M, n2==B(=2^(L-1)), 见书《数字信号处理》丁玉美 page103*/
	{	n1=2*n1;
		n2=n1/2;
		e=3.14159265359/n2;
		c=1.0;
		s=0.0;
		cl=cos(e);
		sl=-sign*sin(e);
		for(j=0;j<n2;j++)
		{	for(i=j;i<n;i+=n1)
			{	k=i+n2;
				tr=c*x[k]-s*y[k];
				ti=c*y[k]+s*x[k];
				x[k]=x[i]-tr;
				y[k]=y[i]-ti;
				x[i]=x[i]+tr;
				y[i]=y[i]+ti;
			}
			t=c;
			c=c*cl-s*sl; /*cos(a+e)=cos(a)*cos(e)-sin(a)*sin(e)*/
			s=t*sl+s*cl; /*sin(a+e)=cos(a)*sin(e)+sin(a)*cos(e)*/
		}
	}
	
	/*若为IDFT,则每个数据都除以n*/
	if(sign==-1)     
	{	for(i=0;i<n;i++)
		{	x[i]/=n;
			y[i]/=n;
		}
	}
}