# HG changeset patch
# User Chris Cannam
# Date 1342091038 -3600
# Node ID fb862b3418f3f6d5a92f50bd3dd88067a4fa30c1
# Parent  c3ce99b666b4f85b9284e81959e498c5285077d1
Update to use Vamp SDK 2.4 FFT implementation

diff -r c3ce99b666b4 -r fb862b3418f3 CepstrumPitchTracker.cpp
--- a/CepstrumPitchTracker.cpp	Tue Jul 10 22:16:00 2012 +0100
+++ b/CepstrumPitchTracker.cpp	Thu Jul 12 12:03:58 2012 +0100
@@ -22,6 +22,8 @@
 
 #include "CepstrumPitchTracker.h"
 
+#include "vamp-sdk/FFT.h"
+
 #include <vector>
 #include <algorithm>
 
@@ -530,7 +532,7 @@
     magmean /= hs;
     double threshold = 0.1; // for magmean
     
-    fft(bs, true, logmag, 0, rawcep, io);
+    Vamp::FFT::inverse(bs, logmag, 0, rawcep, io);
     
     delete[] logmag;
     delete[] io;
@@ -649,99 +651,3 @@
     }
     return fs;
 }
-
-void
-CepstrumPitchTracker::fft(unsigned int n, bool inverse,
-                          double *ri, double *ii, double *ro, double *io)
-{
-    if (!ri || !ro || !io) return;
-
-    unsigned int bits;
-    unsigned int i, j, k, m;
-    unsigned int blockSize, blockEnd;
-
-    double tr, ti;
-
-    if (n < 2) return;
-    if (n & (n-1)) return;
-
-    double angle = 2.0 * M_PI;
-    if (inverse) angle = -angle;
-
-    for (i = 0; ; ++i) {
-	if (n & (1 << i)) {
-	    bits = i;
-	    break;
-	}
-    }
-
-    int table[n];
-
-    for (i = 0; i < n; ++i) {
-        m = i;
-        for (j = k = 0; j < bits; ++j) {
-            k = (k << 1) | (m & 1);
-            m >>= 1;
-        }
-        table[i] = k;
-    }
-
-    if (ii) {
-	for (i = 0; i < n; ++i) {
-	    ro[table[i]] = ri[i];
-	    io[table[i]] = ii[i];
-	}
-    } else {
-	for (i = 0; i < n; ++i) {
-	    ro[table[i]] = ri[i];
-	    io[table[i]] = 0.0;
-	}
-    }
-
-    blockEnd = 1;
-
-    for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
-
-	double delta = angle / (double)blockSize;
-	double sm2 = -sin(-2 * delta);
-	double sm1 = -sin(-delta);
-	double cm2 = cos(-2 * delta);
-	double cm1 = cos(-delta);
-	double w = 2 * cm1;
-	double ar[3], ai[3];
-
-	for (i = 0; i < n; i += blockSize) {
-
-	    ar[2] = cm2;
-	    ar[1] = cm1;
-
-	    ai[2] = sm2;
-	    ai[1] = sm1;
-
-	    for (j = i, m = 0; m < blockEnd; j++, m++) {
-
-		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] * ro[k] - ai[0] * io[k];
-		ti = ar[0] * io[k] + ai[0] * ro[k];
-
-		ro[k] = ro[j] - tr;
-		io[k] = io[j] - ti;
-
-		ro[j] += tr;
-		io[j] += ti;
-	    }
-	}
-
-	blockEnd = blockSize;
-    }
-}
-
-
diff -r c3ce99b666b4 -r fb862b3418f3 CepstrumPitchTracker.h
--- a/CepstrumPitchTracker.h	Tue Jul 10 22:16:00 2012 +0100
+++ b/CepstrumPitchTracker.h	Thu Jul 12 12:03:58 2012 +0100
@@ -130,8 +130,6 @@
     void filter(const double *in, double *out);
     double cubicInterpolate(const double[4], double);
     double findInterpolatedPeak(const double *in, int maxbin);
-    void fft(unsigned int n, bool inverse,
-             double *ri, double *ii, double *ro, double *io);
 };
 
 #endif
diff -r c3ce99b666b4 -r fb862b3418f3 SimpleCepstrum.cpp
--- a/SimpleCepstrum.cpp	Tue Jul 10 22:16:00 2012 +0100
+++ b/SimpleCepstrum.cpp	Thu Jul 12 12:03:58 2012 +0100
@@ -22,6 +22,8 @@
 
 #include "SimpleCepstrum.h"
 
+#include "vamp-sdk/FFT.h"
+
 #include <vector>
 #include <algorithm>
 
@@ -29,6 +31,10 @@
 #include <cmath>
 #include <complex>
 
+#if ( VAMP_SDK_MAJOR_VERSION < 2 || ( VAMP_SDK_MAJOR_VERSION == 2 && VAMP_SDK_MINOR_VERSION < 4 ) )
+#error Vamp SDK version 2.4 or newer required
+#endif
+
 using std::string;
 
 SimpleCepstrum::SimpleCepstrum(float inputSampleRate) :
@@ -654,7 +660,7 @@
     double *io = new double[bs];
 
     //!!! This is only right if the previous transform was an inverse one!
-    fft(bs, false, ecep, 0, env, io);
+    Vamp::FFT::forward(bs, ecep, 0, env, io);
 
     for (int i = 0; i < hs; ++i) {
         env[i] = exp(env[i]);
@@ -724,11 +730,11 @@
         if (m_method == InverseSymmetric ||
             m_method == InverseAsymmetric) {
 
-            fft(bs, true, logmag, 0, rawcep, io);
+            Vamp::FFT::inverse(bs, logmag, 0, rawcep, io);
 
         } else {
 
-            fft(bs, false, logmag, 0, rawcep, io);
+            Vamp::FFT::forward(bs, logmag, 0, rawcep, io);
 
             if (m_method == ForwardDifference) {
                 for (int i = 0; i < hs; ++i) {
@@ -761,7 +767,7 @@
             }
         }
 
-        fft(bs, true, ri, ii, rawcep, io);
+        Vamp::FFT::inverse(bs, ri, ii, rawcep, io);
 
         delete[] ri;
         delete[] ii;
@@ -802,112 +808,3 @@
     FeatureSet fs;
     return fs;
 }
-
-void
-SimpleCepstrum::fft(unsigned int n, bool inverse,
-                    double *ri, double *ii, double *ro, double *io)
-{
-    if (!ri || !ro || !io) return;
-
-    unsigned int bits;
-    unsigned int i, j, k, m;
-    unsigned int blockSize, blockEnd;
-
-    double tr, ti;
-
-    if (n < 2) return;
-    if (n & (n-1)) return;
-
-    double angle = 2.0 * M_PI;
-    if (inverse) angle = -angle;
-
-    for (i = 0; ; ++i) {
-	if (n & (1 << i)) {
-	    bits = i;
-	    break;
-	}
-    }
-
-    static unsigned int tableSize = 0;
-    static int *table = 0;
-
-    if (tableSize != n) {
-
-	delete[] table;
-
-	table = new int[n];
-
-	for (i = 0; i < n; ++i) {
-	
-	    m = i;
-
-	    for (j = k = 0; j < bits; ++j) {
-		k = (k << 1) | (m & 1);
-		m >>= 1;
-	    }
-
-	    table[i] = k;
-	}
-
-	tableSize = n;
-    }
-
-    if (ii) {
-	for (i = 0; i < n; ++i) {
-	    ro[table[i]] = ri[i];
-	    io[table[i]] = ii[i];
-	}
-    } else {
-	for (i = 0; i < n; ++i) {
-	    ro[table[i]] = ri[i];
-	    io[table[i]] = 0.0;
-	}
-    }
-
-    blockEnd = 1;
-
-    for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
-
-	double delta = angle / (double)blockSize;
-	double sm2 = -sin(-2 * delta);
-	double sm1 = -sin(-delta);
-	double cm2 = cos(-2 * delta);
-	double cm1 = cos(-delta);
-	double w = 2 * cm1;
-	double ar[3], ai[3];
-
-	for (i = 0; i < n; i += blockSize) {
-
-	    ar[2] = cm2;
-	    ar[1] = cm1;
-
-	    ai[2] = sm2;
-	    ai[1] = sm1;
-
-	    for (j = i, m = 0; m < blockEnd; j++, m++) {
-
-		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] * ro[k] - ai[0] * io[k];
-		ti = ar[0] * io[k] + ai[0] * ro[k];
-
-		ro[k] = ro[j] - tr;
-		io[k] = io[j] - ti;
-
-		ro[j] += tr;
-		io[j] += ti;
-	    }
-	}
-
-	blockEnd = blockSize;
-    }
-}
-
-
diff -r c3ce99b666b4 -r fb862b3418f3 SimpleCepstrum.h
--- a/SimpleCepstrum.h	Tue Jul 10 22:16:00 2012 +0100
+++ b/SimpleCepstrum.h	Thu Jul 12 12:03:58 2012 +0100
@@ -104,9 +104,6 @@
     void filter(const double *in, double *out);
     double cubicInterpolate(const double[4], double);
     double findInterpolatedPeak(const double *in, int maxbin);
-    void fft(unsigned int n, bool inverse,
-             double *ri, double *ii, double *ro, double *io);
-
     void addStatisticalOutputs(FeatureSet &fs, const double *data);
     void addEnvelopeOutputs(FeatureSet &fs, const float *const *inputBuffers,
                             const double *raw);