/*! \file fftwaff.h * \brief use fftw together with aff containers (prototypes) * * ---------------------------------------------------------------------------- * * \author Thomas Forbriger * \date 11/07/2006 * * use fftw together with aff containers (prototypes) * * link with -lrfftw -lfftw -lm -laff * * ---- * libfourier 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 2 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, write to the Free Software * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA * ---- * * * Copyright (c) 2006 by Thomas Forbriger (BFO Schiltach) * * REVISIONS and CHANGES * - 11/07/2006 V1.0 Thomas Forbriger * - 12/09/2007 V1.1 first running version * - 07/10/2010 V1.2 * - migrate to FFTW3: * - use different fftw header file (fftw3.h) * - type of FFTW plan has changed * - 02/10/2012 V1.3 * - provide size calculation functions * * Migration to FFTW3: * - The location of the arrays in memory are part of the plan. * Consequently we have to either create a new plan for each transform or to * allocate the working array together with the plan and keep it available in * the background. The solution of choice is to keep arrays. We will * introduce a control parameter to the arguments of the constructor which * optionally allows to get rid of the arrays after each transformation. * - FFTW3 uses wisdom by default. * \code * void fftw_forget_wisdom(void); * \endcode * should be called by the destructor of the class. Alternatively * \code * void fftw_cleanup(void); * \endcode * could be called which apparently is even more rigorous. Plans, however, * must be destroyed prior to calling \c fftw_cleanup. * - The FFTW documentation states: * To the extent that this is true, if you have a variable * \code * complex *x, * \endcode * you can pass it directly to FFTW via * \code * reinterpret_cast(x). * \endcode * * ============================================================================ */ // include guard #ifndef TF_FFTWAFF_H_VERSION #define TF_FFTWAFF_H_VERSION \ "TF_FFTWAFF_H V1.3" #include #ifdef FFTWFALLBACK #include #else #include #endif #include namespace fourier { /*! All Fourier transform stuff is collected here. */ namespace fft { /*! A rigid class to do simple transforms using libdrfftw.a. * * uses real double arrays * * How to use this class: * * You may create one instance of this class and use it to transform * several signals in both directions in turn. The class itself takes care * of the transform size and creates a new plan if necessary. FFTs are * invoked by the bracket operators. You should use the class object like * a function. The scaling operators (taking the sampling interval as one * of their arguments) return a series or spectrum scaled appropriately to * match the values of samples from the corresponding Fourier integral * transform (usual convention with \f$dt\f$ and \f$df\f$ * integrals - not \f$d \omega\f$). * * \note * The appropriate number of samples for the time series obtained from a * given set of Fourier coefficients is non-unique, if only the Fourier * coefficients for positive frequencies are given, as is the case here. * A time series with \f$2n\f$ samples and a time series with \f$2n+1\f$ * samples both will results in a set of \f$n+1\f$ Fourier coefficients. * The method to determine the required number of time samples from the * imaginary part of the Nyquist Fourier coefficients only works if the * Nyquist coefficient (at least real part) is finite - which is not the * case for most of our signals. This way we will obtain \f$n+1\f$ Fourier * coefficients from \f$2n\f$ time series samples. When reconstructing the * time series, we will obtain \f$2n+1\f$ samples with a sampling interval * now being \f$2n/(2n+1)\f$ of the original interval. For this reason it * is appropriate to set the size of the expected time series explicitely * either by using the constructor * DRFFTWAFF::DRFFTWAFF(const unsigned int& n) or by * calling DRFFTWAFF::size(const unsigned int& s) prior to * DRFFTWAFF::operator()(const Tspectrum::Tcoc& s). * * \sa \ref page_fftw3 */ class DRFFTWAFF { public: typedef double Tsample; typedef std::complex Tcoeff; typedef aff::Series Tseries; typedef aff::Series Tspectrum; #ifdef FFTWFALLBACK DRFFTWAFF(const unsigned int& n): Msize(n), Mplan_forward(0), Mplan_backward(0) { } DRFFTWAFF(): Msize(0), Mplan_forward(0), Mplan_backward(0) { } #else DRFFTWAFF(const unsigned int& n, const bool& deletearrays=false): Msize(n), Mplan_forward(0), Mplan_backward(0), Mseriesarray(0), Mspectrumarray(0), Mdeletearrays(deletearrays) { } DRFFTWAFF(const bool& deletearrays=false): Msize(0), Mplan_forward(0), Mplan_backward(0), Mseriesarray(0), Mspectrumarray(0), Mdeletearrays(deletearrays) { } #endif ~DRFFTWAFF(); Tspectrum operator()(const Tseries::Tcoc& s, const bool& debug=false) const; Tseries operator()(const Tspectrum::Tcoc& s, const bool& debug=false) const; Tspectrum operator()(const Tseries::Tcoc& s, const Tsample& dt, const bool& debug=false) const; Tseries operator()(const Tspectrum::Tcoc& s, const Tsample& dt, const bool& debug=false) const; Tsample scale_series(const Tsample& dt) const; Tsample scale_spectrum(const Tsample& dt) const; unsigned int size() const { return(Msize); } void size(const unsigned int& s) const { this->set_size(s); } //! return number of coefficients for given number of samples inline static unsigned int spectrumsize(const unsigned int& n) { return(n/2+1); } //! return number of samples for given number of coefficients inline static unsigned int seriessize(const unsigned int& n) { return(n*2-1); } private: void create_plan_forward() const; void create_plan_backward() const; void delete_plans() const; #ifndef FFTWFALLBACK void create_arrays() const; void delete_arrays() const; unsigned int ssize() const { return(DRFFTWAFF::spectrumsize(this->size())); } #endif void set_size(const unsigned int& n) const; mutable unsigned int Msize; #ifdef FFTWFALLBACK mutable rfftw_plan Mplan_forward; mutable rfftw_plan Mplan_backward; #else mutable fftw_plan Mplan_forward; mutable fftw_plan Mplan_backward; mutable double *Mseriesarray; mutable fftw_complex *Mspectrumarray; mutable Tspectrum Mspectrum; mutable Tseries Mseries; bool Mdeletearrays; #endif }; // class DRFFTWAFF } // namespace ftt } // namespace fourier #endif // TF_FFTWAFF_H_VERSION (includeguard) /* ----- END OF fftwaff.h ----- */