simple.c

优化过的128点fft和ifft

/*时间抽选基2 FFT及IFFT算法C语言实现*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/timeb.h>

#include "complex.h"

#define N 128//输入序列最大长度
fft();    /*快速傅里叶变换*/
void ifft();
void initW();  /*初始化变换核*/
void change(); /*变址*/
change2();

void output();                            /*输出结果*/
complex fftout[128],ifftout[128];

int size_x=0;     /*输入序列的大小，在本程序中仅限2的次幂*/
double PI;        /*圆周率*/

unsigned short sequence[128]={0,64,32,96,16,80,48,112,8,72,40,104,24,88,56,120,4,68,36,100,
	20,84,52,116,12,76,44,108,28,92,60,124,2,66,34,98,18,82,50,114,10,74,42,106,26,90,58,
	122,6,70,38,102,22,86,54,118,14,78,46,110,30,94,62,126,1,65,33,97,17,81,49,113,9,73,41,
	105,25,89,57,121,5,69,37,101,21,85,53,117,13,77,45,109,29,93,61,125,3,67,35,99,19,83,
	51,115,11,75,43,107,27,91,59,123,7,71,39,103,23,87,55,119,15,79,47,111,31,95,63,127};
unsigned short * sequence1;
int main()
{	FILE * f2;
	int i,method,k;
	
	double a[128]={0.0998,0.1987,0.2955,0.3894,0.4794,0.5646,0.6442,0.7174,0.7833,0.8415,0.8912,0.932,
		0.9636,0.9854,0.9975,0.9996,0.9917,0.9738,0.9463,0.9093,0.8632,0.8085,0.7457,0.6755,
		0.5985,0.5155,0.4274,0.335,0.2392,0.1411,0.0416,-0.0584,-0.1577,-0.2555,-0.3508,-0.4425,
		-0.5298,-0.6119,-0.6878,-0.7568,-0.8183,-0.8716,-0.9162,-0.9516,-0.9775,-0.9937,-0.9999,
		-0.9962,-0.9825,-0.9589,-0.9258,-0.8835,-0.8323,-0.7728,-0.7055,-0.6313,-0.5507,-0.4646,
		-0.3739,-0.2794,-0.1822,-0.0831,0.0168,0.1165,0.2151,0.3115,0.4048,0.4941,0.5784,0.657,
		0.729,0.7937,0.8504,0.8987,0.938,0.9679,0.9882,0.9985,0.9989,0.9894,0.9699,0.9407,0.9022,
		0.8546,0.7985,0.7344,0.663,0.5849,0.501,0.4121,0.3191,0.2229,0.1245,0.0248,
		-0.0752, -0.1743,-0.2718,-0.3665,-0.4575,-0.544	,-0.6251,-0.6999,-0.7677,
		-0.8278,-0.8797,-0.9228,-0.9566,-0.9809,-0.9954,-1,-0.9946,
		-0.9792,-0.954,-0.9193,-0.8755,-0.8228,-0.762,-0.6935,-0.6181,-0.5366,
		-0.4496,-0.3582,-0.2632,-0.1656,-0.0663,0.0336,	0.1332,	0.2315};
	
	double b[256]={0.3868,0,0.5468,-1.4729,13.7483,-61.7872,-0.5433,2.8661,-0.2901,1.5619,-0.214,
		1.1075,-0.1798,0.8681,-0.1607,0.7175,-0.1493,0.6125,-0.1415,0.5346,-0.1364,0.4745,-0.1316,
		0.4262,-0.1291,0.3865,-0.127,0.3533,-0.1252,0.3256,-0.1235,0.3009,-0.1225,0.2794,-0.1213,
		0.2603,-0.1204,0.244,-0.1202,0.229,-0.1195,0.2153,-0.119,0.2026,-0.1187,0.1918,-0.1184,
		0.1812,-0.1177,0.1713,-0.1175,0.1631,-0.117,0.1545,-0.1175,0.147,-0.1173,0.1397,-0.1168,
		0.1324,-0.117,0.1261,-0.1161,0.1202,-0.1163,0.1147,-0.1161,0.1089,-0.1159,0.1034,-0.1159,
		0.0987,-0.1157,0.0936,-0.1162,0.0892,-0.116,0.085,-0.1152,0.0803,-0.1154,0.0765,-0.1159,
		0.0724,-0.1152,0.0685,-0.1154,0.0649,-0.1155,0.0606,-0.116,0.0572,-0.1158,0.0536,-0.1151,
		0.0506,-0.1155,0.0472,-0.1156,0.0441,-0.115,0.0408,-0.1155,0.0379,-0.1149,0.035,-0.1155,
		0.0318,-0.1151,0.0288,-0.1152,0.0257,-0.1152,0.0225,-0.1152,0.0193,-0.1152,0.0166,-0.1151,
		0.0141,-0.1154,0.0111,-0.115,0.0084,-0.1151,0.0057,-0.1151,0.003,-0.115,0,-0.1151,-0.003,
		-0.1151,-0.0057,-0.115,-0.0084,-0.1154,-0.0111,-0.1151,-0.0141,-0.1152,-0.0166,-0.1152,
		-0.0193,-0.1152,-0.0225,-0.1152,-0.0257,-0.1151,-0.0288,-0.1155,-0.0318,-0.1149,-0.035,
		-0.1155,-0.0379,-0.115,-0.0408,-0.1156,-0.0441,-0.1155,-0.0472,-0.1151,-0.0506,-0.1158,
		-0.0536,-0.116,-0.0572,-0.1155,-0.0606,-0.1154,-0.0649,-0.1152,-0.0685,-0.1159,-0.0724,
		-0.1154,-0.0765,-0.1152,-0.0803,-0.116,-0.085,-0.1162,-0.0892,-0.1157,-0.0936,-0.1159,
		-0.0987,-0.1159,-0.1034,-0.1161,-0.1089,-0.1163,-0.1147,-0.1161,-0.1202,-0.117,-0.1261,
		-0.1168,-0.1324,-0.1173,-0.1397,-0.1175,-0.147,-0.117,-0.1545,-0.1175,-0.1631,-0.1177,
		-0.1713,-0.1184,-0.1812,-0.1187,-0.1918,-0.119,-0.2026,-0.1195,-0.2153,-0.1202,-0.229,
		-0.1204,-0.244,-0.1213,-0.2603,-0.1225,-0.2794,-0.1235,-0.3009,-0.1252,-0.3256,-0.127,
		-0.3533,-0.1291,-0.3865,-0.1316,-0.4262,-0.1364,-0.4745,-0.1415,-0.5346,-0.1493,-0.6125,
		-0.1607,-0.7175,-0.1798,-0.8681,-0.214,-1.1075,-0.2901,-1.5619,-0.5433,-2.8661,13.7483,
		61.7872,0.5468,1.4729};
    
	
	char tmpbuf[128];
	struct _timeb tstruct;
	unsigned short counter;
	
	sequence1=(unsigned short *)malloc(sizeof(unsigned short) * size_x);
	for(i=0;i<128;i++)
		sequence1[i]=sequence[i];
	//free(sequence1);
	size_x=N;	
	PI=atan(1)*4;

	initW();
	printf("Use FFT(0) or IFFT(1)?\n");
	scanf("%d",&method);
	_tzset();
    /*记录起始时间*/
    _strtime(tmpbuf);
    _ftime( &tstruct );
	printf( "start time:\t\t\t\t%s\n", tmpbuf );
	printf( "Plus milliseconds:\t\t\t%u\n", tstruct.millitm );
   
	for(k=0;k<10000;k++)
	{
		/*获取输入序列*/ 
		counter=0;
		/*for(i=0;i<size_x;i++)
		{//获取正变换序列,实序列的奇数和偶数点分别作为新序列的实部和虚部
			x[i].real=a[i];
			x[i].img=0;
		}*/			
		
		for(i=0;i<256;i=i+2)
		{//获取反变换序列
			x[counter].real=b[i];
			x[counter].img=-b[i+1];//虚部取反求共轭，为逆变换调用fft子程序服务
			counter++;
		}	
		
	
		if(method==0)//计算傅里叶变换
			fft(x,fftout);

		else if(method==1)//计算傅里叶逆变换
		{
			fft(x,ifftout);
			for(i=0;i<N;i++)
			{
				ifftout[i].real=ifftout[i].real/N;
				ifftout[i].img=-ifftout[i].img/N;
			}
		}
	}
	
	/*记录结束时间*/
	_strtime(tmpbuf);
	_ftime( &tstruct );	
	printf( "end time:\t\t\t\t%s\n", tmpbuf );
	printf( "Plus milliseconds:\t\t\t%u\n", tstruct.millitm );
	free(W);
	free(sequence1);
	if(method==0)
	{
		output(fftout);
		f2=fopen("fftoutput.txt","w");
		for(i=0;i<N;i++)
			fprintf(f2,"%lf+%lfj ",fftout[i].real,fftout[i].img);
	}
	else if(method==1)
	{
		output(ifftout);
		f2=fopen("ifftoutput.txt","w");
		for(i=0;i<N;i++)
			fprintf(f2,"%lf+%lfj ",ifftout[i].real,ifftout[i].img);
	}
	fclose(f2);
	return 0;
	}

/*快速傅里叶变换*/
fft(complex x1[N],complex x2[N])
{   
	unsigned short i=0,j=0,k=0,l=0;
	complex up,down,product;
	change2(x1,x2);
	for(i=0;i< log(size_x)/log(2) ;i++)
	{  //一级蝶形运算,每级有size_x/2l组
		l=1<<i;//时域抽取间隔，第1级间隔1点抽取，第2级间隔2点抽取...
		for(j=0;j<size_x;j+= 2*l)
		{            //一组蝶形运算，每组l个碟形运算
			for(k=0;k<l;k++)
			{       //一个蝶形运算
				//BC
				product.real=x2[j+k+l].real*W[size_x*k/2/l].real-x2[j+k+l].img*W[size_x*k/2/l].img;
				product.img=x2[j+k+l].img*W[size_x*k/2/l].real+x2[j+k+l].real*W[size_x*k/2/l].img;
				//A+BC
				up.real=x2[j+k].real+product.real;
				up.img=x2[j+k].img+product.img;
				//A-BC
				down.real=x2[j+k].real-product.real;
				down.img=x2[j+k].img-product.img;
				x2[j+k]=up;
				x2[j+k+l]=down;
			}
		}
	}
	
}



/*初始化变换核*/
void initW()
{
	int i;
	W=(complex *)malloc(sizeof(complex) * size_x);
	for(i=0;i<size_x;i++)
	{
		W[i].real=cos(2*PI/size_x*i);
		W[i].img=-1*sin(2*PI/size_x*i);
	}
}



/*二进制倒序法实现时域抽取,yonghua20090310*/
change2(complex in_change[N],complex out_change[N])
{
	int i;
	for(i=0;i<128;i++)
	{	out_change[i]=in_change[sequence[i]];
		//printf("%d,%lf ",sequence[i],out_change[i].real);
	}
}

/*输出傅里叶变换的结果*/
void output(complex x[N])
{
	int i;
	printf("The result are as follows\n");
	for(i=0;i<size_x;i++)
	{
		printf("%.4f",x[i].real);
		if(x[i].img>=0.0001)printf("+%.4fj\n",x[i].img);
		else if(fabs(x[i].img)<0.0001)
			printf("\n");
		else 
			printf("%.4fj\n",x[i].img);
	}
}

/*二进制倒序法实现时域抽取,采用循环移位法获取下标*/
/*change2(complex in_change[N],complex out_change[N])
{
	int i;
	unsigned short j,n,sequence;
	in_change[128];
	n=128;
	j=0;	
	
	//f1=fopen("sequence.txt","w");
	for(i=0;i<128;i++)
	{   j=i;
		sequence=((j>>6|j<<2)&17)|((j>>4|j<<4)&34)|((j>>2|j<<6)&68)|(j&8);
		//printf("%d ",sequence);
		//fprintf(f1,"%d ",a);
		out_change[i]=in_change[sequence];
	}
	//fclose(f1);
}*/