pip/libs/main/math/pifft.h

//! \addtogroup Math
//! \{
//! \file pifft.h
//! \brief
//! \~english Declares \a PIFFT classes
//! \~russian Объявление классов \a PIFFT
//! \~\authors
//! \~english Ivan Pelipenko peri4ko@yandex.ru; Andrey Bychkov work.a.b@yandex.ru
//! \~russian Иван Пелипенко peri4ko@yandex.ru; Андрей Бычков work.a.b@yandex.ru
//! \~\}
/*
    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/>.
*/
//! \defgroup FFTW FFTW
//! \~\brief
//! \~english Optimized FFT support via libfftw3
//! \~russian Оптимизированный БПФ с помощью libfftw3
//!
//! \~\details
//! \~english \section cmake_module_FFTW Building with CMake
//! \~russian \section cmake_module_FFTW Сборка с использованием CMake
//!
//! \~\code
//! find_package(PIP REQUIRED)
//! target_link_libraries([target] PIP::FFTW)
//! \endcode
//!
//! \~english \par Common
//! \~russian \par Общее
//!
//! \~english
//! These files provides FFT using [libfftw3](https://fftw.org/)
//!
//! \~russian
//! Эти файлы обеспечивают БПФ с использованием [libfftw3](https://fftw.org/)
//!
//! \~\authors
//! \~english
//! Ivan Pelipenko peri4ko@yandex.ru;
//! Andrey Bychkov work.a.b@yandex.ru;
//! \~russian
//! Иван Пелипенко peri4ko@yandex.ru;
//! Андрей Бычков work.a.b@yandex.ru;
//!

#ifndef PIFFT_H
#define PIFFT_H

#include "pimathcomplex.h"

#ifndef MICRO_PIP

#  include "pip_fftw_export.h"

//! \addtogroup Math
//! \{
//! \class PIFFT_double
//! \brief
//! \~english Fast Fourier Transform implementation for double precision.
//! \~russian Реализация быстрого преобразования Фурье для двойной точности.
//! \~\}
//! \sa \a PIFFT_float, \a PIFFTW
class PIP_EXPORT PIFFT_double {
public:
	//! \~english Default constructor.
	//! \~russian Конструктор по умолчанию.
	PIFFT_double();

	//! \~english Calculate FFT from complex vector.
	//! \~russian Вычисление БПФ из комплексного вектора.
	PIVector<complexd> * calcFFT(const PIVector<complexd> & val);
	//! \~english Calculate FFT from real vector.
	//! \~russian Вычисление БПФ из вещественного вектора.
	PIVector<complexd> * calcFFT(const PIVector<double> & val);
	//! \~english Calculate inverse FFT.
	//! \~russian Вычисление обратного БПФ.
	PIVector<complexd> * calcFFTinverse(const PIVector<complexd> & val);
	//! \~english Calculate Hilbert transform.
	//! \~russian Вычисление преобразования Гильберта.
	PIVector<complexd> * calcHilbert(const PIVector<double> & val);
	//! \~english Get amplitude spectrum.
	//! \~russian Получить амплитудный спектр.
	PIVector<double> getAmplitude() const;
	//! \~english Get real part.
	//! \~russian Получить действительную часть.
	PIVector<double> getReal() const;
	//! \~english Get imaginary part.
	//! \~russian Получить мнимую часть.
	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);
};

//! \addtogroup Math
//! \{
//! \class PIFFT_float
//! \brief
//! \~english Fast Fourier Transform implementation for single precision.
//! \~russian Реализация быстрого преобразования Фурье для одинарной точности.
//! \~\}
//! \sa \a PIFFT_double, \a PIFFTW
class PIP_EXPORT PIFFT_float {
public:
	//! \~english Default constructor.
	//! \~russian Конструктор по умолчанию.
	PIFFT_float();

	//! \~english Calculate FFT from complex vector.
	//! \~russian Вычисление БПФ из комплексного вектора.
	PIVector<complexf> * calcFFT(const PIVector<complexf> & val);
	//! \~english Calculate FFT from real vector.
	//! \~russian Вычисление БПФ из вещественного вектора.
	PIVector<complexf> * calcFFT(const PIVector<float> & val);
	//! \~english Calculate inverse FFT.
	//! \~russian Вычисление обратного БПФ.
	PIVector<complexf> * calcFFTinverse(const PIVector<complexf> & val);
	//! \~english Calculate Hilbert transform.
	//! \~russian Вычисление преобразования Гильберта.
	PIVector<complexf> * calcHilbert(const PIVector<float> & val);
	//! \~english Get amplitude spectrum.
	//! \~russian Получить амплитудный спектр.
	PIVector<float> getAmplitude() const;
	//! \~english Get real part.
	//! \~russian Получить действительную часть.
	PIVector<float> getReal() const;
	//! \~english Get imaginary part.
	//! \~russian Получить мнимую часть.
	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)

//! \addtogroup FFTW
//! \{
//! \class PIFFTW
//! \brief
//! \~english FFTW wrapper for arbitrary precision types.
//! \~russian Обёртка FFTW для типов произвольной точности.
//! \~\}
//! \note Requires linking against libfftw3
//! \sa \a PIFFT_double, \a PIFFT_float
template<typename T>
class PIFFTW {
public:
	//! \~english Default constructor.
	//! \~russian Конструктор по умолчанию.
	explicit PIFFTW() {
		p = 0;
		newP(p);
	}
	//! \~english Destructor.
	//! \~russian Деструктор.
	~PIFFTW() { deleteP(p); }

	//! \~english Calculate FFT from complex vector.
	//! \~russian Вычисление БПФ из комплексного вектора.
	inline const PIVector<complex<T>> & calcFFT(const PIVector<complex<T>> & in) { return PIVector<complex<T>>().resize(in.size()); }
	//! \~english Calculate FFT from real vector.
	//! \~russian Вычисление БПФ из вещественного вектора.
	inline const PIVector<complex<T>> & calcFFT(const PIVector<T> & in) { return PIVector<complex<T>>().resize(in.size()); }
	//! \~english Calculate inverse FFT.
	//! \~russian Вычисление обратного БПФ.
	inline const PIVector<complex<T>> & calcFFTinverse(const PIVector<complex<T>> & in) { return PIVector<complex<T>>().resize(in.size()); }

	//! \~english FFT operation type.
	//! \~russian Тип операции БПФ.
	enum FFT_Operation {
		foReal,
		foComplex,
		foInverse
	};

	//! \~english Prepare computation plan.
	//! \~russian Подготовить план вычислений.
	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 // MICRO_PIP

#endif // PIFFT_H