diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-11-09 09:08:03 -0500 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-11-09 09:08:03 -0500 |
commit | 92749eed11d000300cfa54654f1043cd52399ed8 (patch) | |
tree | ba227522582b2f9f4280ed1404e74c654e21ccb3 /unsupported | |
parent | 4b366b07be4e409239c61158a23d93e8ebf3811b (diff) | |
parent | 670651e2e0932c5edfe2a2da4b9f3c42af3b7dec (diff) |
* merge
* remove a ctor in QuaternionBase as it gives a strange error with GCC 4.4.2.
Diffstat (limited to 'unsupported')
-rw-r--r-- | unsupported/Eigen/Complex | 75 | ||||
-rw-r--r-- | unsupported/Eigen/FFT | 115 | ||||
-rw-r--r-- | unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h | 2 | ||||
-rw-r--r-- | unsupported/Eigen/src/AutoDiff/AutoDiffVector.h | 128 | ||||
-rw-r--r-- | unsupported/Eigen/src/FFT/ei_fftw_impl.h | 17 | ||||
-rw-r--r-- | unsupported/Eigen/src/FFT/ei_kissfft_impl.h | 702 | ||||
-rw-r--r-- | unsupported/test/CMakeLists.txt | 1 | ||||
-rw-r--r-- | unsupported/test/Complex.cpp | 77 | ||||
-rw-r--r-- | unsupported/test/FFT.cpp | 35 |
9 files changed, 695 insertions, 457 deletions
diff --git a/unsupported/Eigen/Complex b/unsupported/Eigen/Complex index e1c41ab38..04228c95a 100644 --- a/unsupported/Eigen/Complex +++ b/unsupported/Eigen/Complex @@ -33,14 +33,26 @@ namespace Eigen { -template <typename _NativePtr,typename _PunnedPtr> +template <typename _NativeData,typename _PunnedData> struct castable_pointer { - castable_pointer(_NativePtr ptr) : _ptr(ptr) {} - operator _NativePtr () {return _ptr;} - operator _PunnedPtr () {return reinterpret_cast<_PunnedPtr>(_ptr);} + castable_pointer(_NativeData * ptr) : _ptr(ptr) { } + operator _NativeData * () {return _ptr;} + operator _PunnedData * () {return reinterpret_cast<_PunnedData*>(_ptr);} + operator const _NativeData * () const {return _ptr;} + operator const _PunnedData * () const {return reinterpret_cast<_PunnedData*>(_ptr);} private: - _NativePtr _ptr; + _NativeData * _ptr; +}; + +template <typename _NativeData,typename _PunnedData> +struct const_castable_pointer +{ + const_castable_pointer(_NativeData * ptr) : _ptr(ptr) { } + operator const _NativeData * () const {return _ptr;} + operator const _PunnedData * () const {return reinterpret_cast<_PunnedData*>(_ptr);} + private: + _NativeData * _ptr; }; template <typename T> @@ -50,7 +62,8 @@ struct Complex typedef T value_type; // constructors - Complex(const T& re = T(), const T& im = T()) : _re(re),_im(im) { } + Complex() {} + Complex(const T& re, const T& im = T()) : _re(re),_im(im) { } Complex(const Complex&other ): _re(other.real()) ,_im(other.imag()) {} template<class X> @@ -58,40 +71,63 @@ struct Complex template<class X> Complex(const std::complex<X>&other): _re(other.real()) ,_im(other.imag()) {} - // allow binary access to the object as a std::complex - typedef castable_pointer< Complex<T>*, StandardComplex* > pointer_type; - typedef castable_pointer< const Complex<T>*, const StandardComplex* > const_pointer_type; + typedef castable_pointer< Complex<T>, StandardComplex > pointer_type; + typedef const_castable_pointer< Complex<T>, StandardComplex > const_pointer_type; + + inline pointer_type operator & () {return pointer_type(this);} + + inline const_pointer_type operator & () const {return const_pointer_type(this);} + inline operator StandardComplex () const {return std_type();} + inline operator StandardComplex & () {return std_type();} - StandardComplex std_type() const {return StandardComplex(real(),imag());} + inline + const StandardComplex & std_type() const {return *reinterpret_cast<const StandardComplex*>(this);} + + inline StandardComplex & std_type() {return *reinterpret_cast<StandardComplex*>(this);} // every sort of accessor and mutator that has ever been in fashion. // For a brief history, search for "std::complex over-encapsulated" // http://www.open-std.org/jtc1/sc22/wg21/docs/lwg-defects.html#387 + inline const T & real() const {return _re;} + inline const T & imag() const {return _im;} + inline T & real() {return _re;} + inline T & imag() {return _im;} + inline T & real(const T & x) {return _re=x;} + inline T & imag(const T & x) {return _im=x;} + inline void set_real(const T & x) {_re = x;} + inline void set_imag(const T & x) {_im = x;} // *** complex member functions: *** + inline Complex<T>& operator= (const T& val) { _re=val;_im=0;return *this; } + inline Complex<T>& operator+= (const T& val) {_re+=val;return *this;} + inline Complex<T>& operator-= (const T& val) {_re-=val;return *this;} + inline Complex<T>& operator*= (const T& val) {_re*=val;_im*=val;return *this; } + inline Complex<T>& operator/= (const T& val) {_re/=val;_im/=val;return *this; } + inline Complex& operator= (const Complex& rhs) {_re=rhs._re;_im=rhs._im;return *this;} + inline Complex& operator= (const StandardComplex& rhs) {_re=rhs.real();_im=rhs.imag();return *this;} template<class X> Complex<T>& operator= (const Complex<X>& rhs) { _re=rhs._re;_im=rhs._im;return *this;} @@ -105,8 +141,7 @@ struct Complex T _im; }; -template <typename T> -T ei_to_std( const T & x) {return x;} +//template <typename T> T ei_to_std( const T & x) {return x;} template <typename T> std::complex<T> ei_to_std( const Complex<T> & x) {return x.std_type();} @@ -165,7 +200,7 @@ operator<< (std::basic_ostream<charT,traits>& ostr, const Complex<T>& rhs) template<class T> Complex<T> log (const Complex<T>&x){return log(ei_to_std(x));} template<class T> Complex<T> log10 (const Complex<T>&x){return log10(ei_to_std(x));} - template<class T> Complex<T> pow(const Complex<T>&x, int p) {return pow(ei_to_std(x),ei_to_std(p));} + template<class T> Complex<T> pow(const Complex<T>&x, int p) {return pow(ei_to_std(x),p);} template<class T> Complex<T> pow(const Complex<T>&x, const T&p) {return pow(ei_to_std(x),ei_to_std(p));} template<class T> Complex<T> pow(const Complex<T>&x, const Complex<T>&p) {return pow(ei_to_std(x),ei_to_std(p));} template<class T> Complex<T> pow(const T&x, const Complex<T>&p) {return pow(ei_to_std(x),ei_to_std(p));} @@ -175,8 +210,20 @@ operator<< (std::basic_ostream<charT,traits>& ostr, const Complex<T>& rhs) template<class T> Complex<T> sqrt (const Complex<T>&x){return sqrt(ei_to_std(x));} template<class T> Complex<T> tan (const Complex<T>&x){return tan(ei_to_std(x));} template<class T> Complex<T> tanh (const Complex<T>&x){return tanh(ei_to_std(x));} -} + template<typename _Real> struct NumTraits<Complex<_Real> > + { + typedef _Real Real; + typedef Complex<_Real> FloatingPoint; + enum { + IsComplex = 1, + HasFloatingPoint = NumTraits<Real>::HasFloatingPoint, + ReadCost = 2, + AddCost = 2 * NumTraits<Real>::AddCost, + MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost + }; + }; +} #endif /* vim: set filetype=cpp et sw=2 ts=2 ai: */ diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT index 36afdde8d..8f7a2fcae 100644 --- a/unsupported/Eigen/FFT +++ b/unsupported/Eigen/FFT @@ -28,6 +28,7 @@ #include <complex> #include <vector> #include <map> +#include <Eigen/Core> #ifdef EIGEN_FFTW_DEFAULT // FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size @@ -65,49 +66,87 @@ class FFT typedef typename impl_type::Scalar Scalar; typedef typename impl_type::Complex Complex; - FFT(const impl_type & impl=impl_type() ) :m_impl(impl) { } + enum Flag { + Default=0, // goof proof + Unscaled=1, + HalfSpectrum=2, + // SomeOtherSpeedOptimization=4 + Speedy=32767 + }; - template <typename _Input> - void fwd( Complex * dst, const _Input * src, int nfft) + FFT( const impl_type & impl=impl_type() , Flag flags=Default ) :m_impl(impl),m_flag(flags) { } + + inline + bool HasFlag(Flag f) const { return (m_flag & (int)f) == f;} + + inline + void SetFlag(Flag f) { m_flag |= (int)f;} + + inline + void ClearFlag(Flag f) { m_flag &= (~(int)f);} + + inline + void fwd( Complex * dst, const Scalar * src, int nfft) + { + m_impl.fwd(dst,src,nfft); + if ( HasFlag(HalfSpectrum) == false) + ReflectSpectrum(dst,nfft); + } + + inline + void fwd( Complex * dst, const Complex * src, int nfft) { m_impl.fwd(dst,src,nfft); } template <typename _Input> + inline void fwd( std::vector<Complex> & dst, const std::vector<_Input> & src) { - dst.resize( src.size() ); - fwd( &dst[0],&src[0],src.size() ); + if ( NumTraits<_Input>::IsComplex == 0 && HasFlag(HalfSpectrum) ) + dst.resize( (src.size()>>1)+1); + else + dst.resize(src.size()); + fwd(&dst[0],&src[0],src.size()); } template<typename InputDerived, typename ComplexDerived> + inline void fwd( MatrixBase<ComplexDerived> & dst, const MatrixBase<InputDerived> & src) { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived) - EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived) - EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,InputDerived) // size at compile-time - EIGEN_STATIC_ASSERT((ei_is_same_type<typename ComplexDerived::Scalar, Complex>::ret), - YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) - EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit, - THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES) - dst.derived().resize( src.size() ); - fwd( &dst[0],&src[0],src.size() ); + EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived) + EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived) + EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,InputDerived) // size at compile-time + EIGEN_STATIC_ASSERT((ei_is_same_type<typename ComplexDerived::Scalar, Complex>::ret), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit, + 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); + else + dst.derived().resize(src.size()); + fwd( &dst[0],&src[0],src.size() ); } - template <typename _Output> - void inv( _Output * dst, const Complex * src, int nfft) + inline + void inv( Complex * dst, const Complex * src, int nfft) { m_impl.inv( dst,src,nfft ); + if ( HasFlag( Unscaled ) == false) + scale(dst,1./nfft,nfft); } - template <typename _Output> - void inv( std::vector<_Output> & dst, const std::vector<Complex> & src) + inline + void inv( Scalar * dst, const Complex * src, int nfft) { - dst.resize( src.size() ); - inv( &dst[0],&src[0],src.size() ); + m_impl.inv( dst,src,nfft ); + if ( HasFlag( Unscaled ) == false) + scale(dst,1./nfft,nfft); } template<typename OutputDerived, typename ComplexDerived> + inline void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived) @@ -117,18 +156,52 @@ class FFT YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) 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) - dst.derived().resize( src.size() ); + + int nfft = src.size(); + int nout = HasFlag(HalfSpectrum) ? ((nfft>>1)+1) : nfft; + dst.derived().resize( nout ); inv( &dst[0],&src[0],src.size() ); } + template <typename _Output> + inline + void inv( std::vector<_Output> & dst, const std::vector<Complex> & src) + { + if ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) + dst.resize( 2*(src.size()-1) ); + else + dst.resize( src.size() ); + inv( &dst[0],&src[0],dst.size() ); + } + // TODO: multi-dimensional FFTs // TODO: handle Eigen MatrixBase // ---> i added fwd and inv specializations above + unit test, is this enough? (bjacob) + 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++ *= s; + } + + inline + void ReflectSpectrum(Complex * freq,int nfft) + { + // create the implicit right-half spectrum (conjugate-mirror of the left-half) + int nhbins=(nfft>>1)+1; + for (int k=nhbins;k < nfft; ++k ) + freq[k] = conj(freq[nfft-k]); + } + impl_type m_impl; + int m_flag; }; } #endif diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h index 2fb733a99..c4607c2b8 100644 --- a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h +++ b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h @@ -29,7 +29,7 @@ namespace Eigen { template<typename A, typename B> struct ei_make_coherent_impl { - static void run(A& a, B& b) {} + static void run(A&, B&) {} }; // resize a to match b is a.size()==0, and conversely. diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffVector.h b/unsupported/Eigen/src/AutoDiff/AutoDiffVector.h index 69ea9144e..03c82b7e8 100644 --- a/unsupported/Eigen/src/AutoDiff/AutoDiffVector.h +++ b/unsupported/Eigen/src/AutoDiff/AutoDiffVector.h @@ -35,7 +35,7 @@ namespace Eigen { * This class represents a scalar value while tracking its respective derivatives. * * It supports the following list of global math function: - * - std::abs, std::sqrt, std::pow, std::exp, std::log, std::sin, std::cos, + * - std::abs, std::sqrt, std::pow, std::exp, std::log, std::sin, std::cos, * - ei_abs, ei_sqrt, ei_pow, ei_exp, ei_log, ei_sin, ei_cos, * - ei_conj, ei_real, ei_imag, ei_abs2. * @@ -48,130 +48,150 @@ template<typename ValueType, typename JacobianType> class AutoDiffVector { public: - typedef typename ei_traits<ValueType>::Scalar Scalar; - + //typedef typename ei_traits<ValueType>::Scalar Scalar; + typedef typename ei_traits<ValueType>::Scalar BaseScalar; + typedef AutoDiffScalar<Matrix<BaseScalar,JacobianType::RowsAtCompileTime,1> > ActiveScalar; + typedef ActiveScalar Scalar; + typedef AutoDiffScalar<typename JacobianType::ColXpr> CoeffType; + inline AutoDiffVector() {} - + inline AutoDiffVector(const ValueType& values) : m_values(values) { m_jacobian.setZero(); } - + + + CoeffType operator[] (int i) { return CoeffType(m_values[i], m_jacobian.col(i)); } + const CoeffType operator[] (int i) const { return CoeffType(m_values[i], m_jacobian.col(i)); } + + CoeffType operator() (int i) { return CoeffType(m_values[i], m_jacobian.col(i)); } + const CoeffType operator() (int i) const { return CoeffType(m_values[i], m_jacobian.col(i)); } + + CoeffType coeffRef(int i) { return CoeffType(m_values[i], m_jacobian.col(i)); } + const CoeffType coeffRef(int i) const { return CoeffType(m_values[i], m_jacobian.col(i)); } + + int size() const { return m_values.size(); } + + // FIXME here we could return an expression of the sum + Scalar sum() const { /*std::cerr << "sum \n\n";*/ /*std::cerr << m_jacobian.rowwise().sum() << "\n\n";*/ return Scalar(m_values.sum(), m_jacobian.rowwise().sum()); } + + inline AutoDiffVector(const ValueType& values, const JacobianType& jac) : m_values(values), m_jacobian(jac) {} - + template<typename OtherValueType, typename OtherJacobianType> inline AutoDiffVector(const AutoDiffVector<OtherValueType, OtherJacobianType>& other) : m_values(other.values()), m_jacobian(other.jacobian()) {} - + inline AutoDiffVector(const AutoDiffVector& other) : m_values(other.values()), m_jacobian(other.jacobian()) {} - + template<typename OtherValueType, typename OtherJacobianType> - inline AutoDiffScalar& operator=(const AutoDiffVector<OtherValueType, OtherJacobianType>& other) + inline AutoDiffVector& operator=(const AutoDiffVector<OtherValueType, OtherJacobianType>& other) { m_values = other.values(); m_jacobian = other.jacobian(); return *this; } - + inline AutoDiffVector& operator=(const AutoDiffVector& other) { m_values = other.values(); m_jacobian = other.jacobian(); return *this; } - + inline const ValueType& values() const { return m_values; } inline ValueType& values() { return m_values; } - + inline const JacobianType& jacobian() const { return m_jacobian; } inline JacobianType& jacobian() { return m_jacobian; } - + template<typename OtherValueType,typename OtherJacobianType> inline const AutoDiffVector< - CwiseBinaryOp<ei_scalar_sum_op<Scalar>,ValueType,OtherValueType> > - CwiseBinaryOp<ei_scalar_sum_op<Scalar>,JacobianType,OtherJacobianType> > - operator+(const AutoDiffScalar<OtherDerType>& other) const + typename MakeCwiseBinaryOp<ei_scalar_sum_op<BaseScalar>,ValueType,OtherValueType>::Type, + typename MakeCwiseBinaryOp<ei_scalar_sum_op<BaseScalar>,JacobianType,OtherJacobianType>::Type > + operator+(const AutoDiffVector<OtherValueType,OtherJacobianType>& other) const { return AutoDiffVector< - CwiseBinaryOp<ei_scalar_sum_op<Scalar>,ValueType,OtherValueType> > - CwiseBinaryOp<ei_scalar_sum_op<Scalar>,JacobianType,OtherJacobianType> >( + typename MakeCwiseBinaryOp<ei_scalar_sum_op<BaseScalar>,ValueType,OtherValueType>::Type, + typename MakeCwiseBinaryOp<ei_scalar_sum_op<BaseScalar>,JacobianType,OtherJacobianType>::Type >( m_values + other.values(), m_jacobian + other.jacobian()); } - + template<typename OtherValueType, typename OtherJacobianType> inline AutoDiffVector& - operator+=(const AutoDiffVector<OtherValueType,OtherDerType>& other) + operator+=(const AutoDiffVector<OtherValueType,OtherJacobianType>& other) { m_values += other.values(); m_jacobian += other.jacobian(); return *this; } - + template<typename OtherValueType,typename OtherJacobianType> inline const AutoDiffVector< - CwiseBinaryOp<ei_scalar_difference_op<Scalar>,ValueType,OtherValueType> > - CwiseBinaryOp<ei_scalar_difference_op<Scalar>,JacobianType,OtherJacobianType> > - operator-(const AutoDiffScalar<OtherDerType>& other) const + typename MakeCwiseBinaryOp<ei_scalar_difference_op<Scalar>,ValueType,OtherValueType>::Type, + typename MakeCwiseBinaryOp<ei_scalar_difference_op<Scalar>,JacobianType,OtherJacobianType>::Type > + operator-(const AutoDiffVector<OtherValueType,OtherJacobianType>& other) const { return AutoDiffVector< - CwiseBinaryOp<ei_scalar_difference_op<Scalar>,ValueType,OtherValueType> > - CwiseBinaryOp<ei_scalar_difference_op<Scalar>,JacobianType,OtherJacobianType> >( - m_values - other.values(), - m_jacobian - other.jacobian()); + typename MakeCwiseBinaryOp<ei_scalar_difference_op<Scalar>,ValueType,OtherValueType>::Type, + typename MakeCwiseBinaryOp<ei_scalar_difference_op<Scalar>,JacobianType,OtherJacobianType>::Type >( + m_values - other.values(), + m_jacobian - other.jacobian()); } - + template<typename OtherValueType, typename OtherJacobianType> inline AutoDiffVector& - operator-=(const AutoDiffVector<OtherValueType,OtherDerType>& other) + operator-=(const AutoDiffVector<OtherValueType,OtherJacobianType>& other) { m_values -= other.values(); m_jacobian -= other.jacobian(); return *this; } - + inline const AutoDiffVector< - CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, ValueType> - CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, JacobianType> > + typename MakeCwiseUnaryOp<ei_scalar_opposite_op<Scalar>, ValueType>::Type, + typename MakeCwiseUnaryOp<ei_scalar_opposite_op<Scalar>, JacobianType>::Type > operator-() const { return AutoDiffVector< - CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, ValueType> - CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, JacobianType> >( - -m_values, - -m_jacobian); + typename MakeCwiseUnaryOp<ei_scalar_opposite_op<Scalar>, ValueType>::Type, + typename MakeCwiseUnaryOp<ei_scalar_opposite_op<Scalar>, JacobianType>::Type >( + -m_values, + -m_jacobian); } - + inline const AutoDiffVector< - CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType> - CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> > - operator*(const Scalar& other) const + typename MakeCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType>::Type, + typename MakeCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType>::Type> + operator*(const BaseScalar& other) const { return AutoDiffVector< - CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType> - CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> >( + typename MakeCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType>::Type, + typename MakeCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType>::Type >( m_values * other, - (m_jacobian * other)); + m_jacobian * other); } - + friend inline const AutoDiffVector< - CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType> - CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> > + typename MakeCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType>::Type, + typename MakeCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType>::Type > operator*(const Scalar& other, const AutoDiffVector& v) { return AutoDiffVector< - CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType> - CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> >( + typename MakeCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType>::Type, + typename MakeCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType>::Type >( v.values() * other, v.jacobian() * other); } - + // template<typename OtherValueType,typename OtherJacobianType> // inline const AutoDiffVector< // CwiseBinaryOp<ei_scalar_multiple_op<Scalar>, ValueType, OtherValueType> @@ -188,25 +208,25 @@ class AutoDiffVector // m_values.cwise() * other.values(), // (m_jacobian * other.values()).nestByValue() + (m_values * other.jacobian()).nestByValue()); // } - + inline AutoDiffVector& operator*=(const Scalar& other) { m_values *= other; m_jacobian *= other; return *this; } - + template<typename OtherValueType,typename OtherJacobianType> inline AutoDiffVector& operator*=(const AutoDiffVector<OtherValueType,OtherJacobianType>& other) { *this = *this * other; return *this; } - + protected: ValueType m_values; JacobianType m_jacobian; - + }; } diff --git a/unsupported/Eigen/src/FFT/ei_fftw_impl.h b/unsupported/Eigen/src/FFT/ei_fftw_impl.h index e1f67f334..a66b7398c 100644 --- a/unsupported/Eigen/src/FFT/ei_fftw_impl.h +++ b/unsupported/Eigen/src/FFT/ei_fftw_impl.h @@ -166,6 +166,7 @@ m_plans.clear(); } + // complex-to-complex forward FFT inline void fwd( Complex * dst,const Complex *src,int nfft) { @@ -177,9 +178,6 @@ 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); - int nhbins=(nfft>>1)+1; - for (int k=nhbins;k < nfft; ++k ) - dst[k] = conj(dst[nfft-k]); } // inverse complex-to-complex @@ -187,12 +185,6 @@ 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 = Scalar(1.)/nfft; - for (int k=0;k<nfft;++k) - dst[k] *= s; } // half-complex to scalar @@ -200,11 +192,6 @@ 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 ); - - //TODO move scaling to Eigen::FFT - Scalar s = Scalar(1.)/nfft; - for (int k=0;k<nfft;++k) - dst[k] *= s; } protected: @@ -222,3 +209,5 @@ return m_plans[key]; } }; +/* vim: set filetype=cpp et sw=2 ts=2 ai: */ + diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h index c068d8765..5c958d1ec 100644 --- a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h +++ b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h @@ -27,388 +27,384 @@ // This FFT implementation was derived from kissfft http:sourceforge.net/projects/kissfft // Copyright 2003-2009 Mark Borgerding - template <typename _Scalar> - struct ei_kiss_cpx_fft +template <typename _Scalar> +struct ei_kiss_cpx_fft +{ + typedef _Scalar Scalar; + typedef std::complex<Scalar> Complex; + std::vector<Complex> m_twiddles; + std::vector<int> m_stageRadix; + std::vector<int> m_stageRemainder; + std::vector<Complex> m_scratchBuf; + bool m_inverse; + + inline + void make_twiddles(int nfft,bool inverse) { - typedef _Scalar Scalar; - typedef std::complex<Scalar> Complex; - std::vector<Complex> m_twiddles; - std::vector<int> m_stageRadix; - std::vector<int> m_stageRemainder; - std::vector<Complex> m_scratchBuf; - bool m_inverse; - - void make_twiddles(int nfft,bool inverse) - { - m_inverse = inverse; - m_twiddles.resize(nfft); - Scalar phinc = (inverse?2:-2)* acos( (Scalar) -1) / nfft; - for (int i=0;i<nfft;++i) - m_twiddles[i] = exp( Complex(0,i*phinc) ); + m_inverse = inverse; + m_twiddles.resize(nfft); + Scalar phinc = (inverse?2:-2)* acos( (Scalar) -1) / nfft; + for (int i=0;i<nfft;++i) + m_twiddles[i] = exp( Complex(0,i*phinc) ); + } + + void factorize(int nfft) + { + //start factoring out 4's, then 2's, then 3,5,7,9,... + int n= nfft; + int p=4; + do { + while (n % p) { + switch (p) { + case 4: p = 2; break; + case 2: p = 3; break; + default: p += 2; break; + } + if (p*p>n) + p=n;// impossible to have a factor > sqrt(n) } + n /= p; + m_stageRadix.push_back(p); + m_stageRemainder.push_back(n); + if ( p > 5 ) + m_scratchBuf.resize(p); // scratchbuf will be needed in bfly_generic + }while(n>1); + } + + template <typename _Src> + inline + void work( int stage,Complex * xout, const _Src * xin, size_t fstride,size_t in_stride) + { + int p = m_stageRadix[stage]; + int m = m_stageRemainder[stage]; + Complex * Fout_beg = xout; + Complex * Fout_end = xout + p*m; - void factorize(int nfft) - { - //start factoring out 4's, then 2's, then 3,5,7,9,... - int n= nfft; - int p=4; - do { - while (n % p) { - switch (p) { - case 4: p = 2; break; - case 2: p = 3; break; - default: p += 2; break; - } - if (p*p>n) - p=n;// impossible to have a factor > sqrt(n) - } - n /= p; - m_stageRadix.push_back(p); - m_stageRemainder.push_back(n); - if ( p > 5 ) - m_scratchBuf.resize(p); // scratchbuf will be needed in bfly_generic - }while(n>1); + if (m>1) { + do{ + // recursive call: + // DFT of size m*p performed by doing + // p instances of smaller DFTs of size m, + // each one takes a decimated version of the input + work(stage+1, xout , xin, fstride*p,in_stride); + xin += fstride*in_stride; + }while( (xout += m) != Fout_end ); + }else{ + do{ + *xout = *xin; + xin += fstride*in_stride; + }while(++xout != Fout_end ); } - - template <typename _Src> - void work( int stage,Complex * xout, const _Src * xin, size_t fstride,size_t in_stride) - { - int p = m_stageRadix[stage]; - int m = m_stageRemainder[stage]; - Complex * Fout_beg = xout; - Complex * Fout_end = xout + p*m; - - if (m>1) { - do{ - // recursive call: - // DFT of size m*p performed by doing - // p instances of smaller DFTs of size m, - // each one takes a decimated version of the input - work(stage+1, xout , xin, fstride*p,in_stride); - xin += fstride*in_stride; - }while( (xout += m) != Fout_end ); - }else{ - do{ - *xout = *xin; - xin += fstride*in_stride; - }while(++xout != Fout_end ); - } - xout=Fout_beg; - - // recombine the p smaller DFTs - switch (p) { - case 2: bfly2(xout,fstride,m); break; - case 3: bfly3(xout,fstride,m); break; - case 4: bfly4(xout,fstride,m); break; - case 5: bfly5(xout,fstride,m); break; - default: bfly_generic(xout,fstride,m,p); break; - } - } - - inline - void bfly2( Complex * Fout, const size_t fstride, int m) - { - for (int k=0;k<m;++k) { - Complex t = Fout[m+k] * m_twiddles[k*fstride]; - Fout[m+k] = Fout[k] - t; - Fout[k] += t; - } + xout=Fout_beg; + + // recombine the p smaller DFTs + switch (p) { + case 2: bfly2(xout,fstride,m); break; + case 3: bfly3(xout,fstride,m); break; + case 4: bfly4(xout,fstride,m); break; + case 5: bfly5(xout,fstride,m); break; + default: bfly_generic(xout,fstride,m,p); break; } + } - inline - void bfly4( Complex * Fout, const size_t fstride, const size_t m) - { - Complex scratch[6]; - int negative_if_inverse = m_inverse * -2 +1; - for (size_t k=0;k<m;++k) { - scratch[0] = Fout[k+m] * m_twiddles[k*fstride]; - scratch[1] = Fout[k+2*m] * m_twiddles[k*fstride*2]; - scratch[2] = Fout[k+3*m] * m_twiddles[k*fstride*3]; - scratch[5] = Fout[k] - scratch[1]; - - Fout[k] += scratch[1]; - scratch[3] = scratch[0] + scratch[2]; - scratch[4] = scratch[0] - scratch[2]; - scratch[4] = Complex( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse ); - - Fout[k+2*m] = Fout[k] - scratch[3]; - Fout[k] += scratch[3]; - Fout[k+m] = scratch[5] + scratch[4]; - Fout[k+3*m] = scratch[5] - scratch[4]; - } + inline + void bfly2( Complex * Fout, const size_t fstride, int m) + { + for (int k=0;k<m;++k) { + Complex t = Fout[m+k] * m_twiddles[k*fstride]; + Fout[m+k] = Fout[k] - t; + Fout[k] += t; } + } - inline - void bfly3( Complex * Fout, const size_t fstride, const size_t m) - { - size_t k=m; - const size_t m2 = 2*m; - Complex *tw1,*tw2; - Complex scratch[5]; - Complex epi3; - epi3 = m_twiddles[fstride*m]; - - tw1=tw2=&m_twiddles[0]; - - do{ - scratch[1]=Fout[m] * *tw1; - scratch[2]=Fout[m2] * *tw2; - - scratch[3]=scratch[1]+scratch[2]; - scratch[0]=scratch[1]-scratch[2]; - tw1 += fstride; - tw2 += fstride*2; - Fout[m] = Complex( Fout->real() - .5*scratch[3].real() , Fout->imag() - .5*scratch[3].imag() ); - scratch[0] *= epi3.imag(); - *Fout += scratch[3]; - Fout[m2] = Complex( Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() ); - Fout[m] += Complex( -scratch[0].imag(),scratch[0].real() ); - ++Fout; - }while(--k); + inline + void bfly4( Complex * Fout, const size_t fstride, const size_t m) + { + Complex scratch[6]; + int negative_if_inverse = m_inverse * -2 +1; + for (size_t k=0;k<m;++k) { + scratch[0] = Fout[k+m] * m_twiddles[k*fstride]; + scratch[1] = Fout[k+2*m] * m_twiddles[k*fstride*2]; + scratch[2] = Fout[k+3*m] * m_twiddles[k*fstride*3]; + scratch[5] = Fout[k] - scratch[1]; + + Fout[k] += scratch[1]; + scratch[3] = scratch[0] + scratch[2]; + scratch[4] = scratch[0] - scratch[2]; + scratch[4] = Complex( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse ); + + Fout[k+2*m] = Fout[k] - scratch[3]; + Fout[k] += scratch[3]; + Fout[k+m] = scratch[5] + scratch[4]; + Fout[k+3*m] = scratch[5] - scratch[4]; } + } - inline - void bfly5( Complex * Fout, const size_t fstride, const size_t m) - { - Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; - size_t u; - Complex scratch[13]; - Complex * twiddles = &m_twiddles[0]; - Complex *tw; - Complex ya,yb; - ya = twiddles[fstride*m]; - yb = twiddles[fstride*2*m]; - - Fout0=Fout; - Fout1=Fout0+m; - Fout2=Fout0+2*m; - Fout3=Fout0+3*m; - Fout4=Fout0+4*m; - - tw=twiddles; - for ( u=0; u<m; ++u ) { - scratch[0] = *Fout0; - - scratch[1] = *Fout1 * tw[u*fstride]; - scratch[2] = *Fout2 * tw[2*u*fstride]; - scratch[3] = *Fout3 * tw[3*u*fstride]; - scratch[4] = *Fout4 * tw[4*u*fstride]; - - scratch[7] = scratch[1] + scratch[4]; - scratch[10] = scratch[1] - scratch[4]; - scratch[8] = scratch[2] + scratch[3]; - scratch[9] = scratch[2] - scratch[3]; - - *Fout0 += scratch[7]; - *Fout0 += scratch[8]; - - scratch[5] = scratch[0] + Complex( - (scratch[7].real()*ya.real() ) + (scratch[8].real() *yb.real() ), - (scratch[7].imag()*ya.real()) + (scratch[8].imag()*yb.real()) - ); - - scratch[6] = Complex( - (scratch[10].imag()*ya.imag()) + (scratch[9].imag()*yb.imag()), - -(scratch[10].real()*ya.imag()) - (scratch[9].real()*yb.imag()) + inline + void bfly3( Complex * Fout, const size_t fstride, const size_t m) + { + size_t k=m; + const size_t m2 = 2*m; + Complex *tw1,*tw2; + Complex scratch[5]; + Complex epi3; + epi3 = m_twiddles[fstride*m]; + + tw1=tw2=&m_twiddles[0]; + + do{ + scratch[1]=Fout[m] * *tw1; + scratch[2]=Fout[m2] * *tw2; + + scratch[3]=scratch[1]+scratch[2]; + scratch[0]=scratch[1]-scratch[2]; + tw1 += fstride; + tw2 += fstride*2; + Fout[m] = Complex( Fout->real() - .5*scratch[3].real() , Fout->imag() - .5*scratch[3].imag() ); + scratch[0] *= epi3.imag(); + *Fout += scratch[3]; + Fout[m2] = Complex( Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() ); + Fout[m] += Complex( -scratch[0].imag(),scratch[0].real() ); + ++Fout; + }while(--k); + } + + inline + void bfly5( Complex * Fout, const size_t fstride, const size_t m) + { + Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; + size_t u; + Complex scratch[13]; + Complex * twiddles = &m_twiddles[0]; + Complex *tw; + Complex ya,yb; + ya = twiddles[fstride*m]; + yb = twiddles[fstride*2*m]; + + Fout0=Fout; + Fout1=Fout0+m; + Fout2=Fout0+2*m; + Fout3=Fout0+3*m; + Fout4=Fout0+4*m; + + tw=twiddles; + for ( u=0; u<m; ++u ) { + scratch[0] = *Fout0; + + scratch[1] = *Fout1 * tw[u*fstride]; + scratch[2] = *Fout2 * tw[2*u*fstride]; + scratch[3] = *Fout3 * tw[3*u*fstride]; + scratch[4] = *Fout4 * tw[4*u*fstride]; + + scratch[7] = scratch[1] + scratch[4]; + scratch[10] = scratch[1] - scratch[4]; + scratch[8] = scratch[2] + scratch[3]; + scratch[9] = scratch[2] - scratch[3]; + + *Fout0 += scratch[7]; + *Fout0 += scratch[8]; + + scratch[5] = scratch[0] + Complex( + (scratch[7].real()*ya.real() ) + (scratch[8].real() *yb.real() ), + (scratch[7].imag()*ya.real()) + (scratch[8].imag()*yb.real()) + ); + + scratch[6] = Complex( + (scratch[10].imag()*ya.imag()) + (scratch[9].imag()*yb.imag()), + -(scratch[10].real()*ya.imag()) - (scratch[9].real()*yb.imag()) + ); + + *Fout1 = scratch[5] - scratch[6]; + *Fout4 = scratch[5] + scratch[6]; + + scratch[11] = scratch[0] + + Complex( + (scratch[7].real()*yb.real()) + (scratch[8].real()*ya.real()), + (scratch[7].imag()*yb.real()) + (scratch[8].imag()*ya.real()) ); - *Fout1 = scratch[5] - scratch[6]; - *Fout4 = scratch[5] + scratch[6]; - - scratch[11] = scratch[0] + - Complex( - (scratch[7].real()*yb.real()) + (scratch[8].real()*ya.real()), - (scratch[7].imag()*yb.real()) + (scratch[8].imag()*ya.real()) - ); + scratch[12] = Complex( + -(scratch[10].imag()*yb.imag()) + (scratch[9].imag()*ya.imag()), + (scratch[10].real()*yb.imag()) - (scratch[9].real()*ya.imag()) + ); - scratch[12] = Complex( - -(scratch[10].imag()*yb.imag()) + (scratch[9].imag()*ya.imag()), - (scratch[10].real()*yb.imag()) - (scratch[9].real()*ya.imag()) - ); - - *Fout2=scratch[11]+scratch[12]; - *Fout3=scratch[11]-scratch[12]; + *Fout2=scratch[11]+scratch[12]; + *Fout3=scratch[11]-scratch[12]; - ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; - } + ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; } + } + + /* perform the butterfly for one stage of a mixed radix FFT */ + inline + void bfly_generic( + Complex * Fout, + const size_t fstride, + int m, + int p + ) + { + int u,k,q1,q; + Complex * twiddles = &m_twiddles[0]; + Complex t; + int Norig = m_twiddles.size(); + Complex * scratchbuf = &m_scratchBuf[0]; + + for ( u=0; u<m; ++u ) { + k=u; + for ( q1=0 ; q1<p ; ++q1 ) { + scratchbuf[q1] = Fout[ k ]; + k += m; + } - /* perform the butterfly for one stage of a mixed radix FFT */ - inline - void bfly_generic( - Complex * Fout, - const size_t fstride, - int m, - int p - ) - { - int u,k,q1,q; - Complex * twiddles = &m_twiddles[0]; - Complex t; - int Norig = m_twiddles.size(); - Complex * scratchbuf = &m_scratchBuf[0]; - - for ( u=0; u<m; ++u ) { - k=u; - for ( q1=0 ; q1<p ; ++q1 ) { - scratchbuf[q1] = Fout[ k ]; - k += m; - } - - k=u; - for ( q1=0 ; q1<p ; ++q1 ) { - int twidx=0; - Fout[ k ] = scratchbuf[0]; - for (q=1;q<p;++q ) { - twidx += fstride * k; - if (twidx>=Norig) twidx-=Norig; - t=scratchbuf[q] * twiddles[twidx]; - Fout[ k ] += t; - } - k += m; + k=u; + for ( q1=0 ; q1<p ; ++q1 ) { + int twidx=0; + Fout[ k ] = scratchbuf[0]; + for (q=1;q<p;++q ) { + twidx += fstride * k; + if (twidx>=Norig) twidx-=Norig; + t=scratchbuf[q] * twiddles[twidx]; + Fout[ k ] += t; } + k += m; } } - }; - - template <typename _Scalar> - struct ei_kissfft_impl + } +}; + +template <typename _Scalar> +struct ei_kissfft_impl +{ + typedef _Scalar Scalar; + typedef std::complex<Scalar> Complex; + + void clear() + { + m_plans.clear(); + m_realTwiddles.clear(); + } + + inline + void fwd( Complex * dst,const Complex *src,int nfft) { - typedef _Scalar Scalar; - typedef std::complex<Scalar> Complex; - - void clear() - { - m_plans.clear(); - m_realTwiddles.clear(); - } - - template <typename _Src> - inline - void fwd( Complex * dst,const _Src *src,int nfft) - { - get_plan(nfft,false).work(0, dst, src, 1,1); + get_plan(nfft,false).work(0, dst, src, 1,1); + } + + // real-to-complex forward FFT + // 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 ) { + // use generic mode for odd + m_tmpBuf1.resize(nfft); + get_plan(nfft,false).work(0, &m_tmpBuf1[0], src, 1,1); + std::copy(m_tmpBuf1.begin(),m_tmpBuf1.begin()+(nfft>>1)+1,dst ); + }else{ + int ncfft = nfft>>1; + int ncfft2 = nfft>>2; + Complex * rtw = real_twiddles(ncfft2); + + // use optimized mode for even real + fwd( dst, reinterpret_cast<const Complex*> (src), ncfft); + Complex dc = dst[0].real() + dst[0].imag(); + Complex nyquist = dst[0].real() - dst[0].imag(); + int k; + for ( k=1;k <= ncfft2 ; ++k ) { + Complex fpk = dst[k]; + Complex fpnk = conj(dst[ncfft-k]); + Complex f1k = fpk + fpnk; + Complex f2k = fpk - fpnk; + Complex tw= f2k * rtw[k-1]; + dst[k] = (f1k + tw) * Scalar(.5); + dst[ncfft-k] = conj(f1k -tw)*Scalar(.5); } + dst[0] = dc; + dst[ncfft] = nyquist; + } + } - // real-to-complex forward FFT - // 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 ) { - // use generic mode for odd - get_plan(nfft,false).work(0, dst, src, 1,1); - }else{ - int ncfft = nfft>>1; - int ncfft2 = nfft>>2; - Complex * rtw = real_twiddles(ncfft2); - - // use optimized mode for even real - fwd( dst, reinterpret_cast<const Complex*> (src), ncfft); - Complex dc = dst[0].real() + dst[0].imag(); - Complex nyquist = dst[0].real() - dst[0].imag(); - int k; - for ( k=1;k <= ncfft2 ; ++k ) { - Complex fpk = dst[k]; - Complex fpnk = conj(dst[ncfft-k]); - Complex f1k = fpk + fpnk; - Complex f2k = fpk - fpnk; - Complex tw= f2k * rtw[k-1]; - dst[k] = (f1k + tw) * Scalar(.5); - dst[ncfft-k] = conj(f1k -tw)*Scalar(.5); - } + // inverse complex-to-complex + inline + void inv(Complex * dst,const Complex *src,int nfft) + { + get_plan(nfft,true).work(0, dst, src, 1,1); + } - // place conjugate-symmetric half at the end for completeness - // TODO: make this configurable ( opt-out ) - for ( k=1;k < ncfft ; ++k ) - dst[nfft-k] = conj(dst[k]); - dst[0] = dc; - dst[ncfft] = nyquist; + // half-complex to scalar + inline + void inv( Scalar * dst,const Complex * src,int nfft) + { + if (nfft&3) { + m_tmpBuf1.resize(nfft); + m_tmpBuf2.resize(nfft); + std::copy(src,src+(nfft>>1)+1,m_tmpBuf1.begin() ); + for (int k=1;k<(nfft>>1)+1;++k) + m_tmpBuf1[nfft-k] = conj(m_tmpBuf1[k]); + inv(&m_tmpBuf2[0],&m_tmpBuf1[0],nfft); + for (int k=0;k<nfft;++k) + dst[k] = m_tmpBuf2[k].real(); + }else{ + // optimized version for multiple of 4 + int ncfft = nfft>>1; + int ncfft2 = nfft>>2; + Complex * rtw = real_twiddles(ncfft2); + m_tmpBuf1.resize(ncfft); + m_tmpBuf1[0] = Complex( src[0].real() + src[ncfft].real(), src[0].real() - src[ncfft].real() ); + for (int k = 1; k <= ncfft / 2; ++k) { + Complex fk = src[k]; + Complex fnkc = conj(src[ncfft-k]); + Complex fek = fk + fnkc; + Complex tmp = fk - fnkc; + Complex fok = tmp * conj(rtw[k-1]); + m_tmpBuf1[k] = fek + fok; + m_tmpBuf1[ncfft-k] = conj(fek - fok); } + get_plan(ncfft,true).work(0, reinterpret_cast<Complex*>(dst), &m_tmpBuf1[0], 1,1); } + } - // inverse complex-to-complex - inline - void inv(Complex * dst,const Complex *src,int nfft) - { - get_plan(nfft,true).work(0, dst, src, 1,1); - scale(dst, nfft, Scalar(1)/nfft ); - } + protected: + typedef ei_kiss_cpx_fft<Scalar> PlanData; + typedef std::map<int,PlanData> PlanMap; - // half-complex to scalar - inline - void inv( Scalar * dst,const Complex * src,int nfft) - { - if (nfft&3) { - m_tmpBuf.resize(nfft); - inv(&m_tmpBuf[0],src,nfft); - for (int k=0;k<nfft;++k) - dst[k] = m_tmpBuf[k].real(); - }else{ - // optimized version for multiple of 4 - int ncfft = nfft>>1; - int ncfft2 = nfft>>2; - Complex * rtw = real_twiddles(ncfft2); - m_tmpBuf.resize(ncfft); - m_tmpBuf[0] = Complex( src[0].real() + src[ncfft].real(), src[0].real() - src[ncfft].real() ); - for (int k = 1; k <= ncfft / 2; ++k) { - Complex fk = src[k]; - Complex fnkc = conj(src[ncfft-k]); - Complex fek = fk + fnkc; - Complex tmp = fk - fnkc; - Complex fok = tmp * conj(rtw[k-1]); - m_tmpBuf[k] = fek + fok; - m_tmpBuf[ncfft-k] = conj(fek - fok); - } - scale(&m_tmpBuf[0], ncfft, Scalar(1)/nfft ); - get_plan(ncfft,true).work(0, reinterpret_cast<Complex*>(dst), &m_tmpBuf[0], 1,1); - } - } + PlanMap m_plans; + std::map<int, std::vector<Complex> > m_realTwiddles; + std::vector<Complex> m_tmpBuf1; + std::vector<Complex> m_tmpBuf2; - protected: - typedef ei_kiss_cpx_fft<Scalar> PlanData; - typedef std::map<int,PlanData> PlanMap; - - PlanMap m_plans; - 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 - PlanData & pd = m_plans[ PlanKey(nfft,inverse) ]; - if ( pd.m_twiddles.size() == 0 ) { - pd.make_twiddles(nfft,inverse); - pd.factorize(nfft); - } - return pd; - } + inline + int PlanKey(int nfft,bool isinverse) const { return (nfft<<1) | isinverse; } - inline - Complex * real_twiddles(int ncfft2) - { - std::vector<Complex> & twidref = m_realTwiddles[ncfft2];// creates new if not there - if ( (int)twidref.size() != ncfft2 ) { - twidref.resize(ncfft2); - int ncfft= ncfft2<<1; - Scalar pi = acos( Scalar(-1) ); - for (int k=1;k<=ncfft2;++k) - twidref[k-1] = exp( Complex(0,-pi * ((double) (k) / ncfft + .5) ) ); - } - return &twidref[0]; + inline + PlanData & get_plan(int nfft,bool inverse) + { + // TODO look for PlanKey(nfft, ! inverse) and conjugate the twiddles + PlanData & pd = m_plans[ PlanKey(nfft,inverse) ]; + if ( pd.m_twiddles.size() == 0 ) { + pd.make_twiddles(nfft,inverse); + pd.factorize(nfft); } + return pd; + } - // 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; + inline + Complex * real_twiddles(int ncfft2) + { + std::vector<Complex> & twidref = m_realTwiddles[ncfft2];// creates new if not there + if ( (int)twidref.size() != ncfft2 ) { + twidref.resize(ncfft2); + int ncfft= ncfft2<<1; + Scalar pi = acos( Scalar(-1) ); + for (int k=1;k<=ncfft2;++k) + twidref[k-1] = exp( Complex(0,-pi * ((double) (k) / ncfft + .5) ) ); } - }; + return &twidref[0]; + } +}; + +/* vim: set filetype=cpp et sw=2 ts=2 ai: */ + diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index d182c9abf..c75d5eb0b 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -26,3 +26,4 @@ if(FFTW_FOUND) ei_add_test(FFTW "-DEIGEN_FFTW_DEFAULT " "-lfftw3 -lfftw3f -lfftw3l" ) endif(FFTW_FOUND) +ei_add_test(Complex) diff --git a/unsupported/test/Complex.cpp b/unsupported/test/Complex.cpp new file mode 100644 index 000000000..bedeb9f27 --- /dev/null +++ b/unsupported/test/Complex.cpp @@ -0,0 +1,77 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Mark Borgerding mark a borgerding net +// +// Eigen 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. +// +// Alternatively, 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. +// +// Eigen 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 or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. +#ifdef EIGEN_TEST_FUNC +# include "main.h" +#else +# include <iostream> +# define CALL_SUBTEST(x) x +# define VERIFY(x) x +# define test_Complex main +#endif + +#include <unsupported/Eigen/Complex> +#include <vector> + +using namespace std; +using namespace Eigen; + +template <typename T> +void take_std( std::complex<T> * dst, int n ) +{ + cout << dst[n-1] << endl; +} + + +template <typename T> +void syntax() +{ + // this works fine + Matrix< Complex<T>, 9, 1> a; + std::complex<T> * pa = &a[0]; + Complex<T> * pa2 = &a[0]; + take_std( pa,9); + + // this does not work, but I wish it would + // take_std(&a[0];) + // this does + take_std( (std::complex<T> *)&a[0],9); + + // this does not work, but it would be really nice + //vector< Complex<T> > a; + // (on my gcc 4.4.1 ) + // std::vector assumes operator& returns a POD pointer + + // this works fine + Complex<T> b[9]; + std::complex<T> * pb = &b[0]; // this works fine + + take_std( pb,9); +} + +void test_Complex() +{ + CALL_SUBTEST( syntax<float>() ); + CALL_SUBTEST( syntax<double>() ); + CALL_SUBTEST( syntax<long double>() ); +} diff --git a/unsupported/test/FFT.cpp b/unsupported/test/FFT.cpp index cc68f3718..ad0d426e4 100644 --- a/unsupported/test/FFT.cpp +++ b/unsupported/test/FFT.cpp @@ -101,12 +101,34 @@ void test_scalar_generic(int nfft) ComplexVector outbuf; for (int k=0;k<nfft;++k) inbuf[k]= (T)(rand()/(double)RAND_MAX - .5); + + // 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); + VERIFY(outbuf.size() == (nfft>>1)+1); + VERIFY( fft_rmse(outbuf,inbuf) < test_precision<T>() );// gross check + + fft.ClearFlag(fft.HalfSpectrum ); fft.fwd( outbuf,inbuf); VERIFY( fft_rmse(outbuf,inbuf) < test_precision<T>() );// gross check ScalarVector buf3; fft.inv( buf3 , outbuf); VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>() );// gross check + + // verify that the Unscaled flag takes effect + ComplexVector buf4; + fft.SetFlag(fft.Unscaled); + fft.inv( buf4 , outbuf); + for (int k=0;k<nfft;++k) + buf4[k] *= T(1./nfft); + VERIFY( dif_rmse(inbuf,buf4) < test_precision<T>() );// gross check + + // verify that ClearFlag works + fft.ClearFlag(fft.Unscaled); + fft.inv( buf3 , outbuf); + VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>() );// gross check } template <typename T> @@ -136,6 +158,19 @@ void test_complex_generic(int nfft) fft.inv( buf3 , outbuf); VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>() );// gross check + + // verify that the Unscaled flag takes effect + ComplexVector buf4; + fft.SetFlag(fft.Unscaled); + fft.inv( buf4 , outbuf); + for (int k=0;k<nfft;++k) + buf4[k] *= T(1./nfft); + VERIFY( dif_rmse(inbuf,buf4) < test_precision<T>() );// gross check + + // verify that ClearFlag works + fft.ClearFlag(fft.Unscaled); + fft.inv( buf3 , outbuf); + VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>() );// gross check } template <typename T> |