comparison dsp/transforms/FFT.cpp @ 64:6cb2b3cd5356

* Refactor FFT a little bit so as to separate construction and processing
  rather than have a single static method -- will make it easier to use a
  different implementation
* pull in KissFFT implementation (not hooked up yet)

 
 #include <cmath>
 
 #include <iostream>
 
-//////////////////////////////////////////////////////////////////////
-// Construction/Destruction
-//////////////////////////////////////////////////////////////////////
+#define USE_BUILTIN_FFT 1
 
-FFT::FFT()
+#ifdef USE_BUILTIN_FFT
+
+FFT::FFT(unsigned int n) :
+    m_n(n),
+    m_private(0)
 {
-
+    if( !MathUtilities::isPowerOfTwo(m_n) )
+    {
+        std::cerr << "ERROR: FFT: Non-power-of-two FFT size "
+                  << m_n << " not supported in this implementation"
+                  << std::endl;
+	return;
+    }
 }
 
 FFT::~FFT()
 {
 
 }
 
-void FFT::process(unsigned int p_nSamples, bool p_bInverseTransform,
-                  const double *p_lpRealIn, const double *p_lpImagIn,
-                  double *p_lpRealOut, double *p_lpImagOut)
+FFTReal::FFTReal(unsigned int n) :
+    m_n(n),
+    m_private(0)
 {
+    m_private = new FFT(m_n);
+}
 
+FFTReal::~FFTReal()
+{
+    delete (FFT *)m_private;
+}
+
+void
+FFTReal::process(bool inverse,
+                 const double *realIn,
+                 double *realOut, double *imagOut)
+{
+    ((FFT *)m_private)->process(inverse, realIn, 0, realOut, imagOut);
+}
+
+static unsigned int numberOfBitsNeeded(unsigned int p_nSamples)
+{	
+    int i;
+
+    if( p_nSamples < 2 )
+    {
+	return 0;
+    }
+
+    for ( i=0; ; i++ )
+    {
+	if( p_nSamples & (1 << i) ) return i;
+    }
+}
+
+static unsigned int reverseBits(unsigned int p_nIndex, unsigned int p_nBits)
+{
+    unsigned int i, rev;
+
+    for(i=rev=0; i < p_nBits; i++)
+    {
+	rev = (rev << 1) | (p_nIndex & 1);
+	p_nIndex >>= 1;
+    }
+
+    return rev;
+}
+
+void
+FFT::process(bool p_bInverseTransform,
+             const double *p_lpRealIn, const double *p_lpImagIn,
+             double *p_lpRealOut, double *p_lpImagOut)
+{
     if(!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return;
-
 
     unsigned int NumBits;
     unsigned int i, j, k, n;
     unsigned int BlockSize, BlockEnd;
 
     double angle_numerator = 2.0 * M_PI;
     double tr, ti;
 
-    if( !MathUtilities::isPowerOfTwo(p_nSamples) )
+    if( !MathUtilities::isPowerOfTwo(m_n) )
     {
         std::cerr << "ERROR: FFT::process: Non-power-of-two FFT size "
-                  << p_nSamples << " not supported in this implementation"
+                  << m_n << " not supported in this implementation"
                   << std::endl;
 	return;
     }
 
     if( p_bInverseTransform ) angle_numerator = -angle_numerator;
 
-    NumBits = numberOfBitsNeeded ( p_nSamples );
+    NumBits = numberOfBitsNeeded ( m_n );
 
 
-    for( i=0; i < p_nSamples; i++ )
+    for( i=0; i < m_n; i++ )
     {
 	j = reverseBits ( i, NumBits );
 	p_lpRealOut[j] = p_lpRealIn[i];
 	p_lpImagOut[j] = (p_lpImagIn == 0) ? 0.0 : p_lpImagIn[i];
     }
 
 
     BlockEnd = 1;
-    for( BlockSize = 2; BlockSize <= p_nSamples; BlockSize <<= 1 )
+    for( BlockSize = 2; BlockSize <= m_n; BlockSize <<= 1 )
     {
 	double delta_angle = angle_numerator / (double)BlockSize;
 	double sm2 = -sin ( -2 * delta_angle );
 	double sm1 = -sin ( -delta_angle );
 	double cm2 = cos ( -2 * delta_angle );
 	double cm1 = cos ( -delta_angle );
 	double w = 2 * cm1;
 	double ar[3], ai[3];
 
-	for( i=0; i < p_nSamples; i += BlockSize )
+	for( i=0; i < m_n; i += BlockSize )
 	{
 
 	    ar[2] = cm2;
 	    ar[1] = cm1;
 

     }
 
 
     if( p_bInverseTransform )
     {
-	double denom = (double)p_nSamples;
+	double denom = (double)m_n;
 
-	for ( i=0; i < p_nSamples; i++ )
+	for ( i=0; i < m_n; i++ )
 	{
 	    p_lpRealOut[i] /= denom;
 	    p_lpImagOut[i] /= denom;
 	}
     }
 }
 
-unsigned int FFT::numberOfBitsNeeded(unsigned int p_nSamples)
-{	
-    int i;
+#else
 
-    if( p_nSamples < 2 )
-    {
-	return 0;
-    }
+#include "kissfft/kiss_fft.h"
+#include "kissfft/kiss_fftr.h"
 
-    for ( i=0; ; i++ )
-    {
-	if( p_nSamples & (1 << i) ) return i;
-    }
-}
-
-unsigned int FFT::reverseBits(unsigned int p_nIndex, unsigned int p_nBits)
-{
-    unsigned int i, rev;
-
-    for(i=rev=0; i < p_nBits; i++)
-    {
-	rev = (rev << 1) | (p_nIndex & 1);
-	p_nIndex >>= 1;
-    }
-
-    return rev;
-}
+#endif