FFT

// 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"
 

 
 
#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;
 
    typedef int Flag;
    static const Flag Default = 0;
    static const Flag Unscaled = 1;
    static const Flag HalfSpectrum = 2;
    static const Flag Speedy = 32767;
 
    FFT(const impl_type& impl = impl_type(), Flag flags = Default) : m_impl(impl), m_flag(flags)
    {
        eigen_assert((flags == Default || flags == Unscaled || flags == HalfSpectrum || flags == Speedy) && "invalid flags argument");
    }
 
    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