gdft-i.cpp

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

#include "math.h"
#include "malloc.h"
#include "conio.h"
#include "stdio.h"
#include "string.h"

#define pi (double)3.14159265359
#define m 6

//复数定义
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++)
	{
		p=0;
		for (i=0;i<power;i++)
		{
			if (j&(1<<i))
				p+=1<<(power-i-1);
		}
		FD[j]=X1[p];
	}
	
	free(W);
	free(X1);
	free(X2);
}

//逆变换
void IFFT(COMPLEX *FD,COMPLEX *TD,int power)
{
	int i,count;
	COMPLEX *x;

	count=1<<power;
	x=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
	memcpy(x,FD,sizeof(COMPLEX)*count);

	for (i=0;i<count;i++)
	{
		x[i].im=-x[i].im;
	}
	FFT(x,TD,power);

	for(i=0;i<count;i++)
	{
		TD[i].re/=count;
		TD[i].im=-TD[i].im/count;
	}
	free(x);
}

void main()
{
	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=0;
	}
	FFT(TV,FV,m);
	IFFT(FV,TV,m);
	for (i=0;i<count;i++)
	{
		printf("x[%d]=%f+i %f\n",i,TV[i].re,TV[i].im);
	}
	free(TV);
	free(FV);
	getch();
}