390 lines
17 KiB
C++
390 lines
17 KiB
C++
//! \~\file pifft.h
|
||
//! \~\ingroup Math
|
||
//! \~\ingroup FFTW
|
||
//! \~\brief
|
||
//! \~english FFT, IFFT and Hilbert transformations
|
||
//! \~russian БПФ, ОБПФ и преобразования Гильберта
|
||
/*
|
||
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"
|
||
|
||
//! \~\ingroup Math
|
||
//! \~\ingroup FFTW
|
||
//! \~\brief
|
||
//! \~english Double-precision FFT helper that stores the last transform result.
|
||
//! \~russian Вспомогательный класс БПФ двойной точности, сохраняющий результат последнего преобразования.
|
||
class PIP_EXPORT PIFFT_double {
|
||
public:
|
||
//! \~english Constructs an empty transformer.
|
||
//! \~russian Создает пустой объект преобразования.
|
||
PIFFT_double();
|
||
|
||
//! \~english Calculates complex FFT for the input signal.
|
||
//! \~russian Вычисляет комплексное БПФ для входного сигнала.
|
||
PIVector<complexd> * calcFFT(const PIVector<complexd> & val);
|
||
//! \~english Calculates FFT for a real-valued input signal.
|
||
//! \~russian Вычисляет БПФ для вещественного входного сигнала.
|
||
PIVector<complexd> * calcFFT(const PIVector<double> & val);
|
||
//! \~english Calculates inverse FFT for a complex spectrum.
|
||
//! \~russian Вычисляет обратное БПФ для комплексного спектра.
|
||
PIVector<complexd> * calcFFTinverse(const PIVector<complexd> & val);
|
||
//! \~english Calculates the analytic signal using a Hilbert transform.
|
||
//! \~russian Вычисляет аналитический сигнал с помощью преобразования Гильберта.
|
||
PIVector<complexd> * calcHilbert(const PIVector<double> & val);
|
||
//! \~english Returns magnitudes of the last calculated result.
|
||
//! \~russian Возвращает модули последнего вычисленного результата.
|
||
PIVector<double> getAmplitude() const;
|
||
//! \~english Returns real parts of the last calculated result.
|
||
//! \~russian Возвращает действительные части последнего вычисленного результата.
|
||
PIVector<double> getReal() const;
|
||
//! \~english Returns imaginary parts of the last calculated result.
|
||
//! \~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);
|
||
};
|
||
|
||
//! \~\ingroup Math
|
||
//! \~\ingroup FFTW
|
||
//! \~\brief
|
||
//! \~english Single-precision FFT helper that stores the last transform result.
|
||
//! \~russian Вспомогательный класс БПФ одинарной точности, сохраняющий результат последнего преобразования.
|
||
class PIP_EXPORT PIFFT_float {
|
||
public:
|
||
//! \~english Constructs an empty transformer.
|
||
//! \~russian Создает пустой объект преобразования.
|
||
PIFFT_float();
|
||
|
||
//! \~english Calculates complex FFT for the input signal.
|
||
//! \~russian Вычисляет комплексное БПФ для входного сигнала.
|
||
PIVector<complexf> * calcFFT(const PIVector<complexf> & val);
|
||
//! \~english Calculates FFT for a real-valued input signal.
|
||
//! \~russian Вычисляет БПФ для вещественного входного сигнала.
|
||
PIVector<complexf> * calcFFT(const PIVector<float> & val);
|
||
//! \~english Calculates inverse FFT for a complex spectrum.
|
||
//! \~russian Вычисляет обратное БПФ для комплексного спектра.
|
||
PIVector<complexf> * calcFFTinverse(const PIVector<complexf> & val);
|
||
//! \~english Calculates the analytic signal using a Hilbert transform.
|
||
//! \~russian Вычисляет аналитический сигнал с помощью преобразования Гильберта.
|
||
PIVector<complexf> * calcHilbert(const PIVector<float> & val);
|
||
//! \~english Returns magnitudes of the last calculated result.
|
||
//! \~russian Возвращает модули последнего вычисленного результата.
|
||
PIVector<float> getAmplitude() const;
|
||
//! \~english Returns real parts of the last calculated result.
|
||
//! \~russian Возвращает действительные части последнего вычисленного результата.
|
||
PIVector<float> getReal() const;
|
||
//! \~english Returns imaginary parts of the last calculated result.
|
||
//! \~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);
|
||
};
|
||
|
||
//! \~english Default FFT helper alias based on double precision.
|
||
//! \~russian Псевдоним стандартного FFT-помощника на базе двойной точности.
|
||
typedef PIFFT_double PIFFT;
|
||
//! \~english Double-precision FFT helper alias.
|
||
//! \~russian Псевдоним FFT-помощника двойной точности.
|
||
typedef PIFFT_double PIFFTd;
|
||
//! \~english Single-precision FFT helper alias.
|
||
//! \~russian Псевдоним FFT-помощника одинарной точности.
|
||
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)
|
||
|
||
//! \~\ingroup Math
|
||
//! \~\ingroup FFTW
|
||
//! \~\brief
|
||
//! \~english Thin wrapper over libfftw3 plans for a selected scalar type.
|
||
//! \~russian Тонкая обертка над планами libfftw3 для выбранного скалярного типа.
|
||
template<typename T>
|
||
class PIFFTW {
|
||
public:
|
||
//! \~english Constructs the backend object for the selected precision.
|
||
//! \~russian Создает внутренний объект для выбранной точности.
|
||
explicit PIFFTW() {
|
||
p = 0;
|
||
newP(p);
|
||
}
|
||
//! \~english Destroys the backend object and its cached plans.
|
||
//! \~russian Уничтожает внутренний объект и его кэшированные планы.
|
||
~PIFFTW() { deleteP(p); }
|
||
|
||
//! \~english Calculates complex FFT for the input signal.
|
||
//! \~russian Вычисляет комплексное БПФ для входного сигнала.
|
||
inline const PIVector<complex<T>> & calcFFT(const PIVector<complex<T>> & in) { return PIVector<complex<T>>().resize(in.size()); }
|
||
//! \~english Calculates FFT for a real-valued input signal.
|
||
//! \~russian Вычисляет БПФ для вещественного входного сигнала.
|
||
inline const PIVector<complex<T>> & calcFFT(const PIVector<T> & in) { return PIVector<complex<T>>().resize(in.size()); }
|
||
//! \~english Calculates inverse FFT for a complex spectrum.
|
||
//! \~russian Вычисляет обратное БПФ для комплексного спектра.
|
||
inline const PIVector<complex<T>> & calcFFTinverse(const PIVector<complex<T>> & in) { return PIVector<complex<T>>().resize(in.size()); }
|
||
|
||
//! \~english Operation kind for precomputed FFTW plans.
|
||
//! \~russian Тип операции для заранее подготовленных планов FFTW.
|
||
enum FFT_Operation {
|
||
foReal /** \~english Forward transform for real input \~russian Прямое преобразование для вещественного входа */,
|
||
foComplex /** \~english Forward transform for complex input \~russian Прямое преобразование для комплексного входа */,
|
||
foInverse /** \~english Inverse complex transform \~russian Обратное комплексное преобразование */
|
||
};
|
||
|
||
//! \~english Prepares and caches a plan for the specified signal size and operation.
|
||
//! \~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;
|
||
}
|
||
|
||
//! \~english Single-precision FFTW wrapper.
|
||
//! \~russian Обертка FFTW одинарной точности.
|
||
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;
|
||
}
|
||
|
||
//! \~english Double-precision FFTW wrapper.
|
||
//! \~russian Обертка FFTW двойной точности.
|
||
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;
|
||
}
|
||
|
||
//! \~english Extended-precision FFTW wrapper.
|
||
//! \~russian Обертка FFTW расширенной точности.
|
||
typedef PIFFTW<ldouble> PIFFTWld;
|
||
|
||
# endif
|
||
|
||
#endif // MICRO_PIP
|
||
|
||
#endif // PIFFT_H
|