diff options
Diffstat (limited to 'unsupported/Eigen/FFT')
| -rw-r--r-- | unsupported/Eigen/FFT | 101 |
1 files changed, 65 insertions, 36 deletions
diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT index 0cc10bafb..87c31749b 100644 --- a/unsupported/Eigen/FFT +++ b/unsupported/Eigen/FFT @@ -38,16 +38,16 @@ * #include <unsupported/Eigen/FFT> * \endcode * - * This module provides Fast Fourier transformation, either using a built-in implementation - * or as a frontend to various FFT libraries. + * This module provides Fast Fourier transformation, with a configurable backend + * implementation. * - * The build-in implementation is based on kissfft. It is a small, free, and + * The default implementation is based on kissfft. It is a small, free, and * reasonably efficient default. * - * There are currently two frontends: + * There are currently two implementation backend: * * - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size. - * - MLK (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form. + * - MKL (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form. * * \section FFTDesign Design * @@ -152,12 +152,20 @@ class FFT m_impl.fwd(dst,src,nfft); } + /* + inline + void fwd2(Complex * dst, const Complex * src, int n0,int n1) + { + m_impl.fwd2(dst,src,n0,n1); + } + */ + template <typename _Input> inline void fwd( std::vector<Complex> & 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); // half the bins + Nyquist bin else dst.resize(src.size()); fwd(&dst[0],&src[0],static_cast<int>(src.size())); @@ -179,28 +187,34 @@ class FFT dst.derived().resize( (src.size()>>1)+1); else dst.derived().resize(src.size()); - fwd( &dst[0],&src[0],src.size() ); + + if (src.stride() != 1) { + Matrix<typename InputDerived::Scalar,1,Dynamic> tmp = src; + fwd( &dst[0],&tmp[0],src.size() ); + }else{ + fwd( &dst[0],&src[0],src.size() ); + } } inline void inv( Complex * dst, const Complex * src, int nfft) { - m_impl.inv( dst,src,nfft ); - if ( HasFlag( Unscaled ) == false) - scale(dst,_Scalar(1./nfft),nfft); + m_impl.inv( dst,src,nfft ); + if ( HasFlag( Unscaled ) == false) + scale(dst,Scalar(1./nfft_,nfft); // scale the time series } inline void inv( Scalar * dst, const Complex * src, int nfft) { - m_impl.inv( dst,src,nfft ); - if ( HasFlag( Unscaled ) == false) - scale(dst,1./nfft,nfft); + m_impl.inv( dst,src,nfft ); + if ( HasFlag( Unscaled ) == false) + scale(dst,Scalar(1./nfft),nfft); // scale the time series } template<typename OutputDerived, typename ComplexDerived> inline - void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src) + void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src, int nfft=-1) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived) EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived) @@ -210,44 +224,59 @@ class FFT EIGEN_STATIC_ASSERT(int(OutputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit, THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES) - int nfft = src.size(); - int nout = HasFlag(HalfSpectrum) ? ((nfft>>1)+1) : nfft; - dst.derived().resize( nout ); - inv( &dst[0],&src[0], nfft); + if (nfft<1) { + nfft = ( NumTraits<typename OutputDerived::Scalar>::IsComplex == 0 && HasFlag(HalfSpectrum) ) ? 2*(src.size()-1) : src.size(); + } + dst.derived().resize( nfft ); + if (src.stride() != 1) { + // if the vector is strided, then we need to copy it to a packed temporary + Matrix<typename ComplexDerived::Scalar,1,Dynamic> tmp = src; + inv( &dst[0],&tmp[0], nfft); + }else{ + inv( &dst[0],&src[0], nfft); + } } template <typename _Output> inline - void inv( std::vector<_Output> & dst, const std::vector<Complex> & src) + void inv( std::vector<_Output> & dst, const std::vector<Complex> & src,int nfft=-1) { - if ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) - dst.resize( 2*(src.size()-1) ); - else - dst.resize( src.size() ); - inv( &dst[0],&src[0],static_cast<int>(dst.size()) ); + if (nfft<1) + nfft = ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) ? 2*(src.size()-1) : src.size(); + dst.resize( nfft ); + inv( &dst[0],&src[0],nfft); } + + /* // TODO: multi-dimensional FFTs - - // TODO: handle Eigen MatrixBase - // ---> i added fwd and inv specializations above + unit test, is this enough? (bjacob) + inline + void inv2(Complex * dst, const Complex * src, int n0,int n1) + { + m_impl.inv2(dst,src,n0,n1); + if ( HasFlag( Unscaled ) == false) + scale(dst,1./(n0*n1),n0*n1); + } + */ inline impl_type & impl() {return m_impl;} private: - template <typename _It,typename _Val> - inline void scale(_It x,_Val s,int nx) - { - for (int k=0;k<nx;++k) - *x++ *= _Scalar(s); - } - - template <typename _Val> - inline void scale(std::complex<_Val>* x,_Val s,int nx) + template <typename T_Data> + inline + void scale(T_Data * x,Scalar s,int nx) { +#if 1 for (int k=0;k<nx;++k) *x++ *= s; +#else + if ( ((ptrdiff_t)x) & 15 ) + Matrix<T_Data, Dynamic, 1>::Map(x,nx) *= s; + else + Matrix<T_Data, Dynamic, 1>::MapAligned(x,nx) *= s; + //Matrix<T_Data, Dynamic, Dynamic>::Map(x,nx) * s; +#endif } inline |