Move into version control

1 files changed, 419 insertions, 0 deletions

diff --git a/src/EigenUnsupported/FFT b/src/EigenUnsupported/FFT
new file mode 100644
index 0000000..c8c311a
--- /dev/null
+++ b/src/EigenUnsupported/FFT
@@ -0,0 +1,419 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. 
+//
+// Copyright (C) 2009 Mark Borgerding mark a borgerding net
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_FFT_H
+#define EIGEN_FFT_H
+
+#include <complex>
+#include <vector>
+#include <map>
+#include "../../Eigen/Core"
+
+
+/**
+  * \defgroup FFT_Module Fast Fourier Transform module
+  *
+  * \code
+  * #include <unsupported/Eigen/FFT>
+  * \endcode
+  *
+  * This module provides Fast Fourier transformation, with a configurable backend
+  * implementation.
+  *
+  * The default implementation is based on kissfft. It is a small, free, and
+  * reasonably efficient default.
+  *
+  * There are currently two implementation backend:
+  *
+  * - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size.
+  * - MKL (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form.
+  *
+  * \section FFTDesign Design
+  *
+  * The following design decisions were made concerning scaling and
+  * half-spectrum for real FFT.
+  *
+  * The intent is to facilitate generic programming and ease migrating code
+  * from  Matlab/octave.
+  * We think the default behavior of Eigen/FFT should favor correctness and
+  * generality over speed. Of course, the caller should be able to "opt-out" from this
+  * behavior and get the speed increase if they want it.
+  *
+  * 1) %Scaling:
+  * Other libraries (FFTW,IMKL,KISSFFT)  do not perform scaling, so there
+  * is a constant gain incurred after the forward&inverse transforms , so 
+  * IFFT(FFT(x)) = Kx;  this is done to avoid a vector-by-value multiply.  
+  * The downside is that algorithms that worked correctly in Matlab/octave 
+  * don't behave the same way once implemented in C++.
+  *
+  * How Eigen/FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x. 
+  *
+  * 2) Real FFT half-spectrum
+  * Other libraries use only half the frequency spectrum (plus one extra 
+  * sample for the Nyquist bin) for a real FFT, the other half is the 
+  * conjugate-symmetric of the first half.  This saves them a copy and some 
+  * memory.  The downside is the caller needs to have special logic for the 
+  * number of bins in complex vs real.
+  *
+  * How Eigen/FFT differs: The full spectrum is returned from the forward 
+  * transform.  This facilitates generic template programming by obviating 
+  * separate specializations for real vs complex.  On the inverse
+  * transform, only half the spectrum is actually used if the output type is real.
+  */
+ 
+
+#include "../../Eigen/src/Core/util/DisableStupidWarnings.h"
+
+#ifdef EIGEN_FFTW_DEFAULT
+// FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
+#  include <fftw3.h>
+#  include "src/FFT/ei_fftw_impl.h"
+   namespace Eigen {
+     //template <typename T> typedef struct internal::fftw_impl  default_fft_impl; this does not work
+     template <typename T> struct default_fft_impl : public internal::fftw_impl<T> {};
+   }
+#elif defined EIGEN_MKL_DEFAULT
+// TODO 
+// intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form
+#  include "src/FFT/ei_imklfft_impl.h"
+   namespace Eigen {
+     template <typename T> struct default_fft_impl : public internal::imklfft_impl {};
+   }
+#else
+// internal::kissfft_impl:  small, free, reasonably efficient default, derived from kissfft
+//
+# include "src/FFT/ei_kissfft_impl.h"
+  namespace Eigen {
+     template <typename T> 
+       struct default_fft_impl : public internal::kissfft_impl<T> {};
+  }
+#endif
+
+namespace Eigen {
+
+ 
+// 
+template<typename T_SrcMat,typename T_FftIfc> struct fft_fwd_proxy;
+template<typename T_SrcMat,typename T_FftIfc> struct fft_inv_proxy;
+
+namespace internal {
+template<typename T_SrcMat,typename T_FftIfc>
+struct traits< fft_fwd_proxy<T_SrcMat,T_FftIfc> >
+{
+  typedef typename T_SrcMat::PlainObject ReturnType;
+};
+template<typename T_SrcMat,typename T_FftIfc>
+struct traits< fft_inv_proxy<T_SrcMat,T_FftIfc> >
+{
+  typedef typename T_SrcMat::PlainObject ReturnType;
+};
+}
+
+template<typename T_SrcMat,typename T_FftIfc> 
+struct fft_fwd_proxy
+ : public ReturnByValue<fft_fwd_proxy<T_SrcMat,T_FftIfc> >
+{
+  typedef DenseIndex Index;
+
+  fft_fwd_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {}
+
+  template<typename T_DestMat> void evalTo(T_DestMat& dst) const;
+
+  Index rows() const { return m_src.rows(); }
+  Index cols() const { return m_src.cols(); }
+protected:
+  const T_SrcMat & m_src;
+  T_FftIfc & m_ifc;
+  Index m_nfft;
+};
+
+template<typename T_SrcMat,typename T_FftIfc> 
+struct fft_inv_proxy
+ : public ReturnByValue<fft_inv_proxy<T_SrcMat,T_FftIfc> >
+{
+  typedef DenseIndex Index;
+
+  fft_inv_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {}
+
+  template<typename T_DestMat> void evalTo(T_DestMat& dst) const;
+
+  Index rows() const { return m_src.rows(); }
+  Index cols() const { return m_src.cols(); }
+protected:
+  const T_SrcMat & m_src;
+  T_FftIfc & m_ifc;
+  Index m_nfft;
+};
+
+
+template <typename T_Scalar,
+         typename T_Impl=default_fft_impl<T_Scalar> >
+class FFT
+{
+  public:
+    typedef T_Impl impl_type;
+    typedef DenseIndex Index;
+    typedef typename impl_type::Scalar Scalar;
+    typedef typename impl_type::Complex Complex;
+
+    enum Flag {
+      Default=0, // goof proof
+      Unscaled=1,
+      HalfSpectrum=2,
+      // SomeOtherSpeedOptimization=4
+      Speedy=32767
+    };
+
+    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, Index nfft)
+    {
+        m_impl.fwd(dst,src,static_cast<int>(nfft));
+        if ( HasFlag(HalfSpectrum) == false)
+          ReflectSpectrum(dst,nfft);
+    }
+
+    inline
+    void fwd( Complex * dst, const Complex * src, Index nfft)
+    {
+        m_impl.fwd(dst,src,static_cast<int>(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); // half the bins + Nyquist bin
+      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, Index nfft=-1)
+    {
+      typedef typename ComplexDerived::Scalar dst_type;
+      typedef typename InputDerived::Scalar src_type;
+      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((internal::is_same<dst_type, Complex>::value),
+            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 (nfft<1)
+        nfft = src.size();
+
+      if ( NumTraits< src_type >::IsComplex == 0 && HasFlag(HalfSpectrum) )
+        dst.derived().resize( (nfft>>1)+1);
+      else
+        dst.derived().resize(nfft);
+
+      if ( src.innerStride() != 1 || src.size() < nfft ) {
+        Matrix<src_type,1,Dynamic> tmp;
+        if (src.size()<nfft) {
+          tmp.setZero(nfft);
+          tmp.block(0,0,src.size(),1 ) = src;
+        }else{
+          tmp = src;
+        }
+        fwd( &dst[0],&tmp[0],nfft );
+      }else{
+        fwd( &dst[0],&src[0],nfft );
+      }
+    }
+ 
+    template<typename InputDerived>
+    inline
+    fft_fwd_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
+    fwd( const MatrixBase<InputDerived> & src, Index nfft=-1)
+    {
+      return fft_fwd_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
+    }
+
+    template<typename InputDerived>
+    inline
+    fft_inv_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
+    inv( const MatrixBase<InputDerived> & src, Index nfft=-1)
+    {
+      return  fft_inv_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
+    }
+
+    inline
+    void inv( Complex * dst, const Complex * src, Index nfft)
+    {
+      m_impl.inv( dst,src,static_cast<int>(nfft) );
+      if ( HasFlag( Unscaled ) == false)
+        scale(dst,Scalar(1./nfft),nfft); // scale the time series
+    }
+
+    inline
+    void inv( Scalar * dst, const Complex * src, Index nfft)
+    {
+      m_impl.inv( dst,src,static_cast<int>(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, Index nfft=-1)
+    {
+      typedef typename ComplexDerived::Scalar src_type;
+      typedef typename ComplexDerived::RealScalar real_type;
+      typedef typename OutputDerived::Scalar dst_type;
+      const bool realfft= (NumTraits<dst_type>::IsComplex == 0);
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived)
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived)
+      EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,OutputDerived) // size at compile-time
+      EIGEN_STATIC_ASSERT((internal::is_same<src_type, Complex>::value),
+            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)
+
+      if (nfft<1) { //automatic FFT size determination
+        if ( realfft && HasFlag(HalfSpectrum) ) 
+          nfft = 2*(src.size()-1); //assume even fft size
+        else
+          nfft = src.size();
+      }
+      dst.derived().resize( nfft );
+
+      // check for nfft that does not fit the input data size
+      Index resize_input= ( realfft && HasFlag(HalfSpectrum) )
+        ? ( (nfft/2+1) - src.size() )
+        : ( nfft - src.size() );
+
+      if ( src.innerStride() != 1 || resize_input ) {
+        // if the vector is strided, then we need to copy it to a packed temporary
+        Matrix<src_type,1,Dynamic> tmp;
+        if ( resize_input ) {
+          size_t ncopy = (std::min)(src.size(),src.size() + resize_input);
+          tmp.setZero(src.size() + resize_input);
+          if ( realfft && HasFlag(HalfSpectrum) ) {
+            // pad at the Nyquist bin
+            tmp.head(ncopy) = src.head(ncopy);
+            tmp(ncopy-1) = real(tmp(ncopy-1)); // enforce real-only Nyquist bin
+          }else{
+            size_t nhead,ntail;
+            nhead = 1+ncopy/2-1; // range  [0:pi)
+            ntail = ncopy/2-1;   // range (-pi:0)
+            tmp.head(nhead) = src.head(nhead);
+            tmp.tail(ntail) = src.tail(ntail);
+            if (resize_input<0) { //shrinking -- create the Nyquist bin as the average of the two bins that fold into it
+              tmp(nhead) = ( src(nfft/2) + src( src.size() - nfft/2 ) )*real_type(.5);
+            }else{ // expanding -- split the old Nyquist bin into two halves
+              tmp(nhead) = src(nhead) * real_type(.5);
+              tmp(tmp.size()-nhead) = tmp(nhead);
+            }
+          }
+        }else{
+          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,Index nfft=-1)
+    {
+      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
+    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 T_Data>
+    inline
+    void scale(T_Data * x,Scalar s,Index 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
+    void ReflectSpectrum(Complex * freq, Index nfft)
+    {
+      // create the implicit right-half spectrum (conjugate-mirror of the left-half)
+      Index nhbins=(nfft>>1)+1;
+      for (Index k=nhbins;k < nfft; ++k )
+        freq[k] = conj(freq[nfft-k]);
+    }
+
+    impl_type m_impl;
+    int m_flag;
+};
+
+template<typename T_SrcMat,typename T_FftIfc> 
+template<typename T_DestMat> inline 
+void fft_fwd_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const
+{
+    m_ifc.fwd( dst, m_src, m_nfft);
+}
+
+template<typename T_SrcMat,typename T_FftIfc> 
+template<typename T_DestMat> inline 
+void fft_inv_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const
+{
+    m_ifc.inv( dst, m_src, m_nfft);
+}
+
+}
+
+#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
+
+#endif