# HG changeset patch
# User Chris Cannam <c.cannam@qmul.ac.uk>
# Date 1381833498 -3600
# Node ID 3c7338aff6a83c4b9b81f380332a6d4a871eb017
# Parent  8519e43841df27346a53c07b42b5fb3dbd1dc509
Drop in kissfft to replace the "old" fft, and add tests for newly-supported sizes

diff -r 8519e43841df -r 3c7338aff6a8 dsp/transforms/FFT.cpp
--- a/dsp/transforms/FFT.cpp	Tue Oct 15 09:39:58 2013 +0100
+++ b/dsp/transforms/FFT.cpp	Tue Oct 15 11:38:18 2013 +0100
@@ -11,195 +11,171 @@
 
 #include "maths/MathUtilities.h"
 
+#include "kiss_fft.h"
+#include "kiss_fftr.h"
+
 #include <cmath>
 
 #include <iostream>
 
+#include <stdexcept>
+
+class FFT::D
+{
+public:
+    D(int n) : m_n(n) {
+        m_planf = kiss_fft_alloc(m_n, 0, NULL, NULL);
+        m_plani = kiss_fft_alloc(m_n, 1, NULL, NULL);
+        m_kin = new kiss_fft_cpx[m_n];
+        m_kout = new kiss_fft_cpx[m_n];
+    }
+
+    ~D() {
+        kiss_fft_free(m_planf);
+        kiss_fft_free(m_plani);
+        delete[] m_kin;
+        delete[] m_kout;
+    }
+
+    void process(bool inverse,
+                 const double *ri,
+                 const double *ii,
+                 double *ro,
+                 double *io) {
+
+        for (int i = 0; i < m_n; ++i) {
+            m_kin[i].r = ri[i];
+            m_kin[i].i = (ii ? ii[i] : 0.0);
+        }
+
+        if (!inverse) {
+
+            kiss_fft(m_planf, m_kin, m_kout);
+
+            for (int i = 0; i < m_n; ++i) {
+                ro[i] = m_kout[i].r;
+                io[i] = m_kout[i].i;
+            }
+
+        } else {
+
+            kiss_fft(m_plani, m_kin, m_kout);
+
+            double scale = 1.0 / m_n;
+
+            for (int i = 0; i < m_n; ++i) {
+                ro[i] = m_kout[i].r * scale;
+                io[i] = m_kout[i].i * scale;
+            }
+        }
+    }
+    
+private:
+    int m_n;
+    kiss_fft_cfg m_planf;
+    kiss_fft_cfg m_plani;
+    kiss_fft_cpx *m_kin;
+    kiss_fft_cpx *m_kout;
+};        
+
 FFT::FFT(int n) :
-    m_n(n)
+    m_d(new D(n))
 {
-    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()
 {
-
+    delete m_d;
 }
 
+void
+FFT::process(bool inverse,
+             const double *p_lpRealIn, const double *p_lpImagIn,
+             double *p_lpRealOut, double *p_lpImagOut)
+{
+    m_d->process(inverse,
+                 p_lpRealIn, p_lpImagIn,
+                 p_lpRealOut, p_lpImagOut);
+}
+    
+class FFTReal::D
+{
+public:
+    D(int n) : m_n(n) {
+        if (n % 2) {
+            throw std::invalid_argument
+                ("nsamples must be even in FFTReal constructor");
+        }
+        m_planf = kiss_fftr_alloc(m_n, 0, NULL, NULL);
+        m_plani = kiss_fftr_alloc(m_n, 1, NULL, NULL);
+        m_c = new kiss_fft_cpx[m_n];
+    }
+
+    ~D() {
+        kiss_fftr_free(m_planf);
+        kiss_fftr_free(m_plani);
+        delete[] m_c;
+    }
+
+    void forward(const double *ri, double *ro, double *io) {
+
+        kiss_fftr(m_planf, ri, m_c);
+
+        for (int i = 0; i <= m_n/2; ++i) {
+            ro[i] = m_c[i].r;
+            io[i] = m_c[i].i;
+        }
+
+        for (int i = 0; i + 1 < m_n/2; ++i) {
+            ro[m_n - i - 1] =  ro[i + 1];
+            io[m_n - i - 1] = -io[i + 1];
+        }
+    }
+
+    void inverse(const double *ri, const double *ii, double *ro) {
+
+        for (int i = 0; i < m_n; ++i) {
+            m_c[i].r = ri[i];
+            m_c[i].i = ii[i];
+        }
+        
+        kiss_fftri(m_plani, m_c, ro);
+
+        double scale = 1.0 / m_n;
+
+        for (int i = 0; i < m_n; ++i) {
+            ro[i] *= scale;
+        }
+    }
+
+private:
+    int m_n;
+    kiss_fftr_cfg m_planf;
+    kiss_fftr_cfg m_plani;
+    kiss_fft_cpx *m_c;
+};
+
 FFTReal::FFTReal(int n) :
-    m_n(n),
-    m_fft(new FFT(n)),
-    m_r(new double[n]),
-    m_i(new double[n]),
-    m_discard(new double[n])
+    m_d(new D(n)) 
 {
-    for (int i = 0; i < n; ++i) {
-        m_r[i] = 0;
-        m_i[i] = 0;
-        m_discard[i] = 0;
-    }
 }
 
 FFTReal::~FFTReal()
 {
-    delete m_fft;
-    delete[] m_discard;
-    delete[] m_r;
-    delete[] m_i;
+    delete m_d;
 }
 
 void
-FFTReal::forward(const double *realIn,
-                 double *realOut, double *imagOut)
+FFTReal::forward(const double *ri, double *ro, double *io)
 {
-    m_fft->process(false, realIn, 0, realOut, imagOut);
+    m_d->forward(ri, ro, io);
 }
 
 void
-FFTReal::inverse(const double *realIn, const double *imagIn,
-                 double *realOut)
+FFTReal::inverse(const double *ri, const double *ii, double *ro)
 {
-    for (int i = 0; i < m_n/2 + 1; ++i) {
-        m_r[i] = realIn[i];
-        m_i[i] = imagIn[i];
-        if (i > 0 && i < m_n/2) {
-            m_r[m_n - i] = realIn[i];
-            m_i[m_n - i] = -imagIn[i];
-        }
-    }
-    m_fft->process(true, m_r, m_i, realOut, m_discard);
+    m_d->inverse(ri, ii, ro);
 }
 
-static int numberOfBitsNeeded(int p_nSamples)
-{	
-    int i;
 
-    if( p_nSamples < 2 )
-    {
-	return 0;
-    }
-
-    for ( i=0; ; i++ )
-    {
-	if( p_nSamples & (1 << i) ) return i;
-    }
-}
-
-static int reverseBits(int p_nIndex, int p_nBits)
-{
-    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;
-
-//    std::cerr << "FFT::process(" << m_n << "," << p_bInverseTransform << ")" << std::endl;
-
-    int NumBits;
-    int i, j, k, n;
-    int BlockSize, BlockEnd;
-
-    double angle_numerator = 2.0 * M_PI;
-    double tr, ti;
-
-    if( !MathUtilities::isPowerOfTwo(m_n) )
-    {
-        std::cerr << "ERROR: FFT::process: Non-power-of-two FFT size "
-                  << m_n << " not supported in this implementation"
-                  << std::endl;
-	return;
-    }
-
-    if( p_bInverseTransform ) angle_numerator = -angle_numerator;
-
-    NumBits = numberOfBitsNeeded ( m_n );
-
-
-
-    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 <= 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 < m_n; i += BlockSize )
-	{
-
-	    ar[2] = cm2;
-	    ar[1] = cm1;
-
-	    ai[2] = sm2;
-	    ai[1] = sm1;
-
-	    for ( j=i, n=0; n < BlockEnd; j++, n++ )
-	    {
-
-		ar[0] = w*ar[1] - ar[2];
-		ar[2] = ar[1];
-		ar[1] = ar[0];
-
-		ai[0] = w*ai[1] - ai[2];
-		ai[2] = ai[1];
-		ai[1] = ai[0];
-
-		k = j + BlockEnd;
-		tr = ar[0]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k];
-		ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k];
-
-		p_lpRealOut[k] = p_lpRealOut[j] - tr;
-		p_lpImagOut[k] = p_lpImagOut[j] - ti;
-
-		p_lpRealOut[j] += tr;
-		p_lpImagOut[j] += ti;
-
-	    }
-	}
-
-	BlockEnd = BlockSize;
-
-    }
-
-
-    if( p_bInverseTransform )
-    {
-	double denom = (double)m_n;
-
-	for ( i=0; i < m_n; i++ )
-	{
-	    p_lpRealOut[i] /= denom;
-	    p_lpImagOut[i] /= denom;
-	}
-    }
-}
-
+    
diff -r 8519e43841df -r 3c7338aff6a8 dsp/transforms/FFT.h
--- a/dsp/transforms/FFT.h	Tue Oct 15 09:39:58 2013 +0100
+++ b/dsp/transforms/FFT.h	Tue Oct 15 11:38:18 2013 +0100
@@ -14,8 +14,8 @@
 public:
     /**
      * Construct an FFT object to carry out complex-to-complex
-     * transforms of size nsamples. nsamples must be a power of two in
-     * this implementation.
+     * transforms of size nsamples. nsamples does not have to be a
+     * power of two.
      */
     FFT(int nsamples);
     ~FFT();
@@ -39,7 +39,8 @@
                  double *realOut, double *imagOut);
     
 private:
-    int m_n;
+    class D;
+    D *m_d;
 };
 
 class FFTReal
@@ -47,8 +48,9 @@
 public:
     /**
      * Construct an FFT object to carry out real-to-complex transforms
-     * of size nsamples. nsamples must be a power of two in this
-     * implementation.
+     * of size nsamples. nsamples does not have to be a power of two,
+     * but it does have to be even. (A std::invalid_argument exception
+     * will be thrown if nsamples is odd.)
      */
     FFTReal(int nsamples);
     ~FFTReal();
@@ -83,11 +85,8 @@
                  double *realOut);
 
 private:
-    int m_n;
-    FFT *m_fft;
-    double *m_r;
-    double *m_i;
-    double *m_discard;
+    class D;
+    D *m_d;
 };    
 
 #endif
diff -r 8519e43841df -r 3c7338aff6a8 qm-dsp.pro
--- a/qm-dsp.pro	Tue Oct 15 09:39:58 2013 +0100
+++ b/qm-dsp.pro	Tue Oct 15 11:38:18 2013 +0100
@@ -5,6 +5,12 @@
 OBJECTS_DIR = tmp_obj
 MOC_DIR = tmp_moc
 
+# This one is vital
+QMAKE_CFLAGS += -Dkiss_fft_scalar=double
+QMAKE_CXXFLAGS += -Dkiss_fft_scalar=double
+
+INCLUDEPATH += ./ext/kissfft/tools ./ext/kissfft
+
 linux-g++* {
     QMAKE_CXXFLAGS_RELEASE += -DNDEBUG -O3 -fno-exceptions -fPIC -ffast-math -msse -mfpmath=sse -ftree-vectorize -fomit-frame-pointer
     DEFINES += USE_PTHREADS
@@ -82,7 +88,9 @@
            maths/pca/pca.h \
            thread/AsynchronousTask.h \
            thread/BlockAllocator.h \
-           thread/Thread.h
+           thread/Thread.h \
+           ext/kissfft/kiss_fft.h \
+           ext/kissfft/tools/kiss_fftr.h
 SOURCES += base/Pitch.cpp \
            base/KaiserWindow.cpp \
            base/SincWindow.cpp \
@@ -117,4 +125,6 @@
            maths/KLDivergence.cpp \
            maths/MathUtilities.cpp \
            maths/pca/pca.c \
-           thread/Thread.cpp
+           thread/Thread.cpp \
+           ext/kissfft/kiss_fft.c \
+           ext/kissfft/tools/kiss_fftr.c
diff -r 8519e43841df -r 3c7338aff6a8 tests/Makefile
--- a/tests/Makefile	Tue Oct 15 09:39:58 2013 +0100
+++ b/tests/Makefile	Tue Oct 15 11:38:18 2013 +0100
@@ -7,8 +7,10 @@
 
 TESTS	:= test-mathutilities test-window test-fft test-pvoc
 
+#VG	:= valgrind
+
 all: $(TESTS)
-	for t in $(TESTS); do echo "Running $$t"; ./"$$t" || exit 1; done
+	for t in $(TESTS); do echo "Running $$t"; $(VG) ./"$$t" || exit 1; done
 
 test-mathutilities: TestMathUtilities.o $(LIBS)
 	$(CXX) -o $@ $^ $(LDFLAGS)
diff -r 8519e43841df -r 3c7338aff6a8 tests/TestFFT.cpp
--- a/tests/TestFFT.cpp	Tue Oct 15 09:39:58 2013 +0100
+++ b/tests/TestFFT.cpp	Tue Oct 15 11:38:18 2013 +0100
@@ -7,6 +7,8 @@
 
 #include <boost/test/unit_test.hpp>
 
+#include <stdexcept>
+
 BOOST_AUTO_TEST_SUITE(TestFFT)
 
 #define COMPARE_CONST(a, n) \
@@ -316,5 +318,48 @@
     COMPARE_ARRAY(back, in);
 }
 
+BOOST_AUTO_TEST_CASE(sizes) 
+{
+    // Complex supports any size. A single test with an odd size
+    // will do here, without getting too much into our expectations
+    // about supported butterflies etc
+
+    double in[] = { 1, 1, 1 };
+    double re[] = { 999, 999, 999 };
+    double im[] = { 999, 999, 999 };
+    FFT(3).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 3.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[3];
+    double backim[3];
+    FFT(3).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_sizes)
+{
+    // Real supports any even size, but not odd ones
+
+    BOOST_CHECK_THROW(FFTReal r(3), std::invalid_argument);
+
+    double in[] = { 1, 1, 1, 1, 1, 1 };
+    double re[] = { 999, 999, 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999, 999, 999 };
+    FFTReal(6).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 6.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    BOOST_CHECK_EQUAL(re[4], 0.0);
+    BOOST_CHECK_EQUAL(re[5], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[6];
+    FFTReal(6).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
 BOOST_AUTO_TEST_SUITE_END()
 
diff -r 8519e43841df -r 3c7338aff6a8 tests/TestPhaseVocoder.cpp
--- a/tests/TestPhaseVocoder.cpp	Tue Oct 15 09:39:58 2013 +0100
+++ b/tests/TestPhaseVocoder.cpp	Tue Oct 15 11:38:18 2013 +0100
@@ -164,14 +164,16 @@
     //   phase of -pi/2 plus advance of 2pi). It has the same measured
     //   phase as last time around, -pi/2, which is consistent with
     //   the expected phase, so the unwrapped phase is 3pi/2.
-    //!!!
-    // * Bin 3 is a bit of a puzzle -- it has an effectively zero
-    //   magnitude but a non-zero measured phase. Spectral leakage?
+    //
+    // * Bin 3 I don't really know about -- the magnitude here is 0,
+    //   but we get non-zero measured phase whose sign is
+    //   implementation-dependent
     //
     // * Bin 4 (Nyquist) has expected phase 4*pi, measured phase 0,
     //   hence error 0 and unwrapped phase 4*pi.
 
-    double phaseExpected1[] = { 999, 0, -M_PI, -M_PI/2, M_PI, 0, 999 };
+    phase[1+3] = 0.0; // Because we aren't testing for this one
+    double phaseExpected1[] = { 999, 0, -M_PI, -M_PI/2, 0, 0, 999 };
     COMPARE_ARRAY(phase, phaseExpected1);
 
     double unwExpected1[] = { 999, 0, M_PI, 3*M_PI/2, 3*M_PI, 4*M_PI, 999 };