initial pass of FFT module -- includes complex 1-d case only

4 files changed, 466 insertions, 0 deletions

diff --git a/unsupported/Eigen/FFT.h b/unsupported/Eigen/FFT.h
new file mode 100644
index 000000000..03490d2c5
--- /dev/null
+++ b/unsupported/Eigen/FFT.h
@@ -0,0 +1,84 @@
+// 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/>.
+
+#ifndef EIGEN_FFT_H
+#define EIGEN_FFT_H
+
+// simple_fft_traits:  small, free, reasonably efficient default, derived from kissfft
+#include "src/FFT/simple_fft_traits.h"
+#define DEFAULT_FFT_TRAITS simple_fft_traits
+
+// FFTW: faster, GPL-not LGPL, bigger code size
+#ifdef FFTW_PATIENT  // definition of FFTW_PATIENT indicates the caller has included fftw3.h, we can use FFTW routines
+// TODO 
+// #include "src/FFT/fftw_traits.h"
+// #define DEFAULT_FFT_TRAITS fftw_traits
+#endif
+
+// intel Math Kernel Library: fastest, commerical
+#ifdef _MKL_DFTI_H_ // mkl_dfti.h has been included, we can use MKL FFT routines
+// TODO 
+// #include "src/FFT/imkl_traits.h"
+// #define DEFAULT_FFT_TRAITS imkl_traits
+#endif
+
+namespace Eigen {
+
+template <typename _Scalar,
+         typename _Traits=DEFAULT_FFT_TRAITS<_Scalar> 
+         >
+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) { }
+
+    void fwd( Complex * dst, const Complex * src, int nfft)
+    {
+      m_traits.prepare(nfft,false,dst,src);
+      m_traits.exec(dst,src);
+      m_traits.postprocess(dst);
+    }
+
+    void inv( Complex * dst, const Complex * src, int nfft)
+    {
+      m_traits.prepare(nfft,true,dst,src);
+      m_traits.exec(dst,src);
+      m_traits.postprocess(dst);
+    }
+
+    // TODO: fwd,inv for Scalar
+    // TODO: multi-dimensional FFTs
+    // TODO: handle Eigen MatrixBase
+
+    traits_type & traits() {return m_traits;}
+  private:
+    traits_type m_traits;
+};
+#undef DEFAULT_FFT_TRAITS
+}
+#endif
diff --git a/unsupported/Eigen/src/FFT/simple_fft_traits.h b/unsupported/Eigen/src/FFT/simple_fft_traits.h
new file mode 100644
index 000000000..fe9d24b84
--- /dev/null
+++ b/unsupported/Eigen/src/FFT/simple_fft_traits.h
@@ -0,0 +1,296 @@
+// 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/>.
+
+#include <complex>
+#include <vector>
+
+namespace Eigen {
+
+  template <typename _Scalar>
+  struct simple_fft_traits
+  {
+    typedef _Scalar Scalar;
+    typedef std::complex<Scalar> Complex;
+    simple_fft_traits() : m_nfft(0) {} 
+
+    void prepare(int nfft,bool inverse,Complex * dst,const Complex *src)
+    {
+      if (m_nfft == nfft) {
+        // reuse the twiddles, conjugate if necessary
+        if (inverse != m_inverse)
+          for (int i=0;i<nfft;++i)
+            m_twiddles[i] = conj( m_twiddles[i] );
+        m_inverse = inverse;
+        return;
+      }
+      m_nfft = nfft;
+      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_stageRadix.resize(0);
+      m_stageRemainder.resize(0);
+      //factorize
+      //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;// no more factors
+        }
+        n /= p;
+        m_stageRadix.push_back(p);
+        m_stageRemainder.push_back(n);
+      }while(n>1);
+    }
+
+    void exec(Complex * dst, const Complex * src)
+    {
+      work(0, dst, src, 1,1);
+    }
+
+    void postprocess(Complex *dst) 
+    {
+        if (m_inverse) {
+            Scalar scale = 1./m_nfft;
+            for (int k=0;k<m_nfft;++k)
+                dst[k] *= scale;
+        }
+    }
+  private:
+    void work( int stage,Complex * Fout, const Complex * f, size_t fstride,size_t in_stride)
+    {
+      int p = m_stageRadix[stage];
+      int m = m_stageRemainder[stage];
+      Complex * Fout_beg = Fout;
+      Complex * Fout_end = Fout + p*m;
+
+      if (m==1) {
+        do{
+          *Fout = *f;
+          f += fstride*in_stride;
+        }while(++Fout != Fout_end );
+      }else{
+        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, Fout , f, fstride*p,in_stride);
+          f += fstride*in_stride;
+        }while( (Fout += m) != Fout_end );
+      }
+
+      Fout=Fout_beg;
+
+      // recombine the p smaller DFTs 
+      switch (p) {
+        case 2: bfly2(Fout,fstride,m); break;
+        case 3: bfly3(Fout,fstride,m); break;
+        case 4: bfly4(Fout,fstride,m); break;
+        case 5: bfly5(Fout,fstride,m); break;
+        default: bfly_generic(Fout,fstride,m,p); break;
+      }
+    }
+
+    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;
+      }
+    }
+
+    void bfly4( Complex * Fout, const size_t fstride, const size_t m)
+    {
+      Complex scratch[7];
+      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];
+      }
+    }
+
+    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);
+    }
+
+    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())
+              );
+
+        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];
+
+        ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
+      }
+    }
+
+    /* perform the butterfly for one stage of a mixed radix FFT */
+    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_nfft;
+      Complex * scratchbuf = (Complex*)alloca(p*sizeof(Complex) );
+
+      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;
+        }
+      }
+    }
+
+    int m_nfft;
+    bool m_inverse;
+    std::vector<Complex> m_twiddles;
+    std::vector<int> m_stageRadix;
+    std::vector<int> m_stageRemainder;
+  };
+}
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index abfbb0185..49de7dfb2 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -19,3 +19,4 @@ ei_add_test(autodiff)
 ei_add_test(BVH)
 ei_add_test(matrixExponential)
 ei_add_test(alignedvector3)
+ei_add_test(FFT)
diff --git a/unsupported/test/FFT.cpp b/unsupported/test/FFT.cpp
new file mode 100644
index 000000000..b13a7810f
--- /dev/null
+++ b/unsupported/test/FFT.cpp
@@ -0,0 +1,85 @@
+// 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/>.
+
+#include "main.h"
+#include <unsupported/Eigen/FFT.h>
+
+//#include <iostream>
+//#include <cstdlib>
+//#include <typeinfo>
+
+using namespace std;
+
+template <class T>
+void test_fft(int nfft)
+{
+    typedef typename Eigen::FFT<T>::Complex Complex;
+
+    //cout << "type:" << typeid(T).name() << " nfft:" << nfft;
+
+    FFT<T> fft;
+
+    vector<Complex> inbuf(nfft);
+    vector<Complex> buf3(nfft);
+    vector<Complex> outbuf(nfft);
+    for (int k=0;k<nfft;++k)
+        inbuf[k]= Complex( 
+                (T)(rand()/(double)RAND_MAX - .5),
+                (T)(rand()/(double)RAND_MAX - .5) );
+    fft.fwd( &outbuf[0] , &inbuf[0] ,nfft);
+    fft.inv( &buf3[0] , &outbuf[0] ,nfft);
+
+    long double totalpower=0;
+    long double difpower=0;
+    for (int k0=0;k0<nfft;++k0) {
+        complex<long double> acc = 0;
+        long double phinc = 2*k0* M_PIl / nfft;
+        for (int k1=0;k1<nfft;++k1) {
+            complex<long double> x(inbuf[k1].real(),inbuf[k1].imag()); 
+            acc += x * exp( complex<long double>(0,-k1*phinc) );
+        }
+        totalpower += norm(acc);
+        complex<long double> x(outbuf[k0].real(),outbuf[k0].imag()); 
+        complex<long double> dif = acc - x;
+        difpower += norm(dif);
+    }
+    long double rmse = sqrt(difpower/totalpower);
+    VERIFY( rmse < 1e-5 );// gross check
+
+    totalpower=0;
+    difpower=0;
+    for (int k=0;k<nfft;++k) {
+        totalpower += norm( inbuf[k] );
+        difpower += norm(inbuf[k] - buf3[k]);
+    }
+    rmse = sqrt(difpower/totalpower);
+    VERIFY( rmse < 1e-5 );// gross check
+}
+
+void test_FFT()
+{
+  CALL_SUBTEST(( test_fft<float>(32) )); CALL_SUBTEST(( test_fft<double>(32) )); CALL_SUBTEST(( test_fft<long double>(32) ));
+  CALL_SUBTEST(( test_fft<float>(1024) )); CALL_SUBTEST(( test_fft<double>(1024) )); CALL_SUBTEST(( test_fft<long double>(1024) ));
+  CALL_SUBTEST(( test_fft<float>(2*3*4*5*7) )); CALL_SUBTEST(( test_fft<double>(2*3*4*5*7) )); CALL_SUBTEST(( test_fft<long double>(2*3*4*5*7) ));
+}