#include"math.h"

void fft(x,y,n,sign) 
/* x,y 長(zhǎng)度為n的雙精度實(shí)型一維數(shù)組。x: 開始放變換數(shù)據(jù)實(shí)部，最后放變換結(jié)果實(shí)部；
y： 開始放變換數(shù)據(jù)虛部，最后方變換結(jié)果虛部；
n： 數(shù)據(jù)長(zhǎng)度，必須為2^m； sign＝1時(shí)fft()為DFT；sign＝－1時(shí)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的數(shù)值;  MAX值可根據(jù)具體數(shù)據(jù)設(shè)置*/
	MAX=16;
	for(j=1,i=1;i<MAX;i++)
	{	m=i;
		j=2*j;
		if(j==n) break;
	}
	
	/*對(duì)輸入數(shù)據(jù)x[i]進(jìn)行倒序*/
	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;
	}
	
	/*計(jì)算各級(jí)蝶形*/
	n1=1;
	for(l=1;l<=m;l++)  /*l==L(1~M), m==M, n2==B(=2^(L-1)), 見書《數(shù)字信號(hào)處理》丁玉美 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,則每個(gè)數(shù)據(jù)都除以n*/
	if(sign==-1)     
	{	for(i=0;i<n;i++)
		{	x[i]/=n;
			y[i]/=n;
		}
	}
}