From cd7912313dc2477283de767029462d7d0e6ee8ab Mon Sep 17 00:00:00 2001 From: Mark Borgerding Date: Fri, 22 Jan 2010 00:35:03 -0500 Subject: [PATCH] changed FFT function vector and Matrix args to pointer as Benoit suggested implemented 2D Complex FFT for FFTW impl --- unsupported/Eigen/FFT | 45 +++++++++++------ unsupported/Eigen/src/FFT/ei_fftw_impl.h | 61 +++++++++++++++++++++- unsupported/test/FFT.cpp | 18 +++---- unsupported/test/FFTW.cpp | 64 +++++++++++++++++++++--- 4 files changed, 154 insertions(+), 34 deletions(-) diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT index e0841a4e3..caaf79714 100644 --- a/unsupported/Eigen/FFT +++ b/unsupported/Eigen/FFT @@ -152,20 +152,26 @@ class FFT m_impl.fwd(dst,src,nfft); } + inline + void fwd2(Complex * dst, const Complex * src, int nrows,int ncols) + { + m_impl.fwd2(dst,src,nrows,ncols); + } + template inline - void fwd( std::vector & dst, const std::vector<_Input> & src) + void fwd( std::vector * dst, const std::vector<_Input> & src) { if ( NumTraits<_Input>::IsComplex == 0 && HasFlag(HalfSpectrum) ) - dst.resize( (src.size()>>1)+1); + dst->resize( (src.size()>>1)+1); else - dst.resize(src.size()); - fwd(&dst[0],&src[0],static_cast(src.size())); + dst->resize(src.size()); + fwd(&(*dst)[0],&src[0],static_cast(src.size())); } template inline - void fwd( MatrixBase & dst, const MatrixBase & src) + void fwd( MatrixBase * dst, const MatrixBase & src) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived) EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived) @@ -176,10 +182,10 @@ class FFT THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES) if ( NumTraits< typename InputDerived::Scalar >::IsComplex == 0 && HasFlag(HalfSpectrum) ) - dst.derived().resize( (src.size()>>1)+1); + dst->derived().resize( (src.size()>>1)+1); else - dst.derived().resize(src.size()); - fwd( &dst[0],&src[0],src.size() ); + dst->derived().resize(src.size()); + fwd( &(*dst)[0],&src[0],src.size() ); } inline @@ -200,7 +206,7 @@ class FFT template inline - void inv( MatrixBase & dst, const MatrixBase & src) + void inv( MatrixBase * dst, const MatrixBase & src) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived) EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived) @@ -212,19 +218,28 @@ class FFT int nfft = src.size(); int nout = HasFlag(HalfSpectrum) ? ((nfft>>1)+1) : nfft; - dst.derived().resize( nout ); - inv( &dst[0],&src[0], nfft); + dst->derived().resize( nout ); + inv( &(*dst)[0],&src[0], nfft); } template inline - void inv( std::vector<_Output> & dst, const std::vector & src) + void inv( std::vector<_Output> * dst, const std::vector & src) { if ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) - dst.resize( 2*(src.size()-1) ); + dst->resize( 2*(src.size()-1) ); else - dst.resize( src.size() ); - inv( &dst[0],&src[0],static_cast(dst.size()) ); + dst->resize( src.size() ); + inv( &(*dst)[0],&src[0],static_cast(dst->size()) ); + } + + + inline + void inv2(Complex * dst, const Complex * src, int nrows,int ncols) + { + m_impl.inv2(dst,src,nrows,ncols); + if ( HasFlag( Unscaled ) == false) + scale(dst,1./(nrows*ncols),nrows*ncols); } // TODO: multi-dimensional FFTs diff --git a/unsupported/Eigen/src/FFT/ei_fftw_impl.h b/unsupported/Eigen/src/FFT/ei_fftw_impl.h index a66b7398c..411ff7425 100644 --- a/unsupported/Eigen/src/FFT/ei_fftw_impl.h +++ b/unsupported/Eigen/src/FFT/ei_fftw_impl.h @@ -90,6 +90,18 @@ m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE); fftwf_execute_dft_c2r( m_plan, src,dst); } + + inline + void fwd2( complex_type * dst,complex_type * src,int nrows,int ncols) { + if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(ncols,nrows,src,dst,FFTW_FORWARD,FFTW_ESTIMATE); + fftwf_execute_dft( m_plan, src,dst); + } + inline + void inv2( complex_type * dst,complex_type * src,int nrows,int ncols) { + if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(ncols,nrows,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE); + fftwf_execute_dft( m_plan, src,dst); + } + }; template <> struct ei_fftw_plan @@ -121,6 +133,16 @@ m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE); fftw_execute_dft_c2r( m_plan, src,dst); } + inline + void fwd2( complex_type * dst,complex_type * src,int nrows,int ncols) { + if (m_plan==NULL) m_plan = fftw_plan_dft_2d(ncols,nrows,src,dst,FFTW_FORWARD,FFTW_ESTIMATE); + fftw_execute_dft( m_plan, src,dst); + } + inline + void inv2( complex_type * dst,complex_type * src,int nrows,int ncols) { + if (m_plan==NULL) m_plan = fftw_plan_dft_2d(ncols,nrows,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE); + fftw_execute_dft( m_plan, src,dst); + } }; template <> struct ei_fftw_plan @@ -152,6 +174,16 @@ m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE); fftwl_execute_dft_c2r( m_plan, src,dst); } + inline + void fwd2( complex_type * dst,complex_type * src,int nrows,int ncols) { + if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(ncols,nrows,src,dst,FFTW_FORWARD,FFTW_ESTIMATE); + fftwl_execute_dft( m_plan, src,dst); + } + inline + void inv2( complex_type * dst,complex_type * src,int nrows,int ncols) { + if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(ncols,nrows,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE); + fftwl_execute_dft( m_plan, src,dst); + } }; template @@ -180,6 +212,13 @@ get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src) ,nfft); } + // 2-d complex-to-complex + inline + void fwd2(Complex * dst, const Complex * src, int nrows,int ncols) + { + get_plan(nrows,ncols,false,dst,src).fwd2(ei_fftw_cast(dst), ei_fftw_cast(src) ,nrows,ncols); + } + // inverse complex-to-complex inline void inv(Complex * dst,const Complex *src,int nfft) @@ -194,9 +233,18 @@ get_plan(nfft,true,dst,src).inv(ei_fftw_cast(dst), ei_fftw_cast(src),nfft ); } + // 2-d complex-to-complex + inline + void inv2(Complex * dst, const Complex * src, int nrows,int ncols) + { + get_plan(nrows,ncols,true,dst,src).inv2(ei_fftw_cast(dst), ei_fftw_cast(src) ,nrows,ncols); + } + + protected: typedef ei_fftw_plan PlanData; - typedef std::map PlanMap; + + typedef std::map PlanMap; PlanMap m_plans; @@ -205,7 +253,16 @@ { bool inplace = (dst==src); bool aligned = ( (reinterpret_cast(src)&15) | (reinterpret_cast(dst)&15) ) == 0; - int key = (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned; + int64_t key = ( (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1; + return m_plans[key]; + } + + inline + PlanData & get_plan(int nrows,int ncols,bool inverse,void * dst,const void * src) + { + bool inplace = (dst==src); + bool aligned = ( (reinterpret_cast(src)&15) | (reinterpret_cast(dst)&15) ) == 0; + int64_t key = ( ( (((int64_t)ncols) << 30)|(nrows<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1 ) + 1; return m_plans[key]; } }; diff --git a/unsupported/test/FFT.cpp b/unsupported/test/FFT.cpp index 056be2ef3..a2f1d9201 100644 --- a/unsupported/test/FFT.cpp +++ b/unsupported/test/FFT.cpp @@ -106,29 +106,29 @@ void test_scalar_generic(int nfft) // make sure it DOESN'T give the right full spectrum answer // if we've asked for half-spectrum fft.SetFlag(fft.HalfSpectrum ); - fft.fwd( outbuf,inbuf); + fft.fwd( &outbuf,inbuf); VERIFY(outbuf.size() == (size_t)( (nfft>>1)+1) ); VERIFY( fft_rmse(outbuf,inbuf) < test_precision() );// gross check fft.ClearFlag(fft.HalfSpectrum ); - fft.fwd( outbuf,inbuf); + fft.fwd( &outbuf,inbuf); VERIFY( fft_rmse(outbuf,inbuf) < test_precision() );// gross check ScalarVector buf3; - fft.inv( buf3 , outbuf); + fft.inv( &buf3 , outbuf); VERIFY( dif_rmse(inbuf,buf3) < test_precision() );// gross check // verify that the Unscaled flag takes effect ComplexVector buf4; fft.SetFlag(fft.Unscaled); - fft.inv( buf4 , outbuf); + fft.inv( &buf4 , outbuf); for (int k=0;k() );// gross check // verify that ClearFlag works fft.ClearFlag(fft.Unscaled); - fft.inv( buf3 , outbuf); + fft.inv( &buf3 , outbuf); VERIFY( dif_rmse(inbuf,buf3) < test_precision() );// gross check } @@ -152,25 +152,25 @@ void test_complex_generic(int nfft) ComplexVector buf3; for (int k=0;k() );// gross check - fft.inv( buf3 , outbuf); + fft.inv( &buf3 , outbuf); VERIFY( dif_rmse(inbuf,buf3) < test_precision() );// gross check // verify that the Unscaled flag takes effect ComplexVector buf4; fft.SetFlag(fft.Unscaled); - fft.inv( buf4 , outbuf); + fft.inv( &buf4 , outbuf); for (int k=0;k() );// gross check // verify that ClearFlag works fft.ClearFlag(fft.Unscaled); - fft.inv( buf3 , outbuf); + fft.inv( &buf3 , outbuf); VERIFY( dif_rmse(inbuf,buf3) < test_precision() );// gross check } diff --git a/unsupported/test/FFTW.cpp b/unsupported/test/FFTW.cpp index c182cab9d..df38efe64 100644 --- a/unsupported/test/FFTW.cpp +++ b/unsupported/test/FFTW.cpp @@ -26,7 +26,11 @@ #include #include +template +std::complex RandomCpx() { return std::complex( (T)(rand()/(T)RAND_MAX - .5), (T)(rand()/(T)RAND_MAX - .5) ); } + using namespace std; +using namespace Eigen; float norm(float x) {return x*x;} double norm(double x) {return x*x;} @@ -87,11 +91,11 @@ void test_scalar(int nfft) vector outbuf; for (int k=0;k() );// gross check vector buf3; - fft.inv( buf3 , outbuf); + fft.inv( &buf3 , outbuf); VERIFY( dif_rmse(inbuf,buf3) < test_precision() );// gross check } @@ -106,19 +110,65 @@ void test_complex(int nfft) vector outbuf; vector buf3; for (int k=0;k(); + fft.fwd( &outbuf , inbuf); VERIFY( fft_rmse(outbuf,inbuf) < test_precision() );// gross check - fft.inv( buf3 , outbuf); + fft.inv( &buf3 , outbuf); VERIFY( dif_rmse(inbuf,buf3) < test_precision() );// gross check } -void test_FFTW() +template +void test_complex2d() { + typedef typename Eigen::FFT::Complex Complex; + FFT fft; + + Eigen::Matrix src; + Eigen::Matrix dst; + Eigen::Matrix src2; + Eigen::Matrix dst2; + + //src = Eigen::Matrix::Random(); + src = Eigen::Matrix::Identity(); + + for (int k=0;k tmpIn = src.col(k); + Eigen::Matrix tmpOut; + fft.fwd( &tmpOut,tmpIn ); + dst2.col(k) = tmpOut; + } + //cout << "dst2: " << dst2 << "\n\n"; + + for (int k=0;k tmpIn = dst2.row(k); + Eigen::Matrix tmpOut; + fft.fwd( &tmpOut, tmpIn); + dst2.row(k) = tmpOut; + } + +/* +*/ + fft.fwd2(dst.data(),src.data(),nrows,ncols); + fft.inv2(src2.data(),dst.data(),nrows,ncols); + /* + cout << "src: " << src << "\n\n"; + cout << "dst: " << dst << "\n\n"; + cout << "src2: " << src2 << "\n\n"; + cout << "dst2: " << dst2 << "\n\n"; + */ + VERIFY( (src-src2).norm() < test_precision() ); + VERIFY( (dst-dst2).norm() < test_precision() ); +} + +void test_FFTW() +{ + CALL_SUBTEST( ( test_complex2d () ) ); + CALL_SUBTEST( ( test_complex2d () ) ); + //CALL_SUBTEST( ( test_complex2d () ) ); CALL_SUBTEST( test_complex(32) ); CALL_SUBTEST( test_complex(32) ); CALL_SUBTEST( test_complex(32) ); CALL_SUBTEST( test_complex(256) ); CALL_SUBTEST( test_complex(256) ); CALL_SUBTEST( test_complex(256) ); CALL_SUBTEST( test_complex(3*8) ); CALL_SUBTEST( test_complex(3*8) ); CALL_SUBTEST( test_complex(3*8) ); @@ -127,8 +177,6 @@ void test_FFTW() CALL_SUBTEST( test_complex(2*3*4*5) ); CALL_SUBTEST( test_complex(2*3*4*5) ); CALL_SUBTEST( test_complex(2*3*4*5) ); CALL_SUBTEST( test_complex(2*3*4*5*7) ); CALL_SUBTEST( test_complex(2*3*4*5*7) ); CALL_SUBTEST( test_complex(2*3*4*5*7) ); - - CALL_SUBTEST( test_scalar(32) ); CALL_SUBTEST( test_scalar(32) ); CALL_SUBTEST( test_scalar(32) ); CALL_SUBTEST( test_scalar(45) ); CALL_SUBTEST( test_scalar(45) ); CALL_SUBTEST( test_scalar(45) ); CALL_SUBTEST( test_scalar(50) ); CALL_SUBTEST( test_scalar(50) ); CALL_SUBTEST( test_scalar(50) );