fft.cpp

此程序事一个关于快速傅里叶变换的vc程序

/**************************************************/
#include "stdio.h"
#include"dos.h"
#include "stdlib.h"
#include"malloc.h"
#include "math.h"
/***********************************************************
 *Filename:dsp.h                                           *
 *Function:This file includes the function of FFT          *
 **********************************************************/
typedef struct{
	          float real;
			  float imag;
}COMPLEX;
extern void fft(COMPLEX*x,int m);
extern void ifft(COMPLEX*x,int m);
void fft(COMPLEX*x,int m)
{
	static COMPLEX*w;
	static int mstore=0;
	static int n=1;
	
    COMPLEX u,temp,tm;
	COMPLEX *xi,*xip,*xj,*wptr;

	int i,j,k,l,le,windex;
	double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
	if(m!=mstore)
	{
		if(mstore!=0)free(w);
		mstore=m;
		if(m==0)return;

		n=1<<m;
		le=n/2;

		w=(COMPLEX*) calloc(le-1,sizeof(COMPLEX));
		if(!w)
		{
			printf("\nUnable to allocade complex w array\n");
			exit(1);
		}

		arg=4.0*atan(1.0)/le;
		wrecur_real=w_real=cos(arg);
		wrecur_imag=w_imag=-sin(arg);
		xj=w;
		for(j=1;j<le;j++)
		{
			xj->real=(float)wrecur_real;
			xj->imag=(float)wrecur_imag;
			xj++;
			wtemp_real=wrecur_real*w_real-wrecur_imag*w_imag;
			wrecur_imag=wrecur_real*w_imag+wrecur_imag*w_real;
			wrecur_real=wtemp_real;
		}
	}

	le=n;
	windex=1;
	for(l=0;l<m;l++)
	{
		le=le/2;
		for(i=0;i<n;i=i+2*le)
		{
			xi=x+i;
			xip=xi+le;
			temp.real=(xi->real+xip->real);
			temp.imag=(xi->imag+xip->imag);
			xip->real=(xi->real-xip->real);
			xip->imag=(xi->imag-xip->imag);
			*xi=temp;
		}
		wptr=w+windex-1;
		for(j=1;j<le;j++)
		{
			u=*wptr;
				for(i=j;i<n;i=i+2*le)
				{
					xi=x+i;
					xip=xi+le;
					temp.real=(xi->real+xip->real);
					temp.imag=(xi->imag+xip->imag);
					tm.real=xi->real-xip->real;
					tm.imag=xi->imag-xip->imag;
					xip->real=(tm.real*u.real-tm.imag*u.imag);
					xip->imag=(tm.real*u.imag+tm.imag*u.real);
					*xi=temp;
				}
				wptr=wptr+windex;
		}
		windex=2*windex;
	}

	for(i=0;i<n;++i)
	{
		j=0;
		for(k=0;k<m;++k);
		j=(j<<1)|(1&(i>>k));
		if(i<j)
		{
			xi=x+i;
			xj=x+j;
			temp=*xj;
			*xj=*xi;
			*xi=temp;
		}
	}
}

void main()
{
	COMPLEX *x,*xtemp;
	int i,nn,mm;
	double y,q;
	system("cls");
	printf("\n***************************************");
	printf("\n*Filename:fftl.c                      *");
	printf("\n*function:fft analysis for signals    *");
	printf("\n*signal:x(n)=exp(-n)                  *");
	printf("\n***************************************");
	printf("\n Input m=");
	scanf("%d",&mm);
	nn=1<<mm;
	x=(COMPLEX*)calloc(nn,sizeof(COMPLEX));
	xtemp=(COMPLEX*)calloc(nn,sizeof(COMPLEX));
	/*store x(n) in the array*/
	x[0].real=1.0;
	x[0].imag=0.0;
	for(i=1;i<nn;i++)
	{
		x[i].real=(float)exp(-i);
		x[i].imag=0.0;
	}
	for(i=0;i<nn;i++) xtemp[i]=x[i];
	/*FFT analysis for x(n)*/
fft(x,mm);
printf("Finished!\n\n");
/*output results in a format*/
printf("  k       XR(K)          XI(K)         |X(K)|          Q(K)\n");
printf("-----------------------------------------------------------\n");
for(i=0;i<nn;i++)
{
	y=x[i].real*x[i].real+x[i].imag*x[i].imag;
	if(x[i].real!=0)
	{
		if(x[i].imag==0.0)
		{
			if(x[i].real>0) q=0.0;
			if(x[i].real<0) q=4.0*(atan(1.0));
		}
		q=atan(x[i].imag/x[i].real);
	}
	printf("%4d     %9f     %9f     %9f     %9f\n",i,x[i].real,x[i].imag,sqrt(y),q);
}
printf("\n\n");
}