#if !defined(FFT_HEADER_HH_INCLUDED)
#define FFT_HEADER_HH_INCLUDED
//
// FFT header file for epsiron image processor fft-header.hh
//
// Ito Lab. at Naruto University of Education, Japan
//
// Version history:
// Ver 0.5 : 2007-10-28 Original version
// Ver 0.8 : 2007-11-11 Revise file input / output
//
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#include <complex>

using namespace std;

FILE * fft_stdin  = stdin;
FILE * fft_stdout = stdout;
FILE * fft_stderr = stderr;

typedef complex<double> creal8;

const int FFT_MAX_FILENAME_LEN = 2048;		// Maximum length of file name

// "Complex Text" format
const char * COMPLEX_TEXT_ID    = "COMPLEX-TEXT";
enum em_complex_type {
	COMPLEX_SPACE = 1,
	COMPLEX_FREQ
};

const char * COMPLEX_TEXT_SPACE = "SPACE";
const char * COMPLEX_TEXT_FREQ  = "FREQ";

//
// Fourier transform class
//
class fourier
{
public:
	creal8 ** p;		// Complex 2D array p[Line_No][Pixel_No]
	creal8 * z;			// Complex 1D array
						// z = p[0];

	creal8 * w_line;	// Use in FFT
	creal8 * w_pixel;	// Use in FFT

	static const int RET_ERROR =  0;		// Return code with error
	static const int RET_OK    = -1;		// Return code with no error

	static const int FORWARD_FT  = -1;
	static const int BACKWARD_FT = +1;

private:
	em_complex_type	complex_type;
	int		nr_lines;			// Number of lines
	int		nr_pixels;			// Number of pixels

public:
	// Constructors
	fourier();
	fourier(					// 1D
		int nr_pixels_arg);
	fourier(					// 2D
		int nr_lines_arg,
		int nr_pixels_arg);

	// Destructor
	~fourier();

	fourier(const fourier &src);			// Copy constructor
	fourier& operator=(const fourier &src);	// Substitution operator =

	void clear(void);			// Clear all image data

	int get_nr_lines(void);		// Get number of lines
	int get_nr_pixels(void);	// Get number of pixels

	void set_complex_type(
		em_complex_type complex_type_arg);
	em_complex_type get_complex_type(void);

	// Allocate or resize complex array
	void cmake(
		int nr_lines_arg,
		int nr_pixels_arg);

	int Forward_DFT(void);
	int Backward_DFT(void);

	int Forward_FFT(void);
	int Backward_FFT(void);

// "Complex Text" file format is as follows.
//
//  1:COMPLEX-TEXT {SPACE, FREQ} {Number of lines} {Number of pixels}
//  2:               {Pixel number} Re Im} ....
//  3: {Line number} {Real part} {Imaginary part} ....
//
//  Where Real and imaginary parts format is "%+14.6lle".
	void cpx_load(void);		// Load complex from "Complex Text" file
	void cpx_load(char * fn);

	void cpx_save();			// Save complex to "Complex Text" file
	void cpx_save(char * fn);

	void MakeImage(
		int	output_data,
		int	shuffling,
		int	log_scale,
		unsigned char ** pw);

private:
	int alloc(
		int	nr_lines_arg,		// Number of lines in this image
		int nr_pixels_arg);		// Number of pixels in this image
	void unalloc();
	void copy_proc(const fourier &src);	// Use in copy constructor and '='

	int DFT(int dir);

	int FFT(int dir);

	void fft_1d(
		int		dir,
		creal8	* data,
		int		fft_len,
		creal8	* w);

	double GetOutputData(
		int output_data,
		creal8 x);
};

// Constructor of image class
fourier::fourier()
{
	p = NULL;
	z = NULL;
	nr_lines  = 0;
	nr_pixels = 0;
	w_line  = NULL;
	w_pixel = NULL;
	complex_type = COMPLEX_SPACE;
}

fourier::fourier(
	int nr_pixels_arg)		// Number of pixels in this image
{
	if (alloc(1, nr_pixels_arg) == RET_OK) {
		clear();
	}
	complex_type = COMPLEX_SPACE;
}

fourier::fourier(
	int nr_lines_arg,		// Number of lines in this image
	int nr_pixels_arg)		// Number of pixels in this image
{
	if (alloc(nr_lines_arg, nr_pixels_arg) == RET_OK) {
		clear();
	}
	complex_type = COMPLEX_SPACE;
}

// Destructor of image class
fourier::~fourier()
{
	unalloc();
}

// Copy constructor
fourier::fourier(
	const fourier	&src)
{
	copy_proc(src);
}

// Substitution operator =
fourier& fourier::operator=(const fourier &src)
{
	unalloc();

	copy_proc(src);
	return (*this);
}

// Copy members
void fourier::copy_proc(
	const fourier	&src)
{
	if (src.p == NULL) return;

	if (alloc(src.nr_lines, src.nr_pixels) == RET_OK) {
		for (int iy = 0; iy < nr_lines; iy ++) {
			for (int ix = 0; ix < nr_pixels; ix ++) {
				p[iy][ix] = src.p[iy][ix];
			}
		}

		if (w_line != NULL) {
			for (int iy = 0; iy < nr_lines; iy ++) {
				w_line[iy] = src.w_line[iy];
			}
		}

		if (w_pixel != NULL) {
			for (int ix = 0; ix < nr_pixels; ix ++) {
				w_pixel[ix] = src.w_pixel[ix];
			}
		}
	}
}

//
// Memory allocation for fourier array
//
// return code: RET_OK, RET_ERROR
//
int fourier::alloc(
	int	nr_lines_arg,		// Number of lines in this image
	int nr_pixels_arg)		// Number of pixels in this image
{
	int error_flag = RET_OK;
	if (nr_lines_arg  == 0 || nr_lines_arg  < 0 ||
	    nr_pixels_arg == 0 || nr_pixels_arg < 0) {
		error_flag = RET_ERROR;
	} else {
		nr_lines  = nr_lines_arg;
		nr_pixels = nr_pixels_arg;
		p = (creal8 **)malloc(sizeof(creal8 *) * nr_lines);
		if (p == NULL) {
			error_flag = RET_ERROR;
		} else {
			for (int iy = 0; iy < nr_lines; iy ++) {
				p[iy] = (creal8 *)malloc(sizeof(creal8) * nr_pixels);
				if (p[iy] == NULL) {
					error_flag = RET_ERROR;
					break;
				}
			}
		}
	}

	if (nr_pixels >= 2) {
		w_pixel = (creal8 *)malloc(sizeof(creal8) * nr_pixels);
		if (w_pixel != NULL) {
			for (int ix = 0; ix < nr_pixels; ix ++) {
				w_pixel[ix] = exp(creal8(0.0, -1.0) * (2.0 * PI / nr_pixels) * (double)ix);
			}
		} else {
			error_flag = RET_ERROR;
		}
	} else {
		w_pixel = NULL;
	}

	if (nr_lines >= 2) {
		w_line = (creal8 *)malloc(sizeof(creal8) * nr_lines);
		if (w_line != NULL) {
			for (int iy = 0; iy < nr_lines; iy ++) {
				w_line[iy] = exp(creal8(0.0, -1.0) * (2.0 * PI / nr_lines) * (double)iy);
			}
		} else {
			error_flag = RET_ERROR;
		}
	} else {
		w_line = NULL;
	}

	if (error_flag == RET_ERROR) {
		printf("Memory allocation error in image construction\n");
		p = NULL;
		z = NULL;
		nr_lines  = 0;
		nr_pixels = 0;
	} else {
		z = p[0];
	}
	return error_flag;
}

//
// Memory unallocation for image array
//
void fourier::unalloc()
{
	if (p != NULL) {
		for (int iy = 0; iy < nr_lines; iy ++) {
			free(p[iy]);
		}
		free(p);
	}
	p = NULL;
	z = NULL;

	if (w_line  != NULL) free(w_line);
	if (w_pixel != NULL) free(w_pixel);

	nr_lines  = 0;
	nr_pixels = 0;
}

//
// Clear all image data
//
void fourier::clear(void)
{
	if (p != NULL) {
		for (int iy = 0; iy < nr_lines; iy ++) {
			for (int ix = 0; ix < nr_pixels; ix ++) {
				p[iy][ix] = creal8(0.0, 0.0);
			}
		}
	}
}

int fourier::get_nr_lines(void)
{
	return nr_lines;
}

int fourier::get_nr_pixels(void)
{
	return nr_pixels;
}

void fourier::set_complex_type(
		em_complex_type complex_type_arg)
{
	complex_type = complex_type_arg;
}

em_complex_type fourier::get_complex_type(void)
{
	return complex_type;
}

//
// Make memory block for this complex array
//
void fourier:: cmake(
	int nr_lines_arg,
	int nr_pixels_arg)
{
	if (alloc(nr_lines_arg, nr_pixels_arg) == RET_OK) {
		clear();
	}
}

//
// Discrete Fourier Transform
//
int fourier::Forward_DFT(
	void)
{
	return DFT(FORWARD_FT);
}

int fourier::Backward_DFT(
	void)
{
	return DFT(BACKWARD_FT);
}

int fourier::DFT(
	int	dir)		// dir = -1  ==> Forward  DFT
					// dir = +1  ==> Backward DFT
{
	if (nr_pixels >= 2) {
		creal8 * u = (creal8 *)malloc(sizeof(creal8) * nr_pixels);
		if (u == NULL) return RET_ERROR;

		for (int iy = 0; iy < nr_lines; iy++) {
			for (int ix = 0; ix < nr_pixels; ix ++) {
				u[ix] = p[iy][ix];
			}

			for (int ix = 0; ix < nr_pixels; ix ++) {
				creal8 s(0.0, 0.0);
				creal8 w = creal8(0.0, (double)dir) * (2.0 * PI / nr_pixels) * (double)ix;
				for (int j = 0; j < nr_pixels; j ++) {
					s += u[j] * exp(w * (double)j);
				}
				p[iy][ix] = s;
			}
		}
		free(u);
	}

	if (nr_lines >= 2) {
		creal8 * u = (creal8 *)malloc(sizeof(creal8) * nr_lines);
		if (u == NULL) return RET_ERROR;

		for (int ix = 0; ix < nr_pixels; ix++) {
			for (int iy = 0; iy < nr_lines; iy ++) {
				u[iy] = p[iy][ix];
			}

			for (int iy = 0; iy < nr_lines; iy ++) {
				creal8 s(0.0, 0.0);
				creal8 w = creal8(0.0, (double)dir) * (2.0 * PI / nr_lines) * (double)iy;
				for (int j = 0; j < nr_lines; j ++) {
					s += u[j] * exp(w * (double)j);
				}
				p[iy][ix] = s;
			}
		}
		free(u);
	}

	complex_type = (dir == -1) ? COMPLEX_FREQ : COMPLEX_SPACE;
	return RET_OK;
}

//
// Fast Fourier Transform
//
int fourier::Forward_FFT(
	void)
{
	return FFT(FORWARD_FT);
}

int fourier::Backward_FFT(
	void)
{
	return FFT(BACKWARD_FT);
}

int fourier::FFT(
	int	dir)		// dir = -1  ==> Forward  DFT
					// dir = +1  ==> Backward DFT
{
	if (nr_pixels >= 2) {
		for (int iy = 0; iy < nr_lines; iy++) {
			fft_1d(
				dir,
				p[iy],
				nr_pixels,
				w_pixel);
		}
	}

	if (nr_lines >= 2) {
		creal8 * u = (creal8 *)malloc(sizeof(creal8) * nr_lines);
		if (u == NULL) return RET_ERROR;

		for (int ix = 0; ix < nr_pixels; ix++) {
			for (int iy = 0; iy < nr_lines; iy ++) {
				u[iy] = p[iy][ix];
			}

			fft_1d(
				dir,
				u,
				nr_lines,
				w_line);

			for (int iy = 0; iy < nr_lines; iy ++) {
				p[iy][ix] = u[iy];
			}

		}
		free(u);
	}

	complex_type = (dir == -1) ? COMPLEX_FREQ : COMPLEX_SPACE;
	return RET_OK;
}


void fourier::fft_1d(
	int		dir,
	creal8	* data,
	int		fft_len,
	creal8	* w)
{
	// Shuffling
	for (int i = 0, j = 0; i < fft_len - 1; i ++) {
		if (i < j) {
			creal8 tmp     = data[j];
			       data[j] = data[i];
			       data[i] = tmp;
		}
		
		int k = fft_len / 2;
		if (k < j + 1) {
			do {
				j -= k;
				k /= 2;
			} while (k < j + 1);
		}
		j += k;
	}

	int power_of_2 = (int)(log((double)fft_len) / log(2.0) + 0.5);

	for (int i = 1; i <= power_of_2; i ++) {	// Butterfly computing
		int k0 = 1 << i;
		int k1 = k0 / 2;
		int kk = 0;

		for (int j = 0; j < k1; j ++) {
			for (int k = j; k < fft_len; k += k0) {
				creal8 tmp = data[k + k1] * ((dir == FORWARD_FT) ? w[kk] : conj(w[kk]));
				data[k + k1] = data[k] - tmp;
				data[k     ] = data[k] + tmp;
			}
			kk += (fft_len / k0);
		}
	}
}

// Load complex from "Complex Text" file
void fourier::cpx_load(void)
{
	char fn[FFT_MAX_FILENAME_LEN];
	fprintf(fft_stdout, "\nLoad complex data from \"Complex Text\" format file.\n");
	fprintf(fft_stdout, "Input file name including file extension(ex. cpx, txt) : ");
	fgets(fn, FFT_MAX_FILENAME_LEN - 1, fft_stdin);
	if (fn[strlen(fn) - 1] == 0x0d ||
	    fn[strlen(fn) - 1] == 0x0a ) fn[strlen(fn) - 1] = 0;
	if (fn[strlen(fn) - 2] == 0x0d ||
	    fn[strlen(fn) - 2] == 0x0a ) fn[strlen(fn) - 2] = 0;
	cpx_load(fn);
}

void fourier::cpx_load(
	char * fn)
{
	if (fn == NULL || fn[0] == 0) {
		cpx_load();
		return;
	}

	unalloc();

	FILE * fp_cpx;
	fp_cpx = fopen(fn, "r");
	if (fp_cpx == NULL) {
		fprintf(fft_stderr, "Complex Text file: %s cannot open.\n", fn);
		return;
	}

	fprintf(fft_stdout, "\nLoading complex data from Complex Text file: %s\n", fn);

	int nr_lines;	// Number of lines
	int nr_pixels;	// Number of pixels

	char * src_cpx_id = (char *)malloc(20);
	char * src_fmt_id = (char *)malloc(20);
	fscanf(fp_cpx, "%19s %19s %d %d%*[^\n]%*c", src_cpx_id, src_fmt_id, &nr_lines, &nr_pixels);
	if (strcmp(src_cpx_id, COMPLEX_TEXT_ID) != 0 ||
	       (strcmp(src_fmt_id, COMPLEX_TEXT_SPACE) != 0 &&
	        strcmp(src_fmt_id, COMPLEX_TEXT_FREQ ) != 0)    ) {
		fprintf(fft_stderr, "File %s has invalid \"Complex Text\" Format.\n", fn);
		fprintf(fft_stderr, "1:%s %s %d %d\n", src_cpx_id, src_fmt_id, nr_lines, nr_pixels);
		return;
	}

	if (alloc(nr_lines, nr_pixels) == 0) {
		fprintf(fft_stderr, "Memory allocation error in loading Complex Text file\n");
	} else {
		fprintf(fft_stdout, "Content of complex data : %s\n", src_fmt_id);
		fprintf(fft_stdout, "Number of lines  : %8d\n", nr_lines);
		fprintf(fft_stdout, "Number of pixels : %8d\n", nr_pixels);

		int dum_iy;
		fscanf(fp_cpx, "%d %*[^\n]%*c", &dum_iy);

		for (int iy = 0; iy < nr_lines; iy ++) {
			fscanf(fp_cpx, "%d", &dum_iy);

			for (int ix = 0; ix < nr_pixels; ix ++) {
				double re, im;
				fscanf(fp_cpx, "%le %le", &re, &im);
				p[iy][ix] = creal8(re, im);
			}
			fscanf(fp_cpx, "%*[^\n]%*c");
		}
	}

	if (strcmp(src_fmt_id, COMPLEX_TEXT_SPACE) == 0) {
		complex_type = COMPLEX_SPACE;
	} else {
		complex_type = COMPLEX_FREQ;
	}

	fclose (fp_cpx);
}

// Save complex to "Complex Text" file
void fourier::cpx_save(void)
{
	if (p == NULL) {
		fprintf(fft_stderr, "\nNo complex data in this class\n");
		return;
	}

	char fn[FFT_MAX_FILENAME_LEN];
	fprintf(fft_stdout, "\nSave complex data (%s) to \"Complex Text\" format file.\n",
		(complex_type == COMPLEX_SPACE) ? COMPLEX_TEXT_SPACE : COMPLEX_TEXT_FREQ);
	fprintf(fft_stdout, "Input file name including file extension(ex. cpx, txt) : ");
	fgets(fn, FFT_MAX_FILENAME_LEN - 1, fft_stdin);
	if (fn[strlen(fn) - 1] == 0x0d ||
	    fn[strlen(fn) - 1] == 0x0a ) fn[strlen(fn) - 1] = 0;
	if (fn[strlen(fn) - 2] == 0x0d ||
	    fn[strlen(fn) - 2] == 0x0a ) fn[strlen(fn) - 2] = 0;
	cpx_save(fn);
}

void fourier::cpx_save(
	char * fn)
{
	if (p == NULL) {
		fprintf(fft_stderr, "\nNo complex data in this class\n");
		return;
	}
	if (fn == NULL || fn[0] == 0) {
		cpx_save();
		return;
	}

	FILE * fp_cpx;
	fp_cpx = fopen(fn, "w");
	if (fp_cpx == NULL) {
		fprintf(fft_stderr, "Complex Text file: %s cannot open.\n", fn);
		return;
	}

	fprintf(fft_stdout, "\nSaving complex data to Complex Text file: %s\n", fn);
	fprintf(fft_stdout, "Content of complex data : %s\n",
		(complex_type == COMPLEX_SPACE) ? COMPLEX_TEXT_SPACE : COMPLEX_TEXT_FREQ);
	fprintf(fft_stdout, "Number of lines  : %8d\n", nr_lines);
	fprintf(fft_stdout, "Number of pixels : %8d\n", nr_pixels);

	fprintf(fp_cpx, "%s %s %d %d\n", COMPLEX_TEXT_ID,
		(complex_type == COMPLEX_SPACE) ? COMPLEX_TEXT_SPACE : COMPLEX_TEXT_FREQ,
		nr_lines, nr_pixels);

	fprintf(fp_cpx, "        ");
	for (int ix = 0; ix < nr_pixels; ix ++) {
		fprintf(fp_cpx,"%6d     Re             Im   ", ix);
	}
	fprintf(fp_cpx, "\n");

	for (int iy = 0; iy < nr_lines; iy ++) {
		fprintf(fp_cpx, "%8d", iy);

		for (int ix = 0; ix < nr_pixels; ix ++) {
			fprintf(fp_cpx, " %+14.6le %+14.6le", real(p[iy][ix]), imag(p[iy][ix]));
		}
		fprintf(fp_cpx, "\n");
	}

	fclose (fp_cpx);
}

const int OUTPUT_REAL                 =  0;
const int OUTPUT_IMAGE                =  1;
const int OUTPUT_AMP                  =  2;
const int OUTPUT_POWER                =  3;
const int OUTPUT_DB                   =  4;
const int OUTPUT_PHASE                =  5;

const int FREQ_SPACE_SHUFFLING_OFF    =  0;
const int FREQ_SPACE_SHUFFLING_ON     =  1;

const int OUTPUT_SCALE_OFF            =  0;
const int OUTPUT_SCALE_ON             =  1;

void fourier::MakeImage(
	int output_data,
	int	shuffling,
	int	scale,
	unsigned char ** g)
{
	fprintf(fft_stdout, "Making image...\n");
	fprintf(fft_stdout, "Output value is ");
	     if (output_data == OUTPUT_REAL)  fprintf(fft_stdout, "Real part\n");
	else if (output_data == OUTPUT_IMAGE) fprintf(fft_stdout, "Imaginary part\n");
	else if (output_data == OUTPUT_AMP)   fprintf(fft_stdout, "Amplitude\n");
	else if (output_data == OUTPUT_POWER) fprintf(fft_stdout, "Power\n");
	else if (output_data == OUTPUT_DB)    fprintf(fft_stdout, "Power [dB]\n");
	else if (output_data == OUTPUT_PHASE) fprintf(fft_stdout, "Phase\n");

	     if (shuffling == FREQ_SPACE_SHUFFLING_OFF) fprintf(fft_stdout, "No shuffling image\n");
	else if (shuffling == FREQ_SPACE_SHUFFLING_ON)  fprintf(fft_stdout, "Shuffling image\n");

	     if (scale == OUTPUT_SCALE_OFF) fprintf(fft_stdout, "No scaling image\n");
	else if (scale == OUTPUT_SCALE_ON)  fprintf(fft_stdout, "Automatic scaling image\n");

	double min_g, max_g;
	if (scale == OUTPUT_SCALE_ON) {
		if (output_data == OUTPUT_PHASE) {
			min_g = - PI;
			max_g = + PI;
		} else {
			min_g = GetOutputData(output_data, p[0][0]);
			max_g = min_g;
			for (int iy = 0; iy < nr_lines; iy ++) {
				for (int ix = 0; ix < nr_pixels; ix ++) {
					double pp = GetOutputData(output_data, p[iy][ix]);
					if (pp < min_g) min_g = pp;
					if (pp > max_g) max_g = pp;
				}
			}
		}
		fprintf(fft_stdout, "Minimum value : %lf\n", min_g);
		fprintf(fft_stdout, "Maximum value : %lf\n", max_g);
	}

	for (int iy = 0; iy < nr_lines; iy ++) {
		for (int ix = 0; ix < nr_pixels; ix ++) {
			int px;
			double pp = GetOutputData(output_data, p[iy][ix]);

			if (scale == OUTPUT_SCALE_OFF) {
				if (pp >= (double)MAX_GRAY_LEVEL) pp = (double)MAX_GRAY_LEVEL;
				if (pp <= (double)MIN_GRAY_LEVEL) pp = (double)MIN_GRAY_LEVEL;
				px = (int)(pp + 0.5);

			} else {
				px = (int)((pp - min_g) / (max_g - min_g) * MAX_GRAY_LEVEL + 0.5);
				if (px >= MAX_GRAY_LEVEL) px = MAX_GRAY_LEVEL;
				if (px <= MIN_GRAY_LEVEL) px = MIN_GRAY_LEVEL;
			}

			int iiy = iy, iix = ix;
			if (shuffling == FREQ_SPACE_SHUFFLING_ON) {
				if (iy < nr_lines / 2) {
					iiy += (nr_lines / 2);
				} else {
					iiy -= (nr_lines / 2);
				}
				if (ix < nr_pixels / 2) {
					iix += (nr_pixels / 2);
				} else {
					iix -= (nr_pixels / 2);
				}
			}
			g[iiy][iix] = px;
		}
	}
}

double fourier::GetOutputData(
	int output_data,
	creal8 x)
{
	double r = 0.0;
	double re = real(x);
	double im = imag(x);
	if (output_data == OUTPUT_REAL) {
		r = re;
	} else if (output_data == OUTPUT_IMAGE) {
		r = im;
	} else if (output_data == OUTPUT_AMP) {
		r = sqrt(re * re + im * im);
	} else if (output_data == OUTPUT_POWER) {
		r = re * re + im * im;
	} else if (output_data == OUTPUT_DB) {
		if (re == 0.0 && im == 0.0) r = -100.0;
		else r = 10.0 * log10(re * re + im * im);
	} else if (output_data == OUTPUT_PHASE) {
		if (re == 0.0 && im == 0.0) r = 0.0;
		else r = arg(x);
	}
	return r;
}

#endif // FFT_HEADER_HH_INCLUDED