fft.c

实现对信号进行傅立叶变化的程序

#define Pi 3.141592653589793238462643

void FFT_ArraySetValue(double *s1, double *s2, int numdat)
{
	int i;
	for(i=0;i<numdat;i++){s2[i]=s1[i];}
}

void FFT_swap(double *s1, double *s2)
{
	double temp;
	temp=(*s1);
	(*s1)=(*s2);
	(*s2)=temp;
}

double windfunc(double t, double h,int numdat,int PP)
{
	double y;
	switch(PP)
	{
	case 1: 
		y=(1+cos(Pi*t/(double)numdat/h));
		break;
	case 2: 
		y=2./3.*(1+cos(Pi*t/(double)numdat/h))*(1+cos(Pi*t/(double)numdat/h));
		break;
	case 3: 
		y=2./5.*pow((1+cos(Pi*t/(double)numdat/h)),3.);
		break;
	case 4: 
		y=8./35.*pow((1+cos(Pi*t/(double)numdat/h)),4.);
		break;
	case 5: 
		y=8./63.*pow((1+cos(Pi*t/(double)numdat/h)),5.);
		break;
	case 6: 
		y=16./231.*pow((1+cos(Pi*t/(double)numdat/h)),6.);
		break;
	}
	return y;
	
}




void WindData(double *xreal,double *yimag,
			  double h, int numdat,int PP)
{
	int i;
	for(i=0;i<numdat;i++)
	{
		xreal[i]=xreal[i]*windfunc(i*h,h,numdat,PP);
		yimag[i]=yimag[i]*windfunc(i*h,h,numdat,PP);
	}
	
}



void FFT(double *xreal, double *yimag, 
		 double *freq, double *ampl,
		 double h, int numdat, int flag,int PP)
{
	int maxpower,arg,cntr,pnt0,pnt1,i;
	int j,a,b,k;
	double sign,prodreal,prodimag,harm;
	double *cosary;
	double *sinary;
	double *tempx;
	double *tempy;
	
	cosary=(double *) calloc(numdat,sizeof(double));
	if(cosary==NULL) exit(1);
	sinary=(double *) calloc(numdat,sizeof(double));
	if(sinary==NULL) exit(1);
	tempx=(double *) calloc(numdat,sizeof(double));
	if(tempx==NULL) exit(1);
	tempy=(double *) calloc(numdat,sizeof(double));
	if(tempy==NULL) exit(1);
	
	FFT_ArraySetValue(xreal,tempx,numdat);
	FFT_ArraySetValue(yimag,tempy,numdat);
	
	WindData(xreal,yimag,h,numdat,PP);
	
	j=0;
	if(flag!=0)
	{
		sign=1.0;
		for(i=0;i<=numdat-1;i++)
		{
			xreal[i]=xreal[i]/numdat;
			yimag[i]=yimag[i]/numdat;
		}
	}
	else
	{
		sign=-1.0;
	}
	
	for(i=0;i<=numdat-2;i++)
	{
		if(i<j)
		{
			FFT_swap(&xreal[i],&xreal[j]);
			FFT_swap(&yimag[i],&yimag[j]);
		}
		
		k=numdat/2;
		while(k<=j)
		{
			j=j-k;
			k=k/2;
			
		}
		j=j+k;
	}
	maxpower=0;
	i=numdat;
	while(i!=1)
	{
		maxpower=maxpower+1;
		i=i/2;
	}
	harm=2*Pi/numdat;
	for(i=0;i<=numdat-1;i++)
	{
		sinary[i]=sign*sin(harm*i);
		cosary[i]=cos(harm*i);
	}
	a=2;
	b=1;
	for(cntr=1;cntr<=maxpower;cntr++)
	{
		pnt0=numdat/a;
		pnt1=0;
		for(k=0;k<=b-1;k++)
		{
			i=k;
			while(i<numdat)
			{
				arg=i+b;
				if(k==0)
				{
					prodreal=xreal[arg];
					prodimag=yimag[arg];
				}
				else
				{
					prodreal=xreal[arg]*cosary[pnt1]-yimag[arg]*sinary[pnt1];
					prodimag=xreal[arg]*sinary[pnt1]+yimag[arg]*cosary[pnt1];
				}
				xreal[arg]=xreal[i]-prodreal;
				yimag[arg]=yimag[i]-prodimag;
				xreal[i]=xreal[i]+prodreal;
				yimag[i]=yimag[i]+prodimag;
				i=i+a;
			}
			pnt1=pnt1+pnt0;
		}
		a=2*a;
		b=b*2;
	} 
	for(i=0;i<numdat/2;i++)
	{ 
		FFT_swap(&xreal[i],&xreal[i+numdat/2]);
		FFT_swap(&yimag[i],&yimag[i+numdat/2]);
	}
	
	FFT_ArraySetValue(xreal,freq,numdat);
	FFT_ArraySetValue(yimag,ampl,numdat);
	FFT_ArraySetValue(tempx,xreal,numdat);
	FFT_ArraySetValue(tempy,yimag,numdat);
	
	
	free(cosary);
	free(sinary);
	free(tempx);
	free(tempy);
}

void FFTClac(double *xreal, double *yimag,
			 double *freq, double *ampl,
			 double h, int numdat,int PP)
{
	FFT(xreal,yimag,freq,ampl,h,numdat,0,PP);
}

void FFTClacFreqAmplPhase(double *xreal, double *yimag, 
						  double *freq, double *ampl,
						  double h,int numdat,int PP)
{
	int i;
	
	FFTClac(xreal,yimag,freq,ampl,h,numdat,PP);
	for(i=0;i<numdat;i++)
	{
		ampl[i]=sqrt(freq[i]*freq[i]+ampl[i]*ampl[i])/numdat;
		freq[i]=(double)(i-numdat/2.)*2.*Pi/(double)numdat/h; 
	}
	
}