mixfft.cpp

fft混合基算法

#include <fstream>
#include <cmath>

#define r1 3
#define r2 5
#define N  15
#define pi 3.141592653589793
using namespace std;

double xR[N], xI[N];
double x1R[r1][r2],x1I[r1][r2];
int inv=0;
FILE *in, *out, *fopen();

void read_data();
void write_data(char*);
void fft();
void ifft();

void main()
{
	char* out1 = "datafft.txt";
	char* out2 = "dataifft.txt";
	
	read_data();
	fft();
	write_data(out1);
	ifft();
	write_data(out2);
}

void fft()
{
	int k, n0, n1;
	double XR[N],XI[N];
	double WN,wk;
	double tempR, tempI;

	for(n1=0,k=0; n1<r1; n1++)                  //改写成数组形式
		for(n0=0; n0<r2; n0++)
		{
			x1R[n1][n0]=xR[k]; 
			x1I[n1][n0]=xI[k];
			k++;
		}

	for(n0=0; n0<r2; n0++)                  //r2个r1点的DFT
	{
		WN=-2*pi/r1;
		
		if(inv==1)
			WN = -WN;

		for(k=0;k<r1;++k)
		{
		    XR[k]=0;
			XI[k]=0;
			wk=k*WN;

			for(n1=0;n1<r1;n1++)
			{
				XR[k]=XR[k]+x1R[n1][n0]*cos(n1*wk)-x1I[n1][n0]*sin(n1*wk);
				XI[k]=XI[k]+x1R[n1][n0]*sin(n1*wk)+x1I[n1][n0]*cos(n1*wk);
			}
		    if(inv==1)
			{
				XR[k]=XR[k]/r1; 
				XI[k]=XI[k]/r1;
			}
		}

		for(k=0;k<r1;++k)
		{
			x1R[k][n0]=XR[k];
			x1I[k][n0]=XI[k];
		}

	}

	for(n1=0; n1<r1; n1++)                              //乘以旋转因子
	{
		WN = -2*pi/N;
		if(inv==1)
			WN = -WN;
		for(n0=0; n0<r2; n0++)
		{
			tempR = x1R[n1][n0]*cos(WN*n1*n0)-x1I[n1][n0]*sin(WN*n1*n0);
			tempI = x1R[n1][n0]*sin(WN*n1*n0)+x1I[n1][n0]*cos(WN*n1*n0);
			x1R[n1][n0]=tempR;
			x1I[n1][n0]=tempI;
		}
	}

	for(n1=0; n1<r1; n1++)                  //r1个r2点的DFT
	{
		WN=-2*pi/r2;

		if(inv==1)
			WN = -WN;

		for(k=0;k<r2;++k)
		{
		    XR[k]=0;
			XI[k]=0;
			wk=k*WN;

			for(n0=0;n0<r2;n0++)
			{
				XR[k]=XR[k]+x1R[n1][n0]*cos(n0*wk)-x1I[n1][n0]*sin(n0*wk);
				XI[k]=XI[k]+x1R[n1][n0]*sin(n0*wk)+x1I[n1][n0]*cos(n0*wk);
			}
			if(inv==1)
			{
				XR[k]=XR[k]/r2; 
				XI[k]=XI[k]/r2;
			}
		}

		for(k=0;k<r2;++k)
		{
			x1R[n1][k]=XR[k];
			x1I[n1][k]=XI[k];
		}
	}

	for(n1=0; n1<r1; n1++)                     //调整顺序
	{
		for(n0=0; n0<r2; n0++)
		{
			xR[r1*n0+n1] = x1R[n1][n0];
			xI[r1*n0+n1] = x1I[n1][n0];
		}
	}


}

void ifft()
{
	inv=1;
	fft();
}

void read_data()
{
	int n;
	if ((in=fopen("data.txt","rt"))==NULL)
	{
		printf("cannot open the file data.txt\n");
		exit(1);
	}	
	for(n=0;n<N;n++)
	{
		fscanf(in,"%lf",&xR[n]);
		xI[n] = 0;
	}

    fclose(in);
}

void write_data(char* name)
{
	int n;
	if ((out=fopen(name,"w"))==NULL)
	{
		printf("cannot open the file datafft.txt\n");
		exit(1);
	}	
	for(n=0; n<N; n++)
		fprintf(out,"%lf \t %lf \n",xR[n],xI[n]);			
    fclose(out);
}