diff --git a/bench/benchFFT.cpp b/bench/benchFFT.cpp
new file mode 100644
index 000000000..4b6cabb55
--- /dev/null
+++ b/bench/benchFFT.cpp
@@ -0,0 +1,115 @@
+// 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 .
+
+#include
+#include
+#include
+#include
+#ifdef USE_FFTW
+#include
+#endif
+
+#include
+
+using namespace Eigen;
+using namespace std;
+
+
+template
+string nameof();
+
+template <> string nameof() {return "float";}
+template <> string nameof() {return "double";}
+template <> string nameof() {return "long double";}
+
+#ifndef TYPE
+#define TYPE float
+#endif
+
+#ifndef NFFT
+#define NFFT 1024
+#endif
+#ifndef NDATA
+#define NDATA 1000000
+#endif
+
+using namespace Eigen;
+
+template
+void bench(int nfft,bool fwd)
+{
+ typedef typename NumTraits::Real Scalar;
+ typedef typename std::complex Complex;
+ int nits = NDATA/nfft;
+ vector inbuf(nfft);
+ vector outbuf(nfft);
+ FFT< Scalar > fft;
+
+ fft.fwd( outbuf , inbuf);
+
+ BenchTimer timer;
+ timer.reset();
+ for (int k=0;k<8;++k) {
+ timer.start();
+ for(int i = 0; i < nits; i++)
+ if (fwd)
+ fft.fwd( outbuf , inbuf);
+ else
+ fft.inv(inbuf,outbuf);
+ timer.stop();
+ }
+
+ cout << nameof() << " ";
+ double mflops = 5.*nfft*log2((double)nfft) / (1e6 * timer.value() / (double)nits );
+ if ( NumTraits::IsComplex ) {
+ cout << "complex";
+ }else{
+ cout << "real ";
+ mflops /= 2;
+ }
+
+ if (fwd)
+ cout << " fwd";
+ else
+ cout << " inv";
+
+ cout << " NFFT=" << nfft << " " << (double(1e-6*nfft*nits)/timer.value()) << " MS/s " << mflops << "MFLOPS\n";
+}
+
+int main(int argc,char ** argv)
+{
+ bench >(NFFT,true);
+ bench >(NFFT,false);
+ bench(NFFT,true);
+ bench(NFFT,false);
+ bench >(NFFT,true);
+ bench >(NFFT,false);
+ bench(NFFT,true);
+ bench(NFFT,false);
+ bench >(NFFT,true);
+ bench >(NFFT,false);
+ bench(NFFT,true);
+ bench(NFFT,false);
+ return 0;
+}
diff --git a/cmake/FindFFTW.cmake b/cmake/FindFFTW.cmake
new file mode 100644
index 000000000..a56450b17
--- /dev/null
+++ b/cmake/FindFFTW.cmake
@@ -0,0 +1,24 @@
+
+if (FFTW_INCLUDES AND FFTW_LIBRARIES)
+ set(FFTW_FIND_QUIETLY TRUE)
+endif (FFTW_INCLUDES AND FFTW_LIBRARIES)
+
+find_path(FFTW_INCLUDES
+ NAMES
+ fftw3.h
+ PATHS
+ $ENV{FFTWDIR}
+ ${INCLUDE_INSTALL_DIR}
+)
+
+find_library(FFTWF_LIB NAMES fftw3f PATHS $ENV{FFTWDIR} ${LIB_INSTALL_DIR})
+find_library(FFTW_LIB NAMES fftw3 PATHS $ENV{FFTWDIR} ${LIB_INSTALL_DIR})
+find_library(FFTWL_LIB NAMES fftw3l PATHS $ENV{FFTWDIR} ${LIB_INSTALL_DIR})
+set(FFTW_LIBRARIES "${FFTWF_LIB} ${FFTW_LIB} ${FFTWL_LIB}" )
+message(STATUS "FFTW ${FFTW_LIBRARIES}" )
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(FFTW DEFAULT_MSG
+ FFTW_INCLUDES FFTW_LIBRARIES)
+
+mark_as_advanced(FFTW_INCLUDES FFTW_LIBRARIES)
diff --git a/unsupported/Eigen/Complex b/unsupported/Eigen/Complex
new file mode 100644
index 000000000..a0483abdf
--- /dev/null
+++ b/unsupported/Eigen/Complex
@@ -0,0 +1,180 @@
+#ifndef EIGEN_COMPLEX_H
+#define EIGEN_COMPLEX_H
+
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// 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 .
+
+// Eigen::Complex reuses as much as possible from std::complex
+// and allows easy conversion to and from, even at the pointer level.
+
+
+#include
+
+namespace Eigen {
+
+template
+struct castable_pointer
+{
+ castable_pointer(_NativePtr ptr) : _ptr(ptr) {}
+ operator _NativePtr () {return _ptr;}
+ operator _PunnedPtr () {return reinterpret_cast<_PunnedPtr>(_ptr);}
+ private:
+ _NativePtr _ptr;
+};
+
+template
+struct Complex
+{
+ typedef typename std::complex StandardComplex;
+ typedef T value_type;
+
+ // constructors
+ Complex(const T& re = T(), const T& im = T()) : _re(re),_im(im) { }
+ Complex(const Complex&other ): _re(other.real()) ,_im(other.imag()) {}
+
+ template
+ Complex(const Complex&other): _re(other.real()) ,_im(other.imag()) {}
+ template
+ Complex(const std::complex&other): _re(other.real()) ,_im(other.imag()) {}
+
+
+ // allow binary access to the object as a std::complex
+ typedef castable_pointer< Complex*, StandardComplex* > pointer_type;
+ typedef castable_pointer< const Complex*, const StandardComplex* > const_pointer_type;
+ pointer_type operator & () {return pointer_type(this);}
+ const_pointer_type operator & () const {return const_pointer_type(this);}
+
+ operator StandardComplex () const {return std_type();}
+ operator StandardComplex & () {return std_type();}
+
+ StandardComplex std_type() const {return StandardComplex(real(),imag());}
+ StandardComplex & std_type() {return *reinterpret_cast(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
+ const T & real() const {return _re;}
+ const T & imag() const {return _im;}
+ T & real() {return _re;}
+ T & imag() {return _im;}
+ T & real(const T & x) {return _re=x;}
+ T & imag(const T & x) {return _im=x;}
+ void set_real(const T & x) {_re = x;}
+ void set_imag(const T & x) {_im = x;}
+
+ // *** complex member functions: ***
+ Complex& operator= (const T& val) { _re=val;_im=0;return *this; }
+ Complex& operator+= (const T& val) {_re+=val;return *this;}
+ Complex& operator-= (const T& val) {_re-=val;return *this;}
+ Complex& operator*= (const T& val) {_re*=val;_im*=val;return *this; }
+ Complex& operator/= (const T& val) {_re/=val;_im/=val;return *this; }
+
+ Complex& operator= (const Complex& rhs) {_re=rhs._re;_im=rhs._im;return *this;}
+ Complex& operator= (const StandardComplex& rhs) {_re=rhs.real();_im=rhs.imag();return *this;}
+
+ template Complex& operator= (const Complex& rhs) { _re=rhs._re;_im=rhs._im;return *this;}
+ template Complex& operator+= (const Complex& rhs) { _re+=rhs._re;_im+=rhs._im;return *this;}
+ template Complex& operator-= (const Complex& rhs) { _re-=rhs._re;_im-=rhs._im;return *this;}
+ template Complex& operator*= (const Complex& rhs) { this->std_type() *= rhs.std_type(); return *this; }
+ template Complex& operator/= (const Complex& rhs) { this->std_type() /= rhs.std_type(); return *this; }
+
+ private:
+ T _re;
+ T _im;
+};
+
+template
+T ei_to_std( const T & x) {return x;}
+
+template
+std::complex ei_to_std( const Complex & x) {return x.std_type();}
+
+// 26.2.6 operators
+template Complex operator+(const Complex& rhs) {return rhs;}
+template Complex operator-(const Complex& rhs) {return -ei_to_std(rhs);}
+
+template Complex operator+(const Complex& lhs, const Complex& rhs) { return ei_to_std(lhs) + ei_to_std(rhs);}
+template Complex operator-(const Complex& lhs, const Complex& rhs) { return ei_to_std(lhs) - ei_to_std(rhs);}
+template Complex operator*(const Complex& lhs, const Complex& rhs) { return ei_to_std(lhs) * ei_to_std(rhs);}
+template Complex operator/(const Complex& lhs, const Complex& rhs) { return ei_to_std(lhs) / ei_to_std(rhs);}
+template bool operator==(const Complex& lhs, const Complex& rhs) { return ei_to_std(lhs) == ei_to_std(rhs);}
+template bool operator!=(const Complex& lhs, const Complex& rhs) { return ei_to_std(lhs) != ei_to_std(rhs);}
+
+template Complex operator+(const Complex& lhs, const T& rhs) {return ei_to_std(lhs) + ei_to_std(rhs); }
+template Complex operator-(const Complex& lhs, const T& rhs) {return ei_to_std(lhs) - ei_to_std(rhs); }
+template Complex operator*(const Complex& lhs, const T& rhs) {return ei_to_std(lhs) * ei_to_std(rhs); }
+template Complex operator/(const Complex& lhs, const T& rhs) {return ei_to_std(lhs) / ei_to_std(rhs); }
+template bool operator==(const Complex& lhs, const T& rhs) {return ei_to_std(lhs) == ei_to_std(rhs); }
+template bool operator!=(const Complex& lhs, const T& rhs) {return ei_to_std(lhs) != ei_to_std(rhs); }
+
+template Complex operator+(const T& lhs, const Complex& rhs) {return ei_to_std(lhs) + ei_to_std(rhs); }
+template Complex operator-(const T& lhs, const Complex& rhs) {return ei_to_std(lhs) - ei_to_std(rhs); }
+template Complex operator*(const T& lhs, const Complex& rhs) {return ei_to_std(lhs) * ei_to_std(rhs); }
+template Complex operator/(const T& lhs, const Complex& rhs) {return ei_to_std(lhs) / ei_to_std(rhs); }
+template bool operator==(const T& lhs, const Complex& rhs) {return ei_to_std(lhs) == ei_to_std(rhs); }
+template bool operator!=(const T& lhs, const Complex& rhs) {return ei_to_std(lhs) != ei_to_std(rhs); }
+
+template
+std::basic_istream&
+ operator>> (std::basic_istream& istr, Complex& rhs)
+{
+ return istr >> rhs.std_type();
+}
+
+template
+std::basic_ostream&
+operator<< (std::basic_ostream& ostr, const Complex& rhs)
+{
+ return ostr << rhs.std_type();
+}
+
+ // 26.2.7 values:
+ template T real(const Complex&x) {return real(ei_to_std(x));}
+ template T abs(const Complex&x) {return abs(ei_to_std(x));}
+ template T arg(const Complex&x) {return arg(ei_to_std(x));}
+ template T norm(const Complex&x) {return norm(ei_to_std(x));}
+
+ template Complex conj(const Complex&x) { return conj(ei_to_std(x));}
+ template Complex polar(const T& x, const T&y) {return polar(ei_to_std(x),ei_to_std(y));}
+ // 26.2.8 transcendentals:
+ template Complex cos (const Complex&x){return cos(ei_to_std(x));}
+ template Complex cosh (const Complex&x){return cosh(ei_to_std(x));}
+ template Complex exp (const Complex&x){return exp(ei_to_std(x));}
+ template Complex log (const Complex&x){return log(ei_to_std(x));}
+ template Complex log10 (const Complex&x){return log10(ei_to_std(x));}
+
+ template Complex pow(const Complex&x, int p) {return pow(ei_to_std(x),ei_to_std(p));}
+ template Complex pow(const Complex&x, const T&p) {return pow(ei_to_std(x),ei_to_std(p));}
+ template Complex pow(const Complex&x, const Complex&p) {return pow(ei_to_std(x),ei_to_std(p));}
+ template Complex pow(const T&x, const Complex&p) {return pow(ei_to_std(x),ei_to_std(p));}
+
+ template Complex sin (const Complex&x){return sin(ei_to_std(x));}
+ template Complex sinh (const Complex&x){return sinh(ei_to_std(x));}
+ template Complex sqrt (const Complex&x){return sqrt(ei_to_std(x));}
+ template Complex tan (const Complex&x){return tan(ei_to_std(x));}
+ template Complex tanh (const Complex&x){return tanh(ei_to_std(x));}
+}
+
+#endif
diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT
new file mode 100644
index 000000000..dc7e85908
--- /dev/null
+++ b/unsupported/Eigen/FFT
@@ -0,0 +1,95 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// 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 .
+
+#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
+
+// 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, commerical -- incompatible with Eigen in GPL form
+#ifdef _MKL_DFTI_H_ // mkl_dfti.h has been included, we can use MKL FFT routines
+// TODO
+// #include "src/FFT/ei_imkl_impl.h"
+// #define DEFAULT_FFT_IMPL ei_imkl_impl
+#endif
+
+namespace Eigen {
+
+template
+ >
+class FFT
+{
+ public:
+ typedef _Traits traits_type;
+ typedef typename traits_type::Scalar Scalar;
+ typedef typename traits_type::Complex Complex;
+
+ FFT(const traits_type & traits=traits_type() ) :m_traits(traits) { }
+
+ template
+ void fwd( Complex * dst, const _Input * src, int nfft)
+ {
+ m_traits.fwd(dst,src,nfft);
+ }
+
+ template
+ void fwd( std::vector & dst, const std::vector<_Input> & src)
+ {
+ dst.resize( src.size() );
+ fwd( &dst[0],&src[0],src.size() );
+ }
+
+ template
+ void inv( _Output * dst, const Complex * src, int nfft)
+ {
+ m_traits.inv( dst,src,nfft );
+ }
+
+ template
+ void inv( std::vector<_Output> & dst, const std::vector & src)
+ {
+ dst.resize( src.size() );
+ inv( &dst[0],&src[0],src.size() );
+ }
+
+ // TODO: multi-dimensional FFTs
+ // TODO: handle Eigen MatrixBase
+
+ traits_type & traits() {return m_traits;}
+ private:
+ traits_type m_traits;
+};
+#undef DEFAULT_FFT_IMPL
+}
+#endif
diff --git a/unsupported/Eigen/src/FFT/ei_fftw_impl.h b/unsupported/Eigen/src/FFT/ei_fftw_impl.h
new file mode 100644
index 000000000..d592bbb20
--- /dev/null
+++ b/unsupported/Eigen/src/FFT/ei_fftw_impl.h
@@ -0,0 +1,198 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// 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 .
+
+namespace Eigen {
+ // FFTW uses non-const arguments
+ // so we must use ugly const_cast calls for all the args it uses
+ //
+ // This should be safe as long as
+ // 1. we use FFTW_ESTIMATE for all our planning
+ // see the FFTW docs section 4.3.2 "Planner Flags"
+ // 2. fftw_complex is compatible with std::complex
+ // This assumes std::complex layout is array of size 2 with real,imag
+ template
+ T * ei_fftw_cast(const T* p)
+ {
+ return const_cast( p);
+ }
+
+ fftw_complex * ei_fftw_cast( const std::complex * p)
+ {
+ return const_cast( reinterpret_cast(p) );
+ }
+
+ fftwf_complex * ei_fftw_cast( const std::complex * p)
+ {
+ return const_cast( reinterpret_cast(p) );
+ }
+
+ fftwl_complex * ei_fftw_cast( const std::complex * p)
+ {
+ return const_cast( reinterpret_cast(p) );
+ }
+
+ template
+ struct ei_fftw_plan {};
+
+ template <>
+ struct ei_fftw_plan
+ {
+ typedef float scalar_type;
+ typedef fftwf_complex complex_type;
+ fftwf_plan m_plan;
+ ei_fftw_plan() :m_plan(NULL) {}
+ ~ei_fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
+
+ 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);
+ }
+ 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);
+ }
+ 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);
+ }
+ 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);
+ fftwf_execute_dft_c2r( m_plan, src,dst);
+ }
+ };
+ template <>
+ struct ei_fftw_plan
+ {
+ typedef double scalar_type;
+ typedef fftw_complex complex_type;
+ fftw_plan m_plan;
+ ei_fftw_plan() :m_plan(NULL) {}
+ ~ei_fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
+
+ 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);
+ }
+ 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);
+ }
+ 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);
+ }
+ 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);
+ fftw_execute_dft_c2r( m_plan, src,dst);
+ }
+ };
+ template <>
+ struct ei_fftw_plan
+ {
+ typedef long double scalar_type;
+ typedef fftwl_complex complex_type;
+ fftwl_plan m_plan;
+ ei_fftw_plan() :m_plan(NULL) {}
+ ~ei_fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
+
+ 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);
+ }
+ 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);
+ }
+ 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);
+ }
+ 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);
+ fftwl_execute_dft_c2r( m_plan, src,dst);
+ }
+ };
+
+ template
+ struct ei_fftw_impl
+ {
+ typedef _Scalar Scalar;
+ typedef std::complex Complex;
+
+ void clear()
+ {
+ m_plans.clear();
+ }
+
+ 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
+ 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
+ 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 );
+ // scaling
+ Scalar s = 1./nfft;
+ for (int k=0;k PlanData;
+ typedef std::map PlanMap;
+
+ PlanMap m_plans;
+
+ PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
+ {
+ bool inplace = (dst==src);
+ bool aligned = ( (reinterpret_cast(src)&15) | (reinterpret_cast(dst)&15) ) == 0;
+ int key = (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned;
+ return m_plans[key];
+ }
+ };
+}
diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
new file mode 100644
index 000000000..a84ac68a0
--- /dev/null
+++ b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
@@ -0,0 +1,412 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// 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 .
+
+#include
+#include
+#include