fft1d.cpp

This is a simple one dimensional Fast Fourier Transform (and Inverse FFT) test program. FFT(and Inve

#include "FFT1D.H"
#include <math.h>
#define PI 3.141592

FFT1D::FFT1D()
{
	in_1D_re = NULL;
	in_1D_im = NULL;
	out_1D_re = NULL;
	out_1D_im = NULL;
	MAX = 0;
}

FFT1D::~FFT1D()
{
	delete[] in_1D_re;
	in_1D_re = NULL;
	delete[] in_1D_im;
	in_1D_im = NULL;
	delete[] out_1D_re;
	out_1D_re = NULL;
	delete[] out_1D_im;
	out_1D_im = NULL;
}

void FFT1D::set_in(double* in_re, double* in_im, int MAX){
	in_1D_re = new double[MAX];
	in_1D_im = new double[MAX];
	memcpy(in_1D_re, in_re, MAX*sizeof(double));
	memcpy(in_1D_im, in_im, MAX*sizeof(double));
	this->MAX = MAX;
}

void FFT1D::FFT(BOOL inverse){
	out_1D_re = new double[MAX];
	out_1D_im = new double[MAX];
	
	//decimation in time
	int p=MAX, i, q, r;
	int* n;
	int* n_temp;
	n = new int[MAX];
	n_temp = new int[MAX];

	// reordering input index
	for(i=0; i<=MAX-1; i++)
		n_temp[i]=i;
	while(p>2){
		p=p/2; q=0;
		while(q<MAX){
			for(i=q;i<=q+p-1;i++){
				n[i]=n_temp[q+2*(i-q)];
				n[p+i]=n_temp[q+2*(i-q)+1];
			}
			q=q+2*p;
		}
		for(i=0; i<=MAX-1; i++)
			n_temp[i]=n[i];
	}

	for(i=0; i<=MAX-2; i=i+2){
		out_1D_re[i]=in_1D_re[n[i]]+in_1D_re[n[i+1]];
		out_1D_im[i]=in_1D_im[n[i]]+in_1D_im[n[i+1]];
		out_1D_re[i+1]=in_1D_re[n[i]]-in_1D_re[n[i+1]];
		out_1D_im[i+1]=in_1D_im[n[i]]-in_1D_im[n[i+1]];
	}
	
	delete[] n;
	n = NULL;
	delete[] n_temp;
	n_temp = NULL;

	// 2point-DFT + butterfly + twiddle factor
	double* out_1D_re_temp;
	double* out_1D_im_temp;
	out_1D_re_temp = new double[MAX];
	out_1D_im_temp = new double[MAX];

	p=2;
	while(p<MAX){
		p=p*2; 
		q=p/2;
		r=0;
		for(i=0; i<=MAX-1; i++){
			out_1D_re_temp[i]=out_1D_re[i];
			out_1D_im_temp[i]=out_1D_im[i];
		}
		while(r<MAX){
			for(i=r; i<q+r; i++){
				if(inverse==FALSE){
					out_1D_re[i]=out_1D_re_temp[i]+(out_1D_re_temp[q+i]*cos(2*PI/p*(i-r))+out_1D_im_temp[q+i]*sin(2*PI/p*(i-r)));
					out_1D_im[i]=out_1D_im_temp[i]+(out_1D_im_temp[q+i]*cos(2*PI/p*(i-r))-out_1D_re_temp[q+i]*sin(2*PI/p*(i-r)));
					out_1D_re[q+i]=out_1D_re_temp[i]-(out_1D_re_temp[q+i]*cos(2*PI/p*(i-r))+out_1D_im_temp[q+i]*sin(2*PI/p*(i-r)));
					out_1D_im[q+i]=out_1D_im_temp[i]-(out_1D_im_temp[q+i]*cos(2*PI/p*(i-r))-out_1D_re_temp[q+i]*sin(2*PI/p*(i-r)));
				}
				else if(inverse==TRUE){
					out_1D_re[i]=out_1D_re_temp[i]+(out_1D_re_temp[q+i]*cos(2*PI/p*(i-r))-out_1D_im_temp[q+i]*sin(2*PI/p*(i-r)));
					out_1D_im[i]=out_1D_im_temp[i]+(out_1D_im_temp[q+i]*cos(2*PI/p*(i-r))+out_1D_re_temp[q+i]*sin(2*PI/p*(i-r)));
					out_1D_re[q+i]=out_1D_re_temp[i]-(out_1D_re_temp[q+i]*cos(2*PI/p*(i-r))-out_1D_im_temp[q+i]*sin(2*PI/p*(i-r)));
					out_1D_im[q+i]=out_1D_im_temp[i]-(out_1D_im_temp[q+i]*cos(2*PI/p*(i-r))+out_1D_re_temp[q+i]*sin(2*PI/p*(i-r)));
				}
			}
			r=r+p;
		}
	}
	delete[] out_1D_re_temp;
	out_1D_re_temp = NULL;
	delete[] out_1D_im_temp;
	out_1D_im_temp = NULL;
}

double* FFT1D::get_out_re(){
	return out_1D_re;
}

double* FFT1D::get_out_im(){
	return out_1D_im;
}

void FFT1D::close(){
	delete[] in_1D_re;
	in_1D_re = NULL;
	delete[] in_1D_im;
	in_1D_im = NULL;
	delete[] out_1D_re;
	out_1D_re = NULL;
	delete[] out_1D_im;
	out_1D_im = NULL;
}