diff --git a/CMakeLists.txt b/CMakeLists.txt index ae49efc4..1d8c8585 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -150,7 +150,7 @@ if (FFTW) if (WIN32) list(APPEND LIBS fftw3-3 fftw3f-3 fftw3l-3) else () - list(APPEND LIBS fftw3 fftw3f) + list(APPEND LIBS fftw3 fftw3f fftw3l) endif () else () message(STATUS "Building without fftw3 support") diff --git a/main.cpp b/main.cpp index ebf271aa..847c293d 100644 --- a/main.cpp +++ b/main.cpp @@ -31,7 +31,15 @@ const char pult_config[] = connectionmodel = AAAM2Xja7VXLTttAFD12QpsikKjUSixYlLbKEpIUtVIlVEfqhk2FWuiGRRolUYtoHgrmpYiv6IItf8AveMMH9E/YsG6Ph3sde5hGoQoblLGuPHfunTPjc49nADxDA110+LTYC7FrPCAPeAO+vZu+aX7c/8PGd45WCJC0OGcfT6FDnmfSTPtwhZFt3HjgDs/Qtu5jPbZHtI/x50XfIzMQbdwEolbg9INP4ku++myPaUtCHYRaT2j1ldIh3VP60/Qff8vSfXLu9BP6JX9K/0TVH6jqVe22P1X/fao/oddWu/paDs1vBf9Jv/EZ91clbyHqv7BL6sscDOd4v4WTqs6jzaHGJ8QJerxlpJSpdZ7IWFJvDW7I2JxZqIM62k6A57RZmMQGmlyrxdV+WGBnmR01mXPI267hBKwp4FeBeo9VPtssxyb7rzHg1B7T9nCMU45U8BZlWuVWtIcD/CRGOqtsbW09851tXsHN0UTlLIAdASjSXnLyLn+H7L2+xbGYvC63Ezqg543egkLmn8qnRF6USbM4Qp9godkhzI777Ne5bCIt/5UtGz2o/yGby0nKpjqmbOa1ynkjmyzIrzvIZUeBPjvlUmbh32EFJbGyJZhR8YcvlS+3TpjhqeWSyvUkpbI9plSWtcKLcsK05beOJVEnhaEFfHEH+RwpeMcpn1JKGqWMNOL+G6wZyahlpdVOtufKfbDS+guLke9O\n\ "; */ +#include "pifft.h" int main(int argc, char *argv[]) { + PIFFTWf fft; + PIVector in; + for (int i = 0; i < 32; ++i) + in << i%10; + PIVector out = fft.calcFFT(in); + piCout << out; + return 0; /*__S__ s, s1; s.text = "123"; s.i = -1; diff --git a/src/math/pifft.cpp b/src/math/pifft.cpp index 1e53fa61..9e3c487c 100644 --- a/src/math/pifft.cpp +++ b/src/math/pifft.cpp @@ -18,6 +18,7 @@ */ #include "pifft.h" +#include "pifft_p.h" PIFFT_double::PIFFT_double() { @@ -1892,3 +1893,17 @@ void PIFFT_float::ftbase_ffttwcalc(PIVector * a, int aoffset, int n1, int } } + + + +#define _PIFFTW_CPP(type) \ + _PIFFTW_P_##type##_::_PIFFTW_P_##type##_() {impl = new PIFFTW_Private();;} \ + _PIFFTW_P_##type##_::~_PIFFTW_P_##type##_() {delete (PIFFTW_Private*)impl;} \ + const PIVector > & _PIFFTW_P_##type##_::calcFFT(const PIVector > & in) {return ((PIFFTW_Private*)impl)->calcFFT(in);} \ + const PIVector > & _PIFFTW_P_##type##_::calcFFTR(const PIVector & in) {return ((PIFFTW_Private*)impl)->calcFFT(in);} \ + const PIVector > & _PIFFTW_P_##type##_::calcFFTI(const PIVector > & in) {return ((PIFFTW_Private*)impl)->calcFFTinverse(in);} \ + void _PIFFTW_P_##type##_::preparePlan(int size, int op) {return ((PIFFTW_Private*)impl)->preparePlan(size, op);} + +_PIFFTW_CPP(float) +_PIFFTW_CPP(double) +_PIFFTW_CPP(ldouble) diff --git a/src/math/pifft.h b/src/math/pifft.h index cee59c35..81614bec 100644 --- a/src/math/pifft.h +++ b/src/math/pifft.h @@ -122,74 +122,70 @@ typedef PIFFT_double PIFFTd; typedef PIFFT_float PIFFTf; - -// fftw3 realisation -#ifdef PIP_FFTW -#include -#include "fftw3.h" -#include "pimutex.h" - -static PIMutex __pip_fft_plan_mutex; +#define _PIFFTW_H(type) class _PIFFTW_P_##type##_ { \ +public: \ + _PIFFTW_P_##type##_(); \ + ~_PIFFTW_P_##type##_(); \ + const PIVector > & calcFFT(const PIVector > & in); \ + const PIVector > & calcFFTR(const PIVector & in); \ + const PIVector > & calcFFTI(const PIVector > & in); \ + void preparePlan(int size, int op); \ + void * impl; \ +}; +_PIFFTW_H(float) +_PIFFTW_H(double) +_PIFFTW_H(ldouble) template class PIP_EXPORT PIFFTW { public: - const PIVector > & calcFFT(const PIVector > &in) { - result.resize(in.size()); - __pip_fft_plan_mutex.lock(); - p = fftw_plan_dft_1d(in.size_s(), (fftw_complex *)in.data(), (fftw_complex *)result.data(), FFTW_FORWARD, FFTW_ESTIMATE); - __pip_fft_plan_mutex.unlock(); - fftw_execute(p); - fftw_destroy_plan(p); - return result; - } - const PIVector > & calcFFT(const PIVector &in) { - result.resize(in.size()); - __pip_fft_plan_mutex.lock(); - p = fftw_plan_dft_r2c_1d(in.size_s(), (double *)in.data(), (fftw_complex *)result.data(), FFTW_ESTIMATE); - __pip_fft_plan_mutex.unlock(); - fftw_execute(p); - fftw_destroy_plan(p); - return result; - } - const PIVector > & calcFFTinverse(const PIVector > &in) { - result.resize(in.size()); - __pip_fft_plan_mutex.lock(); - p = fftw_plan_dft_1d(in.size_s(), (fftw_complex *)in.data(), (fftw_complex *)result.data(), FFTW_BACKWARD, FFTW_ESTIMATE); - __pip_fft_plan_mutex.unlock(); - fftw_execute(p); - fftw_destroy_plan(p); - return result; - } + explicit PIFFTW() {p = 0; newP(p);} + ~PIFFTW() {deleteP(p);} + const PIVector > & calcFFT(const PIVector > & in) {return PIVector >().resize(in.size());} + const PIVector > & calcFFT(const PIVector & in) {return PIVector >().resize(in.size());} + const PIVector > & calcFFTinverse(const PIVector > & in) {return PIVector >().resize(in.size());} - enum FFT_Operation {fo_real, fo_complex, fo_inverse}; + enum FFT_Operation {foReal, foComplex, foInverse}; - void preparePlan(int size, FFT_Operation op) { - PIVector > in(size), out(size); - PIVector inr(size); - __pip_fft_plan_mutex.lock(); - switch (op) { - case fo_real: - p = fftw_plan_dft_r2c_1d(in.size_s(), (double *)inr.data(), (fftw_complex *)out.data(), FFTW_MEASURE); - break; - case fo_complex: - p = fftw_plan_dft_1d(in.size_s(), (fftw_complex *)in.data(), (fftw_complex *)out.data(), FFTW_FORWARD, FFTW_MEASURE); - break; - case fo_inverse: - p = fftw_plan_dft_1d(in.size_s(), (fftw_complex *)in.data(), (fftw_complex *)out.data(), FFTW_BACKWARD, FFTW_MEASURE); - break; - } - __pip_fft_plan_mutex.unlock(); - } + void preparePlan(int size, FFT_Operation op) {} private: - PIVector > result; - fftw_plan p; + void operator =(const PIFFTW & ); + inline void newP(void *& _p) {} + inline void deleteP(void *& _p) {} + + void * p; }; -typedef PIFFTW PIFFTWd; -#endif // PIP_FFTW +template<> const PIVector > & PIFFTW::calcFFT(const PIVector > & in) {return ((_PIFFTW_P_float_*)p)->calcFFT(in);} +template<> const PIVector > & PIFFTW::calcFFT(const PIVector & in) {return ((_PIFFTW_P_float_*)p)->calcFFTR(in);} +template<> const PIVector > & PIFFTW::calcFFTinverse(const PIVector > & in) {return ((_PIFFTW_P_float_*)p)->calcFFTI(in);} +template<> void PIFFTW::preparePlan(int size, FFT_Operation op) {((_PIFFTW_P_float_*)p)->preparePlan(size, op);} +template<> inline void PIFFTW::newP(void *& _p) {_p = new _PIFFTW_P_float_();} +template<> inline void PIFFTW::deleteP(void *& _p) {if (_p) delete (_PIFFTW_P_float_*)_p; _p = 0;} + + +template<> const PIVector > & PIFFTW::calcFFT(const PIVector > & in) {return ((_PIFFTW_P_double_*)p)->calcFFT(in);} +template<> const PIVector > & PIFFTW::calcFFT(const PIVector & in) {return ((_PIFFTW_P_double_*)p)->calcFFTR(in);} +template<> const PIVector > & PIFFTW::calcFFTinverse(const PIVector > & in) {return ((_PIFFTW_P_double_*)p)->calcFFTI(in);} +template<> void PIFFTW::preparePlan(int size, FFT_Operation op) {((_PIFFTW_P_double_*)p)->preparePlan(size, op);} +template<> inline void PIFFTW::newP(void *& _p) {_p = new _PIFFTW_P_double_();} +template<> inline void PIFFTW::deleteP(void *& _p) {if (_p) delete (_PIFFTW_P_double_*)_p; _p = 0;} + + +template<> const PIVector > & PIFFTW::calcFFT(const PIVector > & in) {return ((_PIFFTW_P_ldouble_*)p)->calcFFT(in);} +template<> const PIVector > & PIFFTW::calcFFT(const PIVector & in) {return ((_PIFFTW_P_ldouble_*)p)->calcFFTR(in);} +template<> const PIVector > & PIFFTW::calcFFTinverse(const PIVector > & in) {return ((_PIFFTW_P_ldouble_*)p)->calcFFTI(in);} +template<> void PIFFTW::preparePlan(int size, FFT_Operation op) {((_PIFFTW_P_ldouble_*)p)->preparePlan(size, op);} +template<> inline void PIFFTW::newP(void *& _p) {_p = new _PIFFTW_P_ldouble_();} +template<> inline void PIFFTW::deleteP(void *& _p) {if (_p) delete (_PIFFTW_P_ldouble_*)_p; _p = 0;} + + +typedef PIFFTW PIFFTWf; +typedef PIFFTW PIFFTWd; +typedef PIFFTW PIFFTWld; + #endif // PIFFT_H diff --git a/src/math/pifft_p.h b/src/math/pifft_p.h new file mode 100644 index 00000000..c43e53a9 --- /dev/null +++ b/src/math/pifft_p.h @@ -0,0 +1,143 @@ +/*! \file pifft_p.h + * \brief Class for FFT, IFFT and Hilbert transformations +*/ +/* + PIP - Platform Independent Primitives + Private header for fftw3 + Copyright (C) 2017 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 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#ifndef PIFFT_P_H +#define PIFFT_P_H + +#include "pivector.h" +#include "pimutex.h" +#include +#ifdef PIP_FFTW +# include "fftw3.h" +#else +# define FFTW_FORWARD 0 +# define FFTW_BACKWARD 0 +# define FFTW_ESTIMATE 0 +# define FFTW_MEASURE 0 +#endif + + +static PIMutex __pip_fft_plan_mutex; + + +template +class PIFFTW_Private +{ +public: + explicit PIFFTW_Private() {plan = 0; newPlan(plan);} + ~PIFFTW_Private() {deletePlan(plan);} + + const PIVector > & calcFFT(const PIVector > & in) { + result.resize(in.size()); + __pip_fft_plan_mutex.lock(); + createPlan_c_1d(plan, in.size_s(), in.data(), result.data(), FFTW_FORWARD, FFTW_ESTIMATE); + __pip_fft_plan_mutex.unlock(); + executePlan(plan); + destroyPlan(plan); + return result; + } + const PIVector > & calcFFT(const PIVector & in) { + result.resize(in.size()); + __pip_fft_plan_mutex.lock(); + createPlan_r2c_1d(plan, in.size_s(), in.data(), result.data(), FFTW_ESTIMATE); + __pip_fft_plan_mutex.unlock(); + executePlan(plan); + destroyPlan(plan); + return result; + } + const PIVector > & calcFFTinverse(const PIVector > & in) { + result.resize(in.size()); + __pip_fft_plan_mutex.lock(); + createPlan_c_1d(plan, in.size_s(), in.data(), result.data(), FFTW_BACKWARD, FFTW_ESTIMATE); + __pip_fft_plan_mutex.unlock(); + executePlan(plan); + destroyPlan(plan); + return result; + } + + enum FFT_Operation {fo_real, fo_complex, fo_inverse}; + + void preparePlan(int size, int op) { + PIVector > in(size), out(size); + PIVector inr(size); + __pip_fft_plan_mutex.lock(); + switch ((FFT_Operation)op) { + case fo_real: + createPlan_r2c_1d(plan, in.size_s(), inr.data(), out.data(), FFTW_MEASURE); + break; + case fo_complex: + createPlan_c_1d(plan, in.size_s(), in.data(), out.data(), FFTW_FORWARD, FFTW_MEASURE); + break; + case fo_inverse: + createPlan_c_1d(plan, in.size_s(), in.data(), out.data(), FFTW_BACKWARD, FFTW_MEASURE); + break; + } + __pip_fft_plan_mutex.unlock(); + } + + inline void createPlan_c_1d(void * plan, int size, const void * in, void * out, int dir, int flags) {} + inline void createPlan_r2c_1d(void * plan, int size, const void * in, void * out, int flags) {} + inline void executePlan(void * plan) {} + inline void destroyPlan(void * plan) {} + inline void newPlan(void *& plan) {} + inline void deletePlan(void *& plan) {} + + PIVector > result; + void * plan; +}; + + +#ifdef PIP_FFTW + +template<> inline void PIFFTW_Private::createPlan_c_1d(void * plan, int size, const void * in, void * out, int dir, int flags) { + *(fftwf_plan*)plan = fftwf_plan_dft_1d(size, (fftwf_complex *)in, (fftwf_complex *)out, dir, flags);} +template<> inline void PIFFTW_Private::createPlan_r2c_1d(void * plan, int size, const void * in, void * out, int flags) { + *(fftwf_plan*)plan = fftwf_plan_dft_r2c_1d(size, (float *)in, (fftwf_complex *)out, flags);} +template<> inline void PIFFTW_Private::executePlan(void * plan) {fftwf_execute(*(fftwf_plan*)plan);} +template<> inline void PIFFTW_Private::destroyPlan(void * plan) {fftwf_destroy_plan(*(fftwf_plan*)plan);} +template<> inline void PIFFTW_Private::newPlan(void *& plan) {plan = new fftwf_plan();} +template<> inline void PIFFTW_Private::deletePlan(void *& plan) {if (plan) delete (fftwf_plan*)plan; plan = 0;} + + +template<> inline void PIFFTW_Private::createPlan_c_1d(void * plan, int size, const void * in, void * out, int dir, int flags) { + *(fftw_plan*)plan = fftw_plan_dft_1d(size, (fftw_complex *)in, (fftw_complex *)out, dir, flags);} +template<> inline void PIFFTW_Private::createPlan_r2c_1d(void * plan, int size, const void * in, void * out, int flags) { + *(fftw_plan*)plan = fftw_plan_dft_r2c_1d(size, (double *)in, (fftw_complex *)out, flags);} +template<> inline void PIFFTW_Private::executePlan(void * plan) {fftw_execute(*(fftw_plan*)plan);} +template<> inline void PIFFTW_Private::destroyPlan(void * plan) {fftw_destroy_plan(*(fftw_plan*)plan);} +template<> inline void PIFFTW_Private::newPlan(void *& plan) {plan = new fftw_plan();} +template<> inline void PIFFTW_Private::deletePlan(void *& plan) {if (plan) delete (fftw_plan*)plan; plan = 0;} + + +template<> inline void PIFFTW_Private::createPlan_c_1d(void * plan, int size, const void * in, void * out, int dir, int flags) { + *(fftwl_plan*)plan = fftwl_plan_dft_1d(size, (fftwl_complex *)in, (fftwl_complex *)out, dir, flags);} +template<> inline void PIFFTW_Private::createPlan_r2c_1d(void * plan, int size, const void * in, void * out, int flags) { + *(fftwl_plan*)plan = fftwl_plan_dft_r2c_1d(size, (ldouble *)in, (fftwl_complex *)out, flags);} +template<> inline void PIFFTW_Private::executePlan(void * plan) {fftwl_execute(*(fftwl_plan*)plan);} +template<> inline void PIFFTW_Private::destroyPlan(void * plan) {fftwl_destroy_plan(*(fftwl_plan*)plan);} +template<> inline void PIFFTW_Private::newPlan(void *& plan) {plan = new fftwl_plan();} +template<> inline void PIFFTW_Private::deletePlan(void *& plan) {if (plan) delete (fftwl_plan*)plan; plan = 0;} + +#endif // PIP_FFTW + + +#endif // PIFFT_H