Mercurial > hg > qm-dsp

changeset 129:6ec45e85ed81 kissfft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Drop in kissfft to replace the "old" fft, and add tests for newly-supported sizes
author Chris Cannam
date Tue, 15 Oct 2013 11:38:18 +0100
parents 5023f521732c
children 2053a308bb4d
files dsp/transforms/FFT.cpp dsp/transforms/FFT.h qm-dsp.pro tests/Makefile tests/TestFFT.cpp tests/TestPhaseVocoder.cpp
diffstat 6 files changed, 214 insertions(+), 180 deletions(-) [+]
line wrap: on
line diff
--- 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;
-	}
-    }
-}
-
+    
--- 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
--- 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
--- 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)
--- 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()
 
--- 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 };