pip/libs/fftw/pifftw_p.h

/*! \file pifftw_p.h
 * \brief Class for FFT, IFFT and Hilbert transformations
 */
/*
    PIP - Platform Independent Primitives
    Private header for fftw3
    Ivan Pelipenko peri4ko@yandex.ru, Andrey Bychkov work.a.b@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 PIFFTW_P_H
#define PIFFTW_P_H

#include "pimathcomplex.h"
#include "pimutex.h"
#include "pip_fftw_export.h"
#include "pivector.h"

template<typename T>
class PIP_FFTW_EXPORT FFTWInterface {
public:
	void * createPlan_c2c_1d(int size, const void * in, void * out, bool hard, bool forward) { return nullptr; }
	void * createPlan_r2c_1d(int size, const void * in, void * out, bool hard) { return nullptr; }
	void executePlan_c2c(void * plan, const void * in, void * out) {}
	void executePlan_r2c(void * plan, const void * in, void * out) {}
	void destroyPlan(void *& plan) {}
};


#define FFTW_DECLARATION(T)                                                                       \
	template<>                                                                                    \
	class PIP_FFTW_EXPORT FFTWInterface<T> {                                                      \
	public:                                                                                       \
		void * createPlan_c2c_1d(int size, const void * in, void * out, bool hard, bool forward); \
		void * createPlan_r2c_1d(int size, const void * in, void * out, bool hard);               \
		void executePlan_c2c(void * plan, const void * in, void * out);                           \
		void executePlan_r2c(void * plan, const void * in, void * out);                           \
		void destroyPlan(void *& plan);                                                           \
	};

FFTW_DECLARATION(float)
FFTW_DECLARATION(double)
FFTW_DECLARATION(ldouble)

#undef FFTW_DECLARATION


class PIFFTWCommon {
public:
	static PIFFTWCommon & instance();
	PIMutex mutex_plan;

private:
	PIFFTWCommon() {}
};


template<typename T>
class FFTW {
public:
	~FFTW() {
		for (int i = 0; i < oOperationsCount; ++i)
			plans[i].destroy(fftw);
	}

	enum Operation {
		oForwardComplex,
		oForwardReal,
		oBackwardComplex,
		oOperationsCount
	};

	struct Plan {
		bool create(FFTWInterface<T> & fftw, size_t size_, int op_, bool hard_) {
			if ((size_ == size) && (op_ == op) && (hard || !hard_)) return true;
			destroy(fftw);
			PIVector<T> p_inr;
			PIVector<complex<T>> p_in, p_out;
			p_out.resize(size_);
			PIFFTWCommon::instance().mutex_plan.lock(); // fftw_plan is not thread-safe
			switch ((Operation)op_) {
			case oForwardComplex:
				p_in.resize(size_);
				plan_ptr = fftw.createPlan_c2c_1d(size_, p_in.data(), p_out.data(), hard_, true);
				break;
			case oForwardReal:
				p_inr.resize(size_);
				plan_ptr = fftw.createPlan_r2c_1d(size_, p_inr.data(), p_out.data(), hard_);
				break;
			case oBackwardComplex:
				p_in.resize(size_);
				plan_ptr = fftw.createPlan_c2c_1d(size_, p_in.data(), p_out.data(), hard_, false);
				break;
			default: break;
			}
			PIFFTWCommon::instance().mutex_plan.unlock();
			if (plan_ptr) {
				size = size_;
				op   = op_;
				hard = hard_;
			}
			return plan_ptr;
		}
		void destroy(FFTWInterface<T> & fftw) {
			fftw.destroyPlan(plan_ptr);
			*this = {};
		}
		size_t size     = 0;
		int op          = oForwardComplex;
		bool hard       = false;
		void * plan_ptr = nullptr;
	};

	void calcFFT(const PIVector<complex<T>> & in, PIVector<complex<T>> & out) {
		Plan & plan(plans[oForwardComplex]);
		if (!plan.create(fftw, in.size(), oForwardComplex, false)) {
			out.clear();
			return;
		}
		out.resize(in.size());
		fftw.executePlan_c2c(plan.plan_ptr, in.data(), out.data());
	}

	void calcFFT(const PIVector<T> & in, PIVector<complex<T>> & out) {
		Plan & plan(plans[oForwardReal]);
		if (!plan.create(fftw, in.size(), oForwardReal, false)) {
			out.clear();
			return;
		}
		out.resize(in.size());
		fftw.executePlan_r2c(plan.plan_ptr, in.data(), out.data());
	}

	void calcFFTinverse(const PIVector<complex<T>> & in, PIVector<complex<T>> & out) {
		Plan & plan(plans[oBackwardComplex]);
		if (!plan.create(fftw, in.size(), oBackwardComplex, false)) {
			out.clear();
			return;
		}
		out.resize(in.size());
		fftw.executePlan_c2c(plan.plan_ptr, in.data(), out.data());
	}

	bool preparePlans(int size) {
		bool ret = true;
		for (int i = 0; i < oOperationsCount; ++i) {
			if (!plans[i].create(fftw, size, i, true)) ret = false;
		}
		return ret;
	}

private:
	Plan plans[oOperationsCount];
	FFTWInterface<T> fftw;
};


#endif // PIFFTW_P_H