From d659fd9b148870e3c1e367139bec388142f2818e Mon Sep 17 00:00:00 2001
From: Mark Borgerding <mark@borgerding.net>
Date: Fri, 30 Oct 2009 20:26:30 -0400
Subject: [PATCH] moved real-half-spectrum reflection into Eigen::FFT

---
 unsupported/Eigen/FFT                       | 32 ++++++++++++++-------
 unsupported/Eigen/src/FFT/ei_fftw_impl.h    |  3 --
 unsupported/Eigen/src/FFT/ei_kissfft_impl.h |  2 --
 3 files changed, 22 insertions(+), 15 deletions(-)

diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT
index e2705abf6..b1d3d9f0e 100644
--- a/unsupported/Eigen/FFT
+++ b/unsupported/Eigen/FFT
@@ -88,6 +88,8 @@ class FFT
     void fwd( Complex * dst, const Scalar * src, int nfft)
     {
         m_impl.fwd(dst,src,nfft);
+        if ( HasFlag(HalfSpectrum) == false)
+          ReflectSpectrum(dst,nfft);
     }
 
     void fwd( Complex * dst, const Complex * src, int nfft)
@@ -108,15 +110,19 @@ class FFT
     template<typename InputDerived, typename ComplexDerived>
     void fwd( MatrixBase<ComplexDerived> & dst, const MatrixBase<InputDerived> & src)
     {
-        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((ei_is_same_type<typename ComplexDerived::Scalar, Complex>::ret),
-                            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)
-        dst.derived().resize( src.size() );
-        fwd( &dst[0],&src[0],src.size() );
+      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((ei_is_same_type<typename ComplexDerived::Scalar, Complex>::ret),
+            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 ( NumTraits< typename InputDerived::Scalar >::IsComplex == 0 && HasFlag(HalfSpectrum) )
+        dst.derived().resize( (src.size()>>1)+1);
+      else
+        dst.derived().resize(src.size());
+      fwd( &dst[0],&src[0],src.size() );
     }
 
     void inv( Complex * dst, const Complex * src, int nfft)
@@ -160,7 +166,6 @@ class FFT
       inv( &dst[0],&src[0],dst.size() );
     }
 
-
     // TODO: multi-dimensional FFTs
     
     // TODO: handle Eigen MatrixBase
@@ -176,6 +181,13 @@ class FFT
         *x++ *= s;
     }
 
+    void ReflectSpectrum(Complex * freq,int nfft)
+    {
+      int nhbins=(nfft>>1)+1;
+      for (int k=nhbins;k < nfft; ++k )
+        freq[k] = conj(freq[nfft-k]);
+    }
+
     impl_type m_impl;
     int m_flag;
 };
diff --git a/unsupported/Eigen/src/FFT/ei_fftw_impl.h b/unsupported/Eigen/src/FFT/ei_fftw_impl.h
index cdb5861e8..18473a29b 100644
--- a/unsupported/Eigen/src/FFT/ei_fftw_impl.h
+++ b/unsupported/Eigen/src/FFT/ei_fftw_impl.h
@@ -177,9 +177,6 @@
       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
diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
index 86ec7a1ca..859e7e6c9 100644
--- a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
+++ b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
@@ -322,8 +322,6 @@
 
           // place conjugate-symmetric half at the end for completeness
           // TODO: make this configurable ( opt-out )
-          for ( k=1;k < ncfft ; ++k )
-            dst[nfft-k] = conj(dst[k]);
           dst[0] = dc;
           dst[ncfft] = nyquist;
         }