diff options
author | Mark Borgerding <mark@borgerding.net> | 2009-10-21 20:53:05 -0400 |
---|---|---|
committer | Mark Borgerding <mark@borgerding.net> | 2009-10-21 20:53:05 -0400 |
commit | e3d08443dc272f740447de0147efc69cf7de1c93 (patch) | |
tree | 1a33e447caafe134731d3388f789d6ed99d50acd | |
parent | 78a53574b7bf386e4786676554b638d917cf8b62 (diff) |
inlining,all namespace declaration moved to FFT, removed preprocessor definitions,
-rw-r--r-- | unsupported/Eigen/FFT | 41 | ||||
-rw-r--r-- | unsupported/Eigen/src/FFT/ei_fftw_impl.h | 38 | ||||
-rw-r--r-- | unsupported/Eigen/src/FFT/ei_kissfft_impl.h | 28 | ||||
-rw-r--r-- | unsupported/test/CMakeLists.txt | 2 |
4 files changed, 73 insertions, 36 deletions
diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT index 97ec8e49b..36afdde8d 100644 --- a/unsupported/Eigen/FFT +++ b/unsupported/Eigen/FFT @@ -25,29 +25,39 @@ #ifndef EIGEN_FFT_H #define EIGEN_FFT_H -// ei_kissfft_impl: small, free, reasonably efficient default, derived from kissfft -#include "src/FFT/ei_kissfft_impl.h" -#define DEFAULT_FFT_IMPL ei_kissfft_impl +#include <complex> +#include <vector> +#include <map> +#ifdef EIGEN_FFTW_DEFAULT // FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size -#ifdef FFTW_ESTIMATE // definition of FFTW_ESTIMATE indicates the caller has included fftw3.h, we can use FFTW routines -#include "src/FFT/ei_fftw_impl.h" -#undef DEFAULT_FFT_IMPL -#define DEFAULT_FFT_IMPL ei_fftw_impl -#endif - -// intel Math Kernel Library: fastest, commercial -- incompatible with Eigen in GPL form -#ifdef _MKL_DFTI_H_ // mkl_dfti.h has been included, we can use MKL FFT routines +# include <fftw3.h> + namespace Eigen { +# include "src/FFT/ei_fftw_impl.h" + //template <typename T> typedef struct ei_fftw_impl default_fft_impl; this does not work + template <typename T> struct default_fft_impl : public ei_fftw_impl<T> {}; + } +#elif defined EIGEN_MKL_DEFAULT // TODO -// #include "src/FFT/ei_imkl_impl.h" -// #define DEFAULT_FFT_IMPL ei_imkl_impl +// intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form + namespace Eigen { +# include "src/FFT/ei_imklfft_impl.h" + template <typename T> struct default_fft_impl : public ei_imklfft_impl {}; + } +#else +// ei_kissfft_impl: small, free, reasonably efficient default, derived from kissfft +// + namespace Eigen { +# include "src/FFT/ei_kissfft_impl.h" + template <typename T> + struct default_fft_impl : public ei_kissfft_impl<T> {}; + } #endif namespace Eigen { template <typename _Scalar, - typename _Impl=DEFAULT_FFT_IMPL<_Scalar> - > + typename _Impl=default_fft_impl<_Scalar> > class FFT { public: @@ -120,7 +130,6 @@ class FFT private: impl_type m_impl; }; -#undef DEFAULT_FFT_IMPL } #endif /* vim: set filetype=cpp et sw=2 ts=2 ai: */ diff --git a/unsupported/Eigen/src/FFT/ei_fftw_impl.h b/unsupported/Eigen/src/FFT/ei_fftw_impl.h index d592bbb20..e1f67f334 100644 --- a/unsupported/Eigen/src/FFT/ei_fftw_impl.h +++ b/unsupported/Eigen/src/FFT/ei_fftw_impl.h @@ -22,7 +22,8 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -namespace Eigen { + + // FFTW uses non-const arguments // so we must use ugly const_cast calls for all the args it uses // @@ -32,21 +33,25 @@ namespace Eigen { // 2. fftw_complex is compatible with std::complex // This assumes std::complex<T> layout is array of size 2 with real,imag template <typename T> + inline T * ei_fftw_cast(const T* p) { return const_cast<T*>( p); } + inline fftw_complex * ei_fftw_cast( const std::complex<double> * p) - { + { return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) ); } + inline fftwf_complex * ei_fftw_cast( const std::complex<float> * p) { return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) ); } + inline fftwl_complex * ei_fftw_cast( const std::complex<long double> * p) { return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) ); @@ -64,18 +69,22 @@ namespace Eigen { ei_fftw_plan() :m_plan(NULL) {} ~ei_fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);} + inline void fwd(complex_type * dst,complex_type * src,int nfft) { if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE); fftwf_execute_dft( m_plan, src,dst); } + inline void inv(complex_type * dst,complex_type * src,int nfft) { if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE); fftwf_execute_dft( m_plan, src,dst); } + inline void fwd(complex_type * dst,scalar_type * src,int nfft) { if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE); fftwf_execute_dft_r2c( m_plan,src,dst); } + inline void inv(scalar_type * dst,complex_type * src,int nfft) { if (m_plan==NULL) m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE); @@ -91,18 +100,22 @@ namespace Eigen { ei_fftw_plan() :m_plan(NULL) {} ~ei_fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);} + inline void fwd(complex_type * dst,complex_type * src,int nfft) { if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE); fftw_execute_dft( m_plan, src,dst); } + inline void inv(complex_type * dst,complex_type * src,int nfft) { if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE); fftw_execute_dft( m_plan, src,dst); } + inline void fwd(complex_type * dst,scalar_type * src,int nfft) { if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE); fftw_execute_dft_r2c( m_plan,src,dst); } + inline void inv(scalar_type * dst,complex_type * src,int nfft) { if (m_plan==NULL) m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE); @@ -118,18 +131,22 @@ namespace Eigen { ei_fftw_plan() :m_plan(NULL) {} ~ei_fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);} + inline void fwd(complex_type * dst,complex_type * src,int nfft) { if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE); fftwl_execute_dft( m_plan, src,dst); } + inline void inv(complex_type * dst,complex_type * src,int nfft) { if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE); fftwl_execute_dft( m_plan, src,dst); } + inline void fwd(complex_type * dst,scalar_type * src,int nfft) { if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE); fftwl_execute_dft_r2c( m_plan,src,dst); } + inline void inv(scalar_type * dst,complex_type * src,int nfft) { if (m_plan==NULL) m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE); @@ -143,17 +160,20 @@ namespace Eigen { typedef _Scalar Scalar; typedef std::complex<Scalar> Complex; + inline void clear() { m_plans.clear(); } + inline void fwd( Complex * dst,const Complex *src,int nfft) { get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src),nfft ); } // real-to-complex forward FFT + inline void fwd( Complex * dst,const Scalar * src,int nfft) { get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src) ,nfft); @@ -163,30 +183,37 @@ namespace Eigen { } // inverse complex-to-complex + inline void inv(Complex * dst,const Complex *src,int nfft) { get_plan(nfft,true,dst,src).inv(ei_fftw_cast(dst), ei_fftw_cast(src),nfft ); + + //TODO move scaling to Eigen::FFT // scaling - Scalar s = 1./nfft; + Scalar s = Scalar(1.)/nfft; for (int k=0;k<nfft;++k) dst[k] *= s; } // half-complex to scalar + inline void inv( Scalar * dst,const Complex * src,int nfft) { get_plan(nfft,true,dst,src).inv(ei_fftw_cast(dst), ei_fftw_cast(src),nfft ); - Scalar s = 1./nfft; + + //TODO move scaling to Eigen::FFT + Scalar s = Scalar(1.)/nfft; for (int k=0;k<nfft;++k) dst[k] *= s; } - private: + protected: typedef ei_fftw_plan<Scalar> PlanData; typedef std::map<int,PlanData> PlanMap; PlanMap m_plans; + inline PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src) { bool inplace = (dst==src); @@ -195,4 +222,3 @@ namespace Eigen { return m_plans[key]; } }; -} diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h index a84ac68a0..c068d8765 100644 --- a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h +++ b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h @@ -22,11 +22,7 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -#include <complex> -#include <vector> -#include <map> -namespace Eigen { // This FFT implementation was derived from kissfft http:sourceforge.net/projects/kissfft // Copyright 2003-2009 Mark Borgerding @@ -51,13 +47,6 @@ namespace Eigen { m_twiddles[i] = exp( Complex(0,i*phinc) ); } - void conjugate() - { - m_inverse = !m_inverse; - for ( size_t i=0;i<m_twiddles.size() ;++i) - m_twiddles[i] = conj( m_twiddles[i] ); - } - void factorize(int nfft) { //start factoring out 4's, then 2's, then 3,5,7,9,... @@ -116,6 +105,7 @@ namespace Eigen { } } + inline void bfly2( Complex * Fout, const size_t fstride, int m) { for (int k=0;k<m;++k) { @@ -125,6 +115,7 @@ namespace Eigen { } } + inline void bfly4( Complex * Fout, const size_t fstride, const size_t m) { Complex scratch[6]; @@ -147,6 +138,7 @@ namespace Eigen { } } + inline void bfly3( Complex * Fout, const size_t fstride, const size_t m) { size_t k=m; @@ -175,6 +167,7 @@ namespace Eigen { }while(--k); } + inline void bfly5( Complex * Fout, const size_t fstride, const size_t m) { Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; @@ -241,6 +234,7 @@ namespace Eigen { } /* perform the butterfly for one stage of a mixed radix FFT */ + inline void bfly_generic( Complex * Fout, const size_t fstride, @@ -290,6 +284,7 @@ namespace Eigen { } template <typename _Src> + inline void fwd( Complex * dst,const _Src *src,int nfft) { get_plan(nfft,false).work(0, dst, src, 1,1); @@ -299,6 +294,7 @@ namespace Eigen { // perform two FFTs of src even and src odd // then twiddle to recombine them into the half-spectrum format // then fill in the conjugate symmetric half + inline void fwd( Complex * dst,const Scalar * src,int nfft) { if ( nfft&3 ) { @@ -334,6 +330,7 @@ namespace Eigen { } // inverse complex-to-complex + inline void inv(Complex * dst,const Complex *src,int nfft) { get_plan(nfft,true).work(0, dst, src, 1,1); @@ -341,6 +338,7 @@ namespace Eigen { } // half-complex to scalar + inline void inv( Scalar * dst,const Complex * src,int nfft) { if (nfft&3) { @@ -369,7 +367,7 @@ namespace Eigen { } } - private: + protected: typedef ei_kiss_cpx_fft<Scalar> PlanData; typedef std::map<int,PlanData> PlanMap; @@ -377,8 +375,10 @@ namespace Eigen { std::map<int, std::vector<Complex> > m_realTwiddles; std::vector<Complex> m_tmpBuf; + inline int PlanKey(int nfft,bool isinverse) const { return (nfft<<1) | isinverse; } + inline PlanData & get_plan(int nfft,bool inverse) { // TODO look for PlanKey(nfft, ! inverse) and conjugate the twiddles @@ -390,6 +390,7 @@ namespace Eigen { return pd; } + inline Complex * real_twiddles(int ncfft2) { std::vector<Complex> & twidref = m_realTwiddles[ncfft2];// creates new if not there @@ -403,10 +404,11 @@ namespace Eigen { return &twidref[0]; } + // TODO move scaling up into Eigen::FFT + inline void scale(Complex *dst,int n,Scalar s) { for (int k=0;k<n;++k) dst[k] *= s; } }; -} diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index f42077bdc..d182c9abf 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -23,6 +23,6 @@ ei_add_test(FFT) find_package(FFTW) if(FFTW_FOUND) - ei_add_test(FFTW " " "-lfftw3 -lfftw3f -lfftw3l" ) + ei_add_test(FFTW "-DEIGEN_FFTW_DEFAULT " "-lfftw3 -lfftw3f -lfftw3l" ) endif(FFTW_FOUND) |