gdft-ii.cpp

贡献几个富丽叶变换的很有用的小程序

#include "math.h"
#include "malloc.h"
#include "conio.h"
#include "stdio.h"
#include "string.h"
#include "stdlib.h"
#include "sys/timeb.h"

#define pi (double)3.14159265359
#define m 16

//复数定义
typedef struct 
{
	double re;
	double im;
}COMPLEX;

//复数加运算
COMPLEX Add(COMPLEX c1,COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re+c2.re;
	c.im=c1.im+c2.im;
	return c;
}

//复数乘运算
COMPLEX Mul(COMPLEX c1,COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re*c2.re-c1.im*c2.im;
	c.im=c1.re*c2.im+c2.re*c1.im;
	return c;
}

//复数减运算
COMPLEX Sub(COMPLEX c1,COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re-c2.re;
	c.im=c1.im-c2.im;
	return c;
}

//快速傅里叶变换
void FFT(COMPLEX *TD,COMPLEX *FD,int power)
{
	int count;
	int i,j,k,bfsize,p;
	double angle;
	COMPLEX *W,*X1,*X2,*X;

	count=1<<power;

	W=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	
	for (i=0;i<count/2;i++)
	{
		angle=-i*pi*2/count;
		W[i].re=cos(angle);
		W[i].im=sin(angle);
	}

	memcpy(X1,TD,sizeof(COMPLEX)*count);

	for (k=0;k<power;k++)
	{
		for (j=0;j<1<<k;j++)
		{
			bfsize=1<<(power-k);
			for (i=0;i<bfsize/2;i++)
			{
				p=j*bfsize;
				X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
				X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
			}
		}
		X=X1;
		X1=X2;
		X2=X;
	}

	for (j=0;j<count;j++)
	{
		angle=-j*pi/count;
		W[j].re=cos(angle);
		W[j].im=sin(angle);
		p=0;
		for (i=0;i<power;i++)
		{
			if (j&(1<<i))
				p+=1<<(power-i-1);
		}
		FD[j]=Mul(X1[p],W[j]);
	}
	
	free(W);
	free(X1);
	free(X2);
}

void main()
{
	double ComMs;
	struct _timeb timebuffer,timebuffer2;
	int count;
	COMPLEX *TV,*FV;
	int i;
	count=1<<m;
	TV=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	FV=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	for (i=0;i<count;i++)
	{
		TV[i].re=i*3/5.0;
		TV[i].im=i-5;
	}
	_ftime(&timebuffer );
	FFT(TV,FV,m);
	_ftime(&timebuffer2 );
	ComMs=(timebuffer2.time-timebuffer.time)+float(timebuffer2.millitm-timebuffer.millitm)/1000;
    /*for (i=0;i<count;i++)
	{
		printf("x[%d]=%f+i %f\n",i,FV[i].re,FV[i].im);
	}*/
	printf("运算时间为%f秒。\n",ComMs);
	free(TV);
	free(FV);
	getch();
}