fft.cpp

快速傅立叶算法

#include "stdio.h"
#include "iostream.h"
#include "complex"
using namespace std;

#define N1 64
#define N2 128
#define N3 256

FILE *fp;

bool reverse(complex<double> *a,unsigned n);
int dit_fft(complex<double> *a,int l)
{
	const double pai=3.141592653589793;
	complex<double> u,w,t;
	unsigned n=1,le,lei,ip;
	unsigned i,j,m;
	double temp;

	n<<=l;
	reverse(a,n);
	le=1;
	for(m=1;m<=l;m++)
	{
		lei=le;
		le<<=1;
		u=complex<double> (1,0);;
		temp=pai/lei;
		w=complex<double> (cos(temp),-sin(temp));
		for(j=0;j<lei;j++)
		{
			for(i=j;i<n;i+=le)
			{
				ip=i+lei;
				t=a[ip]*u;
				a[ip]=a[i]-t;
				a[i]+=t;
			}
			u*=w;
		}
	}
	return 0;
}

int dif_fft(complex<double> *a,int l)
{
	const double pai=3.141592653589793;
	complex<double> u,w,t;
	unsigned n=1,le,lei,ip;
	unsigned i,j,m;
	double temp;

	n<<=l;
	lei=n;
	for(m=l;m>0;m--)
	{
		le=lei;		//le:蝶形之间的距离
		lei>>=1;	//lei：蝶形的距离的一半，碟形两个输入节点之间的距离
		u=complex<double> (1,0);
		temp=pai/lei;
		w=complex<double> (cos(temp),-sin(temp));
		for(j=0;j<lei;j++)
		{
			for(i=j;i<n;i+=le)
			{
				ip=i+lei;
				t=a[i]+a[ip];
				a[ip]*=-1;
				a[ip]+=a[i];
				a[ip]*=u;
				a[i]=t;
			}
			u*=w;
		}
	}
	reverse(a,n);
	return 0;
}

bool reverse(complex<double> *a,unsigned n)
{
	unsigned nv2,nm1,k;
	unsigned i,j;
	complex<double> t;
	nv2=n>>1;
	nm1=n-1;
	j=0;
	for(i=0;i<nm1;i++)
	{
		if(i<j)
		{
			t=a[j];
			a[j]=a[i];
			a[i]=t;
		}
		k=nv2;
		while(k<=j)
		{
			j-=k;
			k>>=1;
		}
		j+=k;
	}
	return true;
}
bool fill(complex<double> *a,unsigned n)
{
	for(int i=0;i<n;i++)
	{
		a[i]=sin(1.6964600329384883487698274269709*i)+
			sin(1.8456856839840035275968029876767*i)+
			sin(1.9713493901275952571353087230079*i);
	}
	return true;
}

bool printResult(complex<double> *a,unsigned n)
{
	for(int i=0;i<n;i++)
	{
		if(i!=0) fprintf(fp,",");
		if(a[i].imag()>0)
		{
			fprintf(fp,"%lf+%lfi",a[i].real(),a[i].imag());
		}
		if(a[i].imag()<0)
		{
			fprintf(fp,"%lf%lfi",a[i].real(),a[i].imag());
		}		
		if(a[i].imag()==0)
		{
			fprintf(fp,"%lf",a[i].real());
		}
		
	}
	return true;
}

int main()
{
	complex<double> x1[N1],x2[N2],x3[N3];
	if((fp=fopen("dit_fft_64.csv","w+"))==NULL) printf("File open error!!!\n");
	fill(x1,N1);
	dit_fft(x1,log10(N1)/log10(2));
	printResult(x1,N1);

	if((fp=fopen("dif_fft_64.csv","w+"))==NULL) printf("File open error!!!\n");
	fill(x1,N1);
	dif_fft(x1,log10(N1)/log10(2));
	printResult(x1,N1);


	if((fp=fopen("dit_fft_128.csv","w+"))==NULL) printf("File open error!!!\n");
	fill(x2,N2);
	dit_fft(x2,log10(N2)/log10(2));
	printResult(x2,N2);

	if((fp=fopen("dif_fft_128.csv","w+"))==NULL) printf("File open error!!!\n");
	fill(x2,N2);
	dif_fft(x2,log10(N2)/log10(2));
	printResult(x2,N2);


	if((fp=fopen("dit_fft_256.csv","w+"))==NULL) printf("File open error!!!\n");
	fill(x3,N3);
	dit_fft(x3,log10(N3)/log10(2));
	printResult(x3,N3);

	if((fp=fopen("dif_fft_256.csv","w+"))==NULL) printf("File open error!!!\n");
	fill(x3,N3);
	dif_fft(x3,log10(N3)/log10(2));
	printResult(x3,N3);
	return 0;
}