1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra. 
3//
4// Copyright (C) 2009 Mark Borgerding mark a borgerding net
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_FFT_H
11#define EIGEN_FFT_H
12
13#include <complex>
14#include <vector>
15#include <map>
16#include <Eigen/Core>
17
18
19/**
20  * \defgroup FFT_Module Fast Fourier Transform module
21  *
22  * \code
23  * #include <unsupported/Eigen/FFT>
24  * \endcode
25  *
26  * This module provides Fast Fourier transformation, with a configurable backend
27  * implementation.
28  *
29  * The default implementation is based on kissfft. It is a small, free, and
30  * reasonably efficient default.
31  *
32  * There are currently two implementation backend:
33  *
34  * - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size.
35  * - MKL (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form.
36  *
37  * \section FFTDesign Design
38  *
39  * The following design decisions were made concerning scaling and
40  * half-spectrum for real FFT.
41  *
42  * The intent is to facilitate generic programming and ease migrating code
43  * from  Matlab/octave.
44  * We think the default behavior of Eigen/FFT should favor correctness and
45  * generality over speed. Of course, the caller should be able to "opt-out" from this
46  * behavior and get the speed increase if they want it.
47  *
48  * 1) %Scaling:
49  * Other libraries (FFTW,IMKL,KISSFFT)  do not perform scaling, so there
50  * is a constant gain incurred after the forward&inverse transforms , so 
51  * IFFT(FFT(x)) = Kx;  this is done to avoid a vector-by-value multiply.  
52  * The downside is that algorithms that worked correctly in Matlab/octave 
53  * don't behave the same way once implemented in C++.
54  *
55  * How Eigen/FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x. 
56  *
57  * 2) Real FFT half-spectrum
58  * Other libraries use only half the frequency spectrum (plus one extra 
59  * sample for the Nyquist bin) for a real FFT, the other half is the 
60  * conjugate-symmetric of the first half.  This saves them a copy and some 
61  * memory.  The downside is the caller needs to have special logic for the 
62  * number of bins in complex vs real.
63  *
64  * How Eigen/FFT differs: The full spectrum is returned from the forward 
65  * transform.  This facilitates generic template programming by obviating 
66  * separate specializations for real vs complex.  On the inverse
67  * transform, only half the spectrum is actually used if the output type is real.
68  */
69 
70
71#ifdef EIGEN_FFTW_DEFAULT
72// FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
73#  include <fftw3.h>
74#  include "src/FFT/ei_fftw_impl.h"
75   namespace Eigen {
76     //template <typename T> typedef struct internal::fftw_impl  default_fft_impl; this does not work
77     template <typename T> struct default_fft_impl : public internal::fftw_impl<T> {};
78   }
79#elif defined EIGEN_MKL_DEFAULT
80// TODO 
81// intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form
82#  include "src/FFT/ei_imklfft_impl.h"
83   namespace Eigen {
84     template <typename T> struct default_fft_impl : public internal::imklfft_impl {};
85   }
86#else
87// internal::kissfft_impl:  small, free, reasonably efficient default, derived from kissfft
88//
89# include "src/FFT/ei_kissfft_impl.h"
90  namespace Eigen {
91     template <typename T> 
92       struct default_fft_impl : public internal::kissfft_impl<T> {};
93  }
94#endif
95
96namespace Eigen {
97
98 
99// 
100template<typename T_SrcMat,typename T_FftIfc> struct fft_fwd_proxy;
101template<typename T_SrcMat,typename T_FftIfc> struct fft_inv_proxy;
102
103namespace internal {
104template<typename T_SrcMat,typename T_FftIfc>
105struct traits< fft_fwd_proxy<T_SrcMat,T_FftIfc> >
106{
107  typedef typename T_SrcMat::PlainObject ReturnType;
108};
109template<typename T_SrcMat,typename T_FftIfc>
110struct traits< fft_inv_proxy<T_SrcMat,T_FftIfc> >
111{
112  typedef typename T_SrcMat::PlainObject ReturnType;
113};
114}
115
116template<typename T_SrcMat,typename T_FftIfc> 
117struct fft_fwd_proxy
118 : public ReturnByValue<fft_fwd_proxy<T_SrcMat,T_FftIfc> >
119{
120  typedef DenseIndex Index;
121
122  fft_fwd_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {}
123
124  template<typename T_DestMat> void evalTo(T_DestMat& dst) const;
125
126  Index rows() const { return m_src.rows(); }
127  Index cols() const { return m_src.cols(); }
128protected:
129  const T_SrcMat & m_src;
130  T_FftIfc & m_ifc;
131  Index m_nfft;
132private:
133  fft_fwd_proxy& operator=(const fft_fwd_proxy&);
134};
135
136template<typename T_SrcMat,typename T_FftIfc> 
137struct fft_inv_proxy
138 : public ReturnByValue<fft_inv_proxy<T_SrcMat,T_FftIfc> >
139{
140  typedef DenseIndex Index;
141
142  fft_inv_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {}
143
144  template<typename T_DestMat> void evalTo(T_DestMat& dst) const;
145
146  Index rows() const { return m_src.rows(); }
147  Index cols() const { return m_src.cols(); }
148protected:
149  const T_SrcMat & m_src;
150  T_FftIfc & m_ifc;
151  Index m_nfft;
152private:
153  fft_inv_proxy& operator=(const fft_inv_proxy&);
154};
155
156
157template <typename T_Scalar,
158         typename T_Impl=default_fft_impl<T_Scalar> >
159class FFT
160{
161  public:
162    typedef T_Impl impl_type;
163    typedef DenseIndex Index;
164    typedef typename impl_type::Scalar Scalar;
165    typedef typename impl_type::Complex Complex;
166
167    enum Flag {
168      Default=0, // goof proof
169      Unscaled=1,
170      HalfSpectrum=2,
171      // SomeOtherSpeedOptimization=4
172      Speedy=32767
173    };
174
175    FFT( const impl_type & impl=impl_type() , Flag flags=Default ) :m_impl(impl),m_flag(flags) { }
176
177    inline
178    bool HasFlag(Flag f) const { return (m_flag & (int)f) == f;}
179
180    inline
181    void SetFlag(Flag f) { m_flag |= (int)f;}
182
183    inline
184    void ClearFlag(Flag f) { m_flag &= (~(int)f);}
185
186    inline
187    void fwd( Complex * dst, const Scalar * src, Index nfft)
188    {
189        m_impl.fwd(dst,src,static_cast<int>(nfft));
190        if ( HasFlag(HalfSpectrum) == false)
191          ReflectSpectrum(dst,nfft);
192    }
193
194    inline
195    void fwd( Complex * dst, const Complex * src, Index nfft)
196    {
197        m_impl.fwd(dst,src,static_cast<int>(nfft));
198    }
199
200    /*
201    inline 
202    void fwd2(Complex * dst, const Complex * src, int n0,int n1)
203    {
204      m_impl.fwd2(dst,src,n0,n1);
205    }
206    */
207
208    template <typename _Input>
209    inline
210    void fwd( std::vector<Complex> & dst, const std::vector<_Input> & src) 
211    {
212      if ( NumTraits<_Input>::IsComplex == 0 && HasFlag(HalfSpectrum) )
213        dst.resize( (src.size()>>1)+1); // half the bins + Nyquist bin
214      else
215        dst.resize(src.size());
216      fwd(&dst[0],&src[0],src.size());
217    }
218
219    template<typename InputDerived, typename ComplexDerived>
220    inline
221    void fwd( MatrixBase<ComplexDerived> & dst, const MatrixBase<InputDerived> & src, Index nfft=-1)
222    {
223      typedef typename ComplexDerived::Scalar dst_type;
224      typedef typename InputDerived::Scalar src_type;
225      EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived)
226      EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived)
227      EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,InputDerived) // size at compile-time
228      EIGEN_STATIC_ASSERT((internal::is_same<dst_type, Complex>::value),
229            YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
230      EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit,
231            THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
232
233      if (nfft<1)
234        nfft = src.size();
235
236      if ( NumTraits< src_type >::IsComplex == 0 && HasFlag(HalfSpectrum) )
237        dst.derived().resize( (nfft>>1)+1);
238      else
239        dst.derived().resize(nfft);
240
241      if ( src.innerStride() != 1 || src.size() < nfft ) {
242        Matrix<src_type,1,Dynamic> tmp;
243        if (src.size()<nfft) {
244          tmp.setZero(nfft);
245          tmp.block(0,0,src.size(),1 ) = src;
246        }else{
247          tmp = src;
248        }
249        fwd( &dst[0],&tmp[0],nfft );
250      }else{
251        fwd( &dst[0],&src[0],nfft );
252      }
253    }
254 
255    template<typename InputDerived>
256    inline
257    fft_fwd_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
258    fwd( const MatrixBase<InputDerived> & src, Index nfft=-1)
259    {
260      return fft_fwd_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
261    }
262
263    template<typename InputDerived>
264    inline
265    fft_inv_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
266    inv( const MatrixBase<InputDerived> & src, Index nfft=-1)
267    {
268      return  fft_inv_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
269    }
270
271    inline
272    void inv( Complex * dst, const Complex * src, Index nfft)
273    {
274      m_impl.inv( dst,src,static_cast<int>(nfft) );
275      if ( HasFlag( Unscaled ) == false)
276        scale(dst,Scalar(1./nfft),nfft); // scale the time series
277    }
278
279    inline
280    void inv( Scalar * dst, const Complex * src, Index nfft)
281    {
282      m_impl.inv( dst,src,static_cast<int>(nfft) );
283      if ( HasFlag( Unscaled ) == false)
284        scale(dst,Scalar(1./nfft),nfft); // scale the time series
285    }
286
287    template<typename OutputDerived, typename ComplexDerived>
288    inline
289    void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src, Index nfft=-1)
290    {
291      typedef typename ComplexDerived::Scalar src_type;
292      typedef typename OutputDerived::Scalar dst_type;
293      const bool realfft= (NumTraits<dst_type>::IsComplex == 0);
294      EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived)
295      EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived)
296      EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,OutputDerived) // size at compile-time
297      EIGEN_STATIC_ASSERT((internal::is_same<src_type, Complex>::value),
298            YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
299      EIGEN_STATIC_ASSERT(int(OutputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit,
300            THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
301
302      if (nfft<1) { //automatic FFT size determination
303        if ( realfft && HasFlag(HalfSpectrum) ) 
304          nfft = 2*(src.size()-1); //assume even fft size
305        else
306          nfft = src.size();
307      }
308      dst.derived().resize( nfft );
309
310      // check for nfft that does not fit the input data size
311      Index resize_input= ( realfft && HasFlag(HalfSpectrum) )
312        ? ( (nfft/2+1) - src.size() )
313        : ( nfft - src.size() );
314
315      if ( src.innerStride() != 1 || resize_input ) {
316        // if the vector is strided, then we need to copy it to a packed temporary
317        Matrix<src_type,1,Dynamic> tmp;
318        if ( resize_input ) {
319          size_t ncopy = (std::min)(src.size(),src.size() + resize_input);
320          tmp.setZero(src.size() + resize_input);
321          if ( realfft && HasFlag(HalfSpectrum) ) {
322            // pad at the Nyquist bin
323            tmp.head(ncopy) = src.head(ncopy);
324            tmp(ncopy-1) = real(tmp(ncopy-1)); // enforce real-only Nyquist bin
325          }else{
326            size_t nhead,ntail;
327            nhead = 1+ncopy/2-1; // range  [0:pi)
328            ntail = ncopy/2-1;   // range (-pi:0)
329            tmp.head(nhead) = src.head(nhead);
330            tmp.tail(ntail) = src.tail(ntail);
331            if (resize_input<0) { //shrinking -- create the Nyquist bin as the average of the two bins that fold into it
332              tmp(nhead) = ( src(nfft/2) + src( src.size() - nfft/2 ) )*src_type(.5);
333            }else{ // expanding -- split the old Nyquist bin into two halves
334              tmp(nhead) = src(nhead) * src_type(.5);
335              tmp(tmp.size()-nhead) = tmp(nhead);
336            }
337          }
338        }else{
339          tmp = src;
340        }
341        inv( &dst[0],&tmp[0], nfft);
342      }else{
343        inv( &dst[0],&src[0], nfft);
344      }
345    }
346
347    template <typename _Output>
348    inline
349    void inv( std::vector<_Output> & dst, const std::vector<Complex> & src,Index nfft=-1)
350    {
351      if (nfft<1)
352        nfft = ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) ? 2*(src.size()-1) : src.size();
353      dst.resize( nfft );
354      inv( &dst[0],&src[0],nfft);
355    }
356
357
358    /*
359    // TODO: multi-dimensional FFTs
360    inline 
361    void inv2(Complex * dst, const Complex * src, int n0,int n1)
362    {
363      m_impl.inv2(dst,src,n0,n1);
364      if ( HasFlag( Unscaled ) == false)
365          scale(dst,1./(n0*n1),n0*n1);
366    }
367  */
368
369    inline
370    impl_type & impl() {return m_impl;}
371  private:
372
373    template <typename T_Data>
374    inline
375    void scale(T_Data * x,Scalar s,Index nx)
376    {
377#if 1
378      for (int k=0;k<nx;++k)
379        *x++ *= s;
380#else
381      if ( ((ptrdiff_t)x) & 15 )
382        Matrix<T_Data, Dynamic, 1>::Map(x,nx) *= s;
383      else
384        Matrix<T_Data, Dynamic, 1>::MapAligned(x,nx) *= s;
385         //Matrix<T_Data, Dynamic, Dynamic>::Map(x,nx) * s;
386#endif  
387    }
388
389    inline
390    void ReflectSpectrum(Complex * freq, Index nfft)
391    {
392      // create the implicit right-half spectrum (conjugate-mirror of the left-half)
393      Index nhbins=(nfft>>1)+1;
394      for (Index k=nhbins;k < nfft; ++k )
395        freq[k] = conj(freq[nfft-k]);
396    }
397
398    impl_type m_impl;
399    int m_flag;
400};
401
402template<typename T_SrcMat,typename T_FftIfc> 
403template<typename T_DestMat> inline 
404void fft_fwd_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const
405{
406    m_ifc.fwd( dst, m_src, m_nfft);
407}
408
409template<typename T_SrcMat,typename T_FftIfc> 
410template<typename T_DestMat> inline 
411void fft_inv_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const
412{
413    m_ifc.inv( dst, m_src, m_nfft);
414}
415
416}
417#endif
418/* vim: set filetype=cpp et sw=2 ts=2 ai: */
419