scalar forward FFT optimized for even size, converts to cpx for odd

unsupported/Eigen/src/FFT/simple_fft_traits.h

@@ -24,6 +24,7 @@
 
 #include <complex>
 #include <vector>
+#include <iostream>
 
 namespace Eigen {
 
@@ -39,51 +40,54 @@ namespace Eigen {
     {
         prepare(nfft,false);
         work(0, dst, src, 1,1);
-        scale(dst);
     }
 
+    // real-to-complex forward FFT
+    // perform two FFTs of src even and src odd
+    // then twiddle to recombine them into the half-spectrum format
+    // then fill in the conjugate symmetric half
     void fwd( Complex * dst,const Scalar * src,int nfft) 
     {
         if ( nfft&1 ) {
             // use generic mode for odd
             prepare(nfft,false);
             work(0, dst, src, 1,1);
-            scale(dst);
         }else{
             int ncfft = nfft>>1;
             // use optimized mode for even real
-            prepare(nfft,false);
+            fwd( dst, reinterpret_cast<const Complex*> (src),ncfft);
-            work(0,dst, reinterpret_cast<const Complex*> (src),2,1);
+            make_real_twiddles(nfft);
+            std::cerr << "dst[0] = " << dst[0] << "\n";
             Complex dc = dst[0].real() +  dst[0].imag();
             Complex nyquist = dst[0].real() -  dst[0].imag();
-            
             int k;
-            for ( k=1;k <= ncfft/2 ; ++k ) {
+#if 0
-/**
+            using namespace std;
-        fpk    = st->tmpbuf[k];
+            cerr << "desired:\n";
-        fpnk.r =   st->tmpbuf[ncfft-k].r;
+            for ( k=1;k <= (ncfft>>1) ; ++k ) {
-        fpnk.i = - st->tmpbuf[ncfft-k].i;
+                Complex x = exp( Complex(0,-3.14159265358979323846264338327 * ((double) (k) / ncfft + .5) ) );
-
+                cerr << k << " " << x << "angle(x):" << arg(x) << "\n";
-        C_ADD( f1k, fpk , fpnk );
+            }
-        C_SUB( f2k, fpk , fpnk );
+            dc=0;
-        C_MUL( tw , f2k , st->super_twiddles[k-1]);
+            cerr << "twiddles:\n";
-
+            for (k=0;k<ncfft;++k) {
-        freqdata[k].r = HALF_OF(f1k.r + tw.r);
+                Complex x = m_twiddles[k];
-        freqdata[k].i = HALF_OF(f1k.i + tw.i);
+                cerr << k << " " << x << "angle(-x):" << arg(-x) << "\n";
-        freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r);
+            }
-        freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i);
+#endif
- */
+            for ( k=1;k <= (ncfft>>1) ; ++k ) {
                 Complex fpk = dst[k];
                 Complex fpnk = conj(dst[ncfft-k]);
-
                 Complex f1k = fpk + fpnk;
-                Complex f2k = fpnk - fpk;
+                Complex f2k = fpk - fpnk;
-                //Complex tw = f2k * exp( Complex(0,-3.14159265358979323846264338327 * ((double)k / ncfft + 1) ) ); // TODO repalce this with index into twiddles
+                //Complex tw = f2k * exp( Complex(0,-3.14159265358979323846264338327 * ((double) (k) / ncfft + .5) ) );
-                Complex tw = f2k * m_twiddles[2*k];;
+                Complex tw= f2k * m_realTwiddles[k-1];
-                
+
                 dst[k] =  (f1k + tw) * Scalar(.5);
                 dst[ncfft-k] =  conj(f1k -tw)*Scalar(.5);
             }
+ 
+
             // place conjugate-symmetric half at the end for completeness
             // TODO: make this configurable ( opt-out )
             for ( k=1;k < ncfft ; ++k )
@@ -98,55 +102,74 @@ namespace Eigen {
     {
         prepare(nfft,true);
         work(0, dst, src, 1,1);
-        scale(dst);
+        scale(dst, Scalar(1)/m_nfft );
     }
 
     void prepare(int nfft,bool inverse)
     {
-      if (m_nfft == nfft) {
+        make_twiddles(nfft,inverse);
-        // reuse the twiddles, conjugate if necessary
+        factorize(nfft);
-        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 scale(Complex *dst) 
+    void make_real_twiddles(int nfft)
+    {
+        int ncfft2 = nfft>>2;
+        if ( m_realTwiddles.size() != ncfft2) {
+            m_realTwiddles.resize(ncfft2);
+            int ncfft= nfft>>1;
+            for (int k=1;k<=ncfft2;++k) 
+                m_realTwiddles[k-1] = exp( Complex(0,-3.14159265358979323846264338327 * ((double) (k) / ncfft + .5) ) );
+        }
+    }
+
+    void make_twiddles(int nfft,bool inverse)
+    {
+        if ( m_twiddles.size() == 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] );
+        }else{
+            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_inverse = inverse;
+    }
+
+    void factorize(int nfft)
+    {
+        if (m_stageRadix.size()==0 || m_stageRadix[0] * m_stageRemainder[0] != nfft)
+        {
+            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);
+        }
+        m_nfft = nfft;
+    }
+
+    void scale(Complex *dst,Scalar s) 
     {
-      if (m_inverse) {
-        Scalar s = 1./m_nfft;
         for (int k=0;k<m_nfft;++k)
-          dst[k] *= s;
+            dst[k] *= s;
-      }
     }
 
     private:
@@ -349,6 +372,7 @@ namespace Eigen {
     int m_nfft;
     bool m_inverse;
     std::vector<Complex> m_twiddles;
+    std::vector<Complex> m_realTwiddles;
     std::vector<int> m_stageRadix;
     std::vector<int> m_stageRemainder;
   };

unsupported/test/FFT.cpp

@@ -41,6 +41,7 @@ complex<long double>  promote(long double x) { return complex<long double>( x);
     {
         long double totalpower=0;
         long double difpower=0;
+        cerr <<"idx\ttruth\t\tvalue\n";
         for (size_t k0=0;k0<fftbuf.size();++k0) {
             complex<long double> acc = 0;
             long double phinc = -2.*k0* M_PIl / timebuf.size();
@@ -51,7 +52,7 @@ complex<long double>  promote(long double x) { return complex<long double>( x);
             complex<long double> x = promote(fftbuf[k0]); 
             complex<long double> dif = acc - x;
             difpower += norm(dif);
-            cerr << k0 << ":" << acc << " " <<  x << endl;
+            cerr << k0 << "\t" << acc << "\t" <<  x << endl;
         }
         cerr << "rmse:" << sqrt(difpower/totalpower) << endl;
         return sqrt(difpower/totalpower);
@@ -108,6 +109,7 @@ void test_complex(int nfft)
 
 void test_FFT()
 {
+#if 0
   CALL_SUBTEST( test_complex<float>(32) ); CALL_SUBTEST( test_complex<double>(32) ); CALL_SUBTEST( test_complex<long double>(32) );
   CALL_SUBTEST( test_complex<float>(1024) ); CALL_SUBTEST( test_complex<double>(1024) ); CALL_SUBTEST( test_complex<long double>(1024) );
   CALL_SUBTEST( test_complex<float>(3*8) ); CALL_SUBTEST( test_complex<double>(3*8) ); CALL_SUBTEST( test_complex<long double>(3*8) );
@@ -115,9 +117,11 @@ void test_FFT()
   CALL_SUBTEST( test_complex<float>(2*3*4) ); CALL_SUBTEST( test_complex<double>(2*3*4) ); CALL_SUBTEST( test_complex<long double>(2*3*4) );
   CALL_SUBTEST( test_complex<float>(2*3*4*5) ); CALL_SUBTEST( test_complex<double>(2*3*4*5) ); CALL_SUBTEST( test_complex<long double>(2*3*4*5) );
   CALL_SUBTEST( test_complex<float>(2*3*4*5*7) ); CALL_SUBTEST( test_complex<double>(2*3*4*5*7) ); CALL_SUBTEST( test_complex<long double>(2*3*4*5*7) );
-/*
+#endif
+
+#if 1
   CALL_SUBTEST( test_scalar<float>(32) ); CALL_SUBTEST( test_scalar<double>(32) ); CALL_SUBTEST( test_scalar<long double>(32) );
   CALL_SUBTEST( test_scalar<float>(1024) ); CALL_SUBTEST( test_scalar<double>(1024) ); CALL_SUBTEST( test_scalar<long double>(1024) );
   CALL_SUBTEST( test_scalar<float>(2*3*4*5*7) ); CALL_SUBTEST( test_scalar<double>(2*3*4*5*7) ); CALL_SUBTEST( test_scalar<long double>(2*3*4*5*7) );
-  */
+#endif
 }