#include "QFFTW.h"

// Usage:
// QFFTW FFT;
// FFT.DetermineFastestFFT(); // optinal 
// FFT.Initialize(FFTW_FORWARD,m_DataSize.x, m_DataSize.y, m_NumberOfThreads);
// or FFT.Initialize(FFTW_BACKWARD,m_DataSize.x, m_DataSize.y, m_NumberOfThreads);
// calculation
// FFT.FFT1D(Input,Output);	
// or FFT.FFT2D(Input,Output);
// or FFT.FFT3D(Input,Output);
// deconstruct
// FFT.Destroy();
//
// FFT can handle arrays of type 
// - vector<CComplex>
// - CComplex*
// - complex<double>*

vector<FFTW_wisdom_list> fftw_wisdom;

QFFTW::QFFTW(void)
{
	m_Init = false;
	m_FlagOptimze = FFTW_ESTIMATE;
}

QFFTW::~QFFTW(void)
{
	Destroy();
}

void QFFTW::Destroy()
{
	if (m_Init==true)
	{
		if (in != NULL){
			fftw_destroy_plan(plan);
			// do never use delete!
			fftw_free(in); 
			fftw_free(out);
			in = NULL;
			out = NULL;
			m_Init=false;			
		}
	}
}

void QFFTW::DetermineFastestFFT()
{
	m_FlagOptimze = FFTW_MEASURE;
}

void QFFTW::DisableDetermineFastestFFT()
{
	m_FlagOptimze = FFTW_ESTIMATE;
}


void QFFTW::Initialize(const int iDirection, const int SizeX, const int iThread)
{
	SetSize(SizeX, 0, 0);
	int N = SizeX;

	m_Dimension = 1;
	m_Direction = iDirection;
	m_Threads = iThread;
	// The data is an array of type fftw_complex, which is by default 
	// a double[2] composed of the real (in[i][0]) and 
	// imaginary (in[i][1]) parts of a complex number.		

	// Arrays dimensionieren
	in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
	out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

	int fftw_init_threads();	
	fftw_plan_with_nthreads(iThread); 

	// Direction: FFTW_FORWARD or FFTW_BACKWARD
	// FFTW_ESTIMATE - set up but dont test and optimize
	// FFTW_MEASURE - set up, run test and optimize, but that takes several seconds!
	
	loadWisdomIfAvailable();

	plan = fftw_plan_dft_1d(N, in, out, iDirection , m_FlagOptimze);

	m_Init=true;

}


void QFFTW::Initialize(const int iDirection, const int SizeX, const int SizeY, const int iThread)
{
	SetSize(SizeX, SizeY);
	int N = SizeX * SizeY;

	m_Dimension = 2;
	m_Direction = iDirection;
	m_Threads = iThread;
	// The data is an array of type fftw_complex, which is by default 
	// a double[2] composed of the real (in[i][0]) and 
	// imaginary (in[i][1]) parts of a complex number.		

	// Arrays dimensionieren
	in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
	out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

	int fftw_init_threads();	
	fftw_plan_with_nthreads(iThread); 

	// Direction: FFTW_FORWARD or FFTW_BACKWARD
	// FFTW_ESTIMATE - set up but dont test and optimize
	// FFTW_MEASURE - set up, run test and optimize, but that takes several seconds!
	
	loadWisdomIfAvailable();
	plan = fftw_plan_dft_2d(SizeX, SizeY, in, out, iDirection , m_FlagOptimze);

	m_Init=true;

}
void QFFTW::Initialize(const int iDirection, const int SizeX, const int SizeY, const int SizeZ, const int iThread)
{
	int N = SizeX * SizeY * SizeZ;

	m_Dimension = 3;
	m_Direction = iDirection;
	m_Threads = iThread;
	// The data is an array of type fftw_complex, which is by default 
	// a double[2] composed of the real (in[i][0]) and 
	// imaginary (in[i][1]) parts of a complex number.		

	// Arrays dimensionieren
	in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
	out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

	int fftw_init_threads();	
	fftw_plan_with_nthreads(iThread); 

	// Direction: FFTW_FORWARD or FFTW_BACKWARD
	// FFTW_ESTIMATE - set up but dont test and optimize
	// FFTW_MEASURE - set up, run test and optimize, but that takes several seconds!
	
	loadWisdomIfAvailable();
	plan = fftw_plan_dft_3d(SizeX, SizeY, SizeZ, in, out, iDirection , m_FlagOptimze);	

	m_Init=true;

}



void QFFTW::FFT1D(complex<double>* InArray, complex<double>* OutArray)
{		
	if (m_Init != true)
		return;		

	StartClock();
	memcpy(in,InArray, m_size.x  * sizeof(complex<double>));
	fftw_execute(plan); 
	memcpy(OutArray,out, m_size.x * sizeof(complex<double>));
	

	EndClock();
}


void QFFTW::FFT2D(complex<double>* InArray, complex<double>* OutArray)
{		
	if (m_Init != true)
		return;		

	StartClock();
	memcpy(in,InArray, m_size.x * m_size.y * sizeof(complex<double>));
	fftw_execute(plan); 
	memcpy(OutArray,out, m_size.x * m_size.y * sizeof(complex<double>));
		
	EndClock();
}

void QFFTW::FFT3D(complex<double>* InArray, complex<double>* OutArray)
{		
	if (m_Init != true)
		return;		

	StartClock();
	memcpy(in,InArray, m_size.x * m_size.y * m_size.z * sizeof(complex<double>));
	fftw_execute(plan); 
	memcpy(OutArray,out, m_size.x * m_size.y * m_size.z * sizeof(complex<double>));

	EndClock();
}


inline int QFFTW::ArrPos(const int x, const int y)
{
	return y + m_size.y * x;
}

inline int QFFTW::ArrPos(const int x, const int y, const int z)
{
	return z + m_size.z * (y + m_size.y * x);
}

void QFFTW::scaleAmplitude(const double amplitude , complex<double> * OutArray)
{		
	
	double CurrentAmplitude = TotalAmplitude(OutArray);
	double changefactor = amplitude/CurrentAmplitude;
	if (CurrentAmplitude==0)
		return;
	int size = m_size.y * m_size.x;
	for (int i=0; i < size; i++) 
	{			
			OutArray[i] = OutArray[i] * changefactor;
	}
}

double QFFTW::TotalAmplitude(complex<double> * OutArray)
{	
	double Amplitude = 0;
	int size = m_size.y * m_size.x;
	for (int i=0; i < size; i++) 
	{
		Amplitude = Amplitude  + abs(OutArray[i]);
	}
	return Amplitude;
}

void QFFTW::SetSize(const int x)
{
	m_size.x = x;
	m_size.y = 0;
	m_size.z = 0;
}

void QFFTW::SetSize(const int x, const int y)
{
	m_size.x = x;
	m_size.y = y;
	m_size.z = 0;
}
void QFFTW::SetSize(const int x, const int y, const int z)
{
	m_size.x = x;
	m_size.y = y;
	m_size.z = z;
}
void QFFTW::StartClock()
{
	start_clock=clock();
}
void QFFTW::EndClock()
{
	end_clock=clock();
}
double QFFTW::GetSpeed()
{
	return (double)(end_clock - start_clock) / CLOCKS_PER_SEC;
}

// *********************************
// Cycle changes 012 345 to 345 012
// *********************************
void QFFTW::Cycle(complex<double> *A1,long N)
{
	Cycle(A1,N,N/2);
}


void QFFTW::Cycle(complex<double> *A1,long N,long del)
{
long i;
complex<double> *D;

	D=new complex<double>[N];

	if (del>=N) del-=N;
	if (del<=-N) del+=N;

	if(del>0 && del<N)
	{
		for (i=N-del;i<N;i++)	D[i-N+del]=A1[i];
		for (i=N-1;i>=del;i--) A1[i]=A1[i-del];
		for (i=0;i<del;i++)	A1[i]=D[i];
	}
	else if(del<0	&&	del > -N)
	{
		del=-del;
		for (i=0;i<del;i++)	D[i]=A1[i];
		for (i=0;i<N-del;i++) A1[i]=A1[i+del];
		for (i=N-del;i<N;i++) A1[i]=D[i-N+del];
	}
	delete [] D;
}

// **************************************
// 1 2      4 3
// 3 4  ->  2 1
// **************************************
// A muss als complex<double> A[N][N] deklariert worden sein!!!
//  A[x][y] = A[x*N+y]    A[0][0]  ...  A[N][0]
//                          ...           ...
//                        A[0][N]  ...  A[N][N]
void QFFTW::Cycle2D(complex<double> *A,long N)  
{												
	complex<double> swap;							

	// Spalten shiften 
	for (int ix=0;ix<N;ix++)	
		Cycle(A+ix*N,N);	
	// Zeilen einzeln shiften
	for (int iy=0;iy<N;iy++)	
		for (int ix=0;ix<N/2;ix++)
		{
			swap=A[(ix+N/2)*N+iy];
			A[(ix+N/2)*N+iy]=A[ix*N+iy];
			A[ix*N+iy]=swap;
		}
}

void QFFTW::runWisdomCreation()
{
	for (unsigned int i = 0; i < fftw_wisdom.size(); i++)
	{
		DetermineFastestFFT();
		switch (fftw_wisdom[i].dimensions)
		{
		case 1:
			Initialize(fftw_wisdom[i].direction, fftw_wisdom[i].size_x, fftw_wisdom[i].threads);
			break;
		case 2:
			Initialize(fftw_wisdom[i].direction, fftw_wisdom[i].size_x,fftw_wisdom[i].size_y,  fftw_wisdom[i].threads);
			break;
		case 3:
			Initialize(fftw_wisdom[i].direction, fftw_wisdom[i].size_x,fftw_wisdom[i].size_y, fftw_wisdom[i].size_z, fftw_wisdom[i].threads);
			break;
		}		
		fftw_wisdom[i].wisdom = QString().fromAscii(fftw_export_wisdom_to_string());
		Destroy();
	}
}

void QFFTW::loadWisdomIfAvailable()
{
	for (unsigned int i = 0; i < fftw_wisdom.size(); i++) 
	{
		if (m_Dimension == fftw_wisdom[i].dimensions)
		{
			if ((m_size.x == fftw_wisdom[i].size_x)
				&& (m_size.y == fftw_wisdom[i].size_y)
				&& (m_size.z == fftw_wisdom[i].size_z)
				&& (m_Direction == fftw_wisdom[i].direction)			
				&& (m_Threads == fftw_wisdom[i].threads))
			{
				restoreWisdow(fftw_wisdom[i].wisdom);
			}
		}
	}
}

void QFFTW::restoreWisdow(QString Wisdom)
{
	int result = fftw_import_wisdom_from_string(Wisdom.toAscii());
}

void QFFTW::saveWisdowListToFile(QString Filename)
{
	QSettings settings(Filename, QSettings::IniFormat); 

	settings.beginGroup("FFTW Wisdom List");
		settings.beginWriteArray("WisdomList");
		for (unsigned int i = 0; i < fftw_wisdom.size(); i++) {
		 settings.setArrayIndex(i);
		 settings.setValue("Dimension", fftw_wisdom[i].dimensions);
		 settings.setValue("Direction", fftw_wisdom[i].direction);
		 settings.setValue("Size_X", fftw_wisdom[i].size_x);
		 settings.setValue("Size_Y", fftw_wisdom[i].size_y);
		 settings.setValue("Size_Z", fftw_wisdom[i].size_z);
		 settings.setValue("Threads", fftw_wisdom[i].threads);
		 settings.setValue("Wisdom", fftw_wisdom[i].wisdom);
		}
		settings.endArray();	
	settings.endGroup();
}

void QFFTW::readWisdowListFromFile(QString Filename)
{
	fftw_forget_wisdom();
	QSettings settings(Filename, QSettings::IniFormat); 
	settings.beginGroup("FFTW Wisdom List");
		fftw_wisdom.clear();
		int size = settings.beginReadArray("WisdomList");
		for (int i = 0; i < size; i++) {
			settings.setArrayIndex(i);
			FFTW_wisdom_list newWisdom;
			newWisdom.dimensions = settings.value("Dimension").toInt();
			newWisdom.direction= settings.value("Direction").toInt();
			newWisdom.size_x = settings.value("Size_X").toInt();
			newWisdom.size_y = settings.value("Size_Y").toInt();
			newWisdom.size_z = settings.value("Size_Z").toInt();
			newWisdom.threads = settings.value("Threads").toInt();
			newWisdom.wisdom = settings.value("Wisdom").toString();
			fftw_wisdom.push_back(newWisdom);
		}
		settings.endArray();				
    settings.endGroup();
}