/*! \file pifft.h
 * \brief Class for FFT, IFFT and Hilbert transformations
*/
/*
	PIP - Platform Independent Primitives
	Class for FFT, IFFT and Hilbert transformations
	Andrey Bychkov work.a.b@yandex.ru, Ivan Pelipenko peri4ko@yandex.ru

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU Lesser General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU Lesser General Public License for more details.

	You should have received a copy of the GNU Lesser General Public License
	along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef PIFFT_H
#define PIFFT_H

#include "pip_fftw_export.h"
#include "pimathcomplex.h"

class PIP_EXPORT PIFFT_double
{
public:
	PIFFT_double();
	
	PIVector<complexd> * calcFFT(const PIVector<complexd> &val);
	PIVector<complexd> * calcFFT(const PIVector<double> &val);
	PIVector<complexd> * calcFFTinverse(const PIVector<complexd> &val);
	PIVector<complexd> * calcHilbert(const PIVector<double> &val);
	PIVector<double> getAmplitude() const;
	PIVector<double> getReal() const;
	PIVector<double> getImag() const;
	
private:
	PIVector<complexd> result;
	typedef ptrdiff_t ae_int_t;

	struct ftplan {
		PIVector<int> plan;
		PIVector<double> precomputed;
		PIVector<double> tmpbuf;
		PIVector<double> stackbuf;
	};

	ftplan curplan;

	void fftc1d(const PIVector<complexd> &a, uint n);
	void fftc1r(const PIVector<double> &a, uint n);
	void fftc1dinv(const PIVector<complexd> &a, uint n);

	void createPlan(uint n);
	void ftbasegeneratecomplexfftplan(uint n, ftplan *plan);
	void ftbase_ftbasegenerateplanrec(int n, int tasktype, ftplan *plan, int *plansize, int *precomputedsize, int *planarraysize, int *tmpmemsize, int *stackmemsize, ae_int_t stackptr, int debugi=0);
	void ftbase_ftbaseprecomputeplanrec(ftplan *plan, int entryoffset, ae_int_t stackptr);
	void ftbasefactorize(int n, int *n1, int *n2);
	void ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int *best);
	int ftbasefindsmooth(int n);
	void ftbaseexecuteplan(PIVector<double> *a, int aoffset, int n, ftplan *plan);
	void ftbaseexecuteplanrec(PIVector<double> *a, int aoffset, ftplan *plan, int entryoffset, ae_int_t stackptr);
	void ftbase_internalcomplexlintranspose(PIVector<double> *a, int m, int n, int astart, PIVector<double> *buf);
	void ftbase_ffticltrec(PIVector<double> *a, int astart, int astride, PIVector<double> *b, int bstart, int bstride, int m, int n);
	void ftbase_internalreallintranspose(PIVector<double> *a, int m, int n, int astart, PIVector<double> *buf);
	void ftbase_fftirltrec(PIVector<double> *a, int astart, int astride, PIVector<double> *b, int bstart, int bstride, int m, int n);
	void ftbase_ffttwcalc(PIVector<double> *a, int aoffset, int n1, int n2);
};

class PIP_EXPORT PIFFT_float
{
public:
	PIFFT_float();

	PIVector<complexf> * calcFFT(const PIVector<complexf> &val);
	PIVector<complexf> * calcFFT(const PIVector<float> &val);
	PIVector<complexf> * calcFFTinverse(const PIVector<complexf> &val);
	PIVector<complexf> * calcHilbert(const PIVector<float> &val);
	PIVector<float> getAmplitude() const;
	PIVector<float> getReal() const;
	PIVector<float> getImag() const;

private:
	PIVector<complexf> result;
	typedef ptrdiff_t ae_int_t;

	struct ftplan {
		PIVector<int> plan;
		PIVector<float> precomputed;
		PIVector<float> tmpbuf;
		PIVector<float> stackbuf;
	};

	ftplan curplan;

	void fftc1d(const PIVector<complexf> &a, uint n);
	void fftc1r(const PIVector<float> &a, uint n);
	void fftc1dinv(const PIVector<complexf> &a, uint n);

	void createPlan(uint n);
	void ftbasegeneratecomplexfftplan(uint n, ftplan *plan);
	void ftbase_ftbasegenerateplanrec(int n, int tasktype, ftplan *plan, int *plansize, int *precomputedsize, int *planarraysize, int *tmpmemsize, int *stackmemsize, ae_int_t stackptr, int debugi=0);
	void ftbase_ftbaseprecomputeplanrec(ftplan *plan, int entryoffset, ae_int_t stackptr);
	void ftbasefactorize(int n, int *n1, int *n2);
	void ftbase_ftbasefindsmoothrec(int n, int seed, int leastfactor, int *best);
	int ftbasefindsmooth(int n);
	void ftbaseexecuteplan(PIVector<float> *a, int aoffset, int n, ftplan *plan);
	void ftbaseexecuteplanrec(PIVector<float> *a, int aoffset, ftplan *plan, int entryoffset, ae_int_t stackptr);
	void ftbase_internalcomplexlintranspose(PIVector<float> *a, int m, int n, int astart, PIVector<float> *buf);
	void ftbase_ffticltrec(PIVector<float> *a, int astart, int astride, PIVector<float> *b, int bstart, int bstride, int m, int n);
	void ftbase_internalreallintranspose(PIVector<float> *a, int m, int n, int astart, PIVector<float> *buf);
	void ftbase_fftirltrec(PIVector<float> *a, int astart, int astride, PIVector<float> *b, int bstart, int bstride, int m, int n);
	void ftbase_ffttwcalc(PIVector<float> *a, int aoffset, int n1, int n2);
};

typedef PIFFT_double PIFFT;
typedef PIFFT_double PIFFTd;
typedef PIFFT_float PIFFTf;

#ifndef CC_VC

#define _PIFFTW_H(type) class PIP_FFTW_EXPORT _PIFFTW_P_##type##_ { \
public: \
	_PIFFTW_P_##type##_(); \
	~_PIFFTW_P_##type##_(); \
	const PIVector<complex<type> > & calcFFT(const PIVector<complex<type> > & in); \
	const PIVector<complex<type> > & calcFFTR(const PIVector<type> & in); \
	const PIVector<complex<type> > & calcFFTI(const PIVector<complex<type> > & in); \
	void preparePlan(int size, int op); \
	void * impl; \
};
_PIFFTW_H(float)
_PIFFTW_H(double)
_PIFFTW_H(ldouble)

template <typename T>
class PIFFTW
{
public:
	explicit PIFFTW() {p = 0; newP(p);}
	~PIFFTW() {deleteP(p);}

	inline const PIVector<complex<T> > & calcFFT(const PIVector<complex<T> > & in) {return PIVector<complex<T> >().resize(in.size());}
	inline const PIVector<complex<T> > & calcFFT(const PIVector<T> & in) {return PIVector<complex<T> >().resize(in.size());}
	inline const PIVector<complex<T> > & calcFFTinverse(const PIVector<complex<T> > & in) {return PIVector<complex<T> >().resize(in.size());}

	enum FFT_Operation {foReal, foComplex, foInverse};

	inline void preparePlan(int size, FFT_Operation op) {}

private:
	void operator =(const PIFFTW & );
	PIFFTW(const PIFFTW &);
	inline void newP(void *& _p) {}
	inline void deleteP(void *& _p) {}

	void * p;
};


template<> inline const PIVector<complex<float> > & PIFFTW<float>::calcFFT(const PIVector<complex<float> > & in) {return ((_PIFFTW_P_float_*)p)->calcFFT(in);}
template<> inline const PIVector<complex<float> > & PIFFTW<float>::calcFFT(const PIVector<float> & in) {return ((_PIFFTW_P_float_*)p)->calcFFTR(in);}
template<> inline const PIVector<complex<float> > & PIFFTW<float>::calcFFTinverse(const PIVector<complex<float> > & in) {return ((_PIFFTW_P_float_*)p)->calcFFTI(in);}
template<> inline void PIFFTW<float>::preparePlan(int size, FFT_Operation op) {((_PIFFTW_P_float_*)p)->preparePlan(size, op);}
template<> inline void PIFFTW<float>::newP(void *& _p) {_p = new _PIFFTW_P_float_();}
template<> inline void PIFFTW<float>::deleteP(void *& _p) {if (_p) delete (_PIFFTW_P_float_*)_p; _p = 0;}

typedef PIFFTW<float>   PIFFTWf;


template<> inline const PIVector<complex<double> > & PIFFTW<double>::calcFFT(const PIVector<complex<double> > & in) {return ((_PIFFTW_P_double_*)p)->calcFFT(in);}
template<> inline const PIVector<complex<double> > & PIFFTW<double>::calcFFT(const PIVector<double> & in) {return ((_PIFFTW_P_double_*)p)->calcFFTR(in);}
template<> inline const PIVector<complex<double> > & PIFFTW<double>::calcFFTinverse(const PIVector<complex<double> > & in) {return ((_PIFFTW_P_double_*)p)->calcFFTI(in);}
template<> inline void PIFFTW<double>::preparePlan(int size, FFT_Operation op) {((_PIFFTW_P_double_*)p)->preparePlan(size, op);}
template<> inline void PIFFTW<double>::newP(void *& _p) {_p = new _PIFFTW_P_double_();}
template<> inline void PIFFTW<double>::deleteP(void *& _p) {if (_p) delete (_PIFFTW_P_double_*)_p; _p = 0;}

typedef PIFFTW<double>  PIFFTWd;


template<> inline const PIVector<complex<ldouble> > & PIFFTW<ldouble>::calcFFT(const PIVector<complex<ldouble> > & in) {return ((_PIFFTW_P_ldouble_*)p)->calcFFT(in);}
template<> inline const PIVector<complex<ldouble> > & PIFFTW<ldouble>::calcFFT(const PIVector<ldouble> & in) {return ((_PIFFTW_P_ldouble_*)p)->calcFFTR(in);}
template<> inline const PIVector<complex<ldouble> > & PIFFTW<ldouble>::calcFFTinverse(const PIVector<complex<ldouble> > & in) {return ((_PIFFTW_P_ldouble_*)p)->calcFFTI(in);}
template<> inline void PIFFTW<ldouble>::preparePlan(int size, FFT_Operation op) {((_PIFFTW_P_ldouble_*)p)->preparePlan(size, op);}
template<> inline void PIFFTW<ldouble>::newP(void *& _p) {_p = new _PIFFTW_P_ldouble_();}
template<> inline void PIFFTW<ldouble>::deleteP(void *& _p) {if (_p) delete (_PIFFTW_P_ldouble_*)_p; _p = 0;}

typedef PIFFTW<ldouble> PIFFTWld;

#endif

#endif // PIFFT_H