# HG changeset patch
# User Chris Cannam
# Date 1479742383 0
# Node ID 604f369c247a7af9448bd17a060315379d01ba66
# Parent  f50c0bbe90960b87f6ab3fd1661f5a9fb04c877a# Parent  bac86d3fc6c95e4af0ef95b664538f2517b46861
Merge from branch bqfft

diff -r f50c0bbe9096 -r 604f369c247a base/Window.h
--- a/base/Window.h	Fri Nov 18 23:30:15 2016 +0000
+++ b/base/Window.h	Mon Nov 21 15:33:03 2016 +0000
@@ -22,6 +22,9 @@
 #include <map>
 #include <cstdlib>
 
+#include <bqvec/VectorOps.h>
+#include <bqvec/Allocators.h>
+
 #include "system/System.h"
 
 enum WindowType {
@@ -47,8 +50,12 @@
      * than symmetrical. (A window of size N is equivalent to a
      * symmetrical window of size N+1 with the final element missing.)
      */
-    Window(WindowType type, int size) : m_type(type), m_size(size) { encache(); }
-    Window(const Window &w) : m_type(w.m_type), m_size(w.m_size) { encache(); }
+    Window(WindowType type, int size) : m_type(type), m_size(size), m_cache(0) {
+        encache();
+    }
+    Window(const Window &w) : m_type(w.m_type), m_size(w.m_size), m_cache(0) {
+        encache();
+    }
     Window &operator=(const Window &w) {
 	if (&w == this) return *this;
 	m_type = w.m_type;
@@ -56,11 +63,16 @@
 	encache();
 	return *this;
     }
-    virtual ~Window() { delete[] m_cache; }
+    virtual ~Window() {
+        breakfastquay::deallocate(m_cache);
+    }
     
-    void cut(T *src) const { cut(src, src); }
-    void cut(T *src, T *dst) const {
-	for (int i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i];
+    inline void cut(T *const BQ_R__ block) const {
+        breakfastquay::v_multiply(block, m_cache, m_size);
+    }
+
+    inline void cut(const T *const BQ_R__ src, T *const BQ_R__ dst) const {
+        breakfastquay::v_multiply(dst, src, m_cache, m_size);
     }
 
     T getArea() { return m_area; }
@@ -78,7 +90,7 @@
 protected:
     WindowType m_type;
     int m_size;
-    T *m_cache;
+    T *BQ_R__ m_cache;
     T m_area;
     
     void encache();
@@ -88,41 +100,42 @@
 template <typename T>
 void Window<T>::encache()
 {
+    if (!m_cache) m_cache = breakfastquay::allocate<T>(m_size);
+
     const int n = m_size;
-    T *mult = new T[n];
+    breakfastquay::v_set(m_cache, T(1.0), n);
     int i;
-    for (i = 0; i < n; ++i) mult[i] = 1.0;
 
     switch (m_type) {
 		
     case RectangularWindow:
 	for (i = 0; i < n; ++i) {
-	    mult[i] *= T(0.5);
+	    m_cache[i] *= T(0.5);
 	}
 	break;
 	    
     case BartlettWindow:
 	for (i = 0; i < n/2; ++i) {
-	    mult[i] *= T(i) / T(n/2);
-	    mult[i + n/2] *= T(1.0) - T(i) / T(n/2);
+	    m_cache[i] *= T(i) / T(n/2);
+	    m_cache[i + n/2] *= T(1.0) - T(i) / T(n/2);
 	}
 	break;
 	    
     case HammingWindow:
-        cosinewin(mult, 0.54, 0.46, 0.0, 0.0);
+        cosinewin(m_cache, 0.54, 0.46, 0.0, 0.0);
 	break;
 	    
     case HanningWindow:
-        cosinewin(mult, 0.50, 0.50, 0.0, 0.0);
+        cosinewin(m_cache, 0.50, 0.50, 0.0, 0.0);
 	break;
 	    
     case BlackmanWindow:
-        cosinewin(mult, 0.42, 0.50, 0.08, 0.0);
+        cosinewin(m_cache, 0.42, 0.50, 0.08, 0.0);
 	break;
 	    
     case GaussianWindow:
 	for (i = 0; i < n; ++i) {
-            mult[i] *= T(pow(2, - pow((i - (n-1)/2.0) / ((n-1)/2.0 / 3), 2)));
+            m_cache[i] *= T(pow(2, - pow((i - (n-1)/2.0) / ((n-1)/2.0 / 3), 2)));
 	}
 	break;
 	    
@@ -131,29 +144,27 @@
         int N = n-1;
         for (i = 0; i < N/4; ++i) {
             T m = T(2 * pow(1.0 - (T(N)/2 - T(i)) / (T(N)/2), 3));
-            mult[i] *= m;
-            mult[N-i] *= m;
+            m_cache[i] *= m;
+            m_cache[N-i] *= m;
         }
         for (i = N/4; i <= N/2; ++i) {
             int wn = i - N/2;
             T m = T(1.0 - 6 * pow(T(wn) / (T(N)/2), 2) * (1.0 - T(abs(wn)) / (T(N)/2)));
-            mult[i] *= m;
-            mult[N-i] *= m;
+            m_cache[i] *= m;
+            m_cache[N-i] *= m;
         }            
         break;
     }
 
     case NuttallWindow:
-        cosinewin(mult, 0.3635819, 0.4891775, 0.1365995, 0.0106411);
+        cosinewin(m_cache, 0.3635819, 0.4891775, 0.1365995, 0.0106411);
 	break;
 
     case BlackmanHarrisWindow:
-        cosinewin(mult, 0.35875, 0.48829, 0.14128, 0.01168);
+        cosinewin(m_cache, 0.35875, 0.48829, 0.14128, 0.01168);
         break;
     }
 	
-    m_cache = mult;
-
     m_area = 0;
     for (int i = 0; i < n; ++i) {
         m_area += m_cache[i];
diff -r f50c0bbe9096 -r 604f369c247a data/fft/FFTapi.cpp
--- a/data/fft/FFTapi.cpp	Fri Nov 18 23:30:15 2016 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,226 +0,0 @@
-/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
-
-/*
-    Sonic Visualiser
-    An audio file viewer and annotation editor.
-    Centre for Digital Music, Queen Mary, University of London.
-    This file copyright 2006 Chris Cannam and QMUL.
-    FFT code from Don Cross's public domain FFT implementation.
-    
-    This program is free software; you can redistribute it and/or
-    modify it under the terms of the GNU General Public License as
-    published by the Free Software Foundation; either version 2 of the
-    License, or (at your option) any later version.  See the file
-    COPYING included with this distribution for more information.
-*/
-
-#include "FFTapi.h"
-
-std::mutex FFTForward::m_mutex;
-
-#ifndef HAVE_FFTW3F
-
-#include <cmath>
-#include <iostream>
-
-void
-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 = 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;
-    }
-
-    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;
-    }
-
-/* fftw doesn't normalise, so nor will we
-
-    if (inverse) {
-
-	double denom = (double)n;
-
-	for (i = 0; i < n; i++) {
-	    ro[i] /= denom;
-	    io[i] /= denom;
-	}
-    }
-*/
-    delete[] table;
-}
-
-struct fftf_plan_ {
-    int size;
-    int inverse;
-    float *real;
-    fftf_complex *cplx;
-};
-
-fftf_plan
-fftf_plan_dft_r2c_1d(int n, float *in, fftf_complex *out, unsigned)
-{
-    if (n < 2) return 0;
-    if (n & (n-1)) return 0;
-    
-    fftf_plan_ *plan = new fftf_plan_;
-    plan->size = n;
-    plan->inverse = 0;
-    plan->real = in;
-    plan->cplx = out;
-    return plan;
-}
-
-fftf_plan
-fftf_plan_dft_c2r_1d(int n, fftf_complex *in, float *out, unsigned)
-{
-    if (n < 2) return 0;
-    if (n & (n-1)) return 0;
-    
-    fftf_plan_ *plan = new fftf_plan_;
-    plan->size = n;
-    plan->inverse = 1;
-    plan->real = out;
-    plan->cplx = in;
-    return plan;
-}
-
-void
-fftf_destroy_plan(fftf_plan p)
-{
-    delete p;
-}
-
-void
-fftf_execute(const fftf_plan p)
-{
-    float *real = p->real;
-    fftf_complex *cplx = p->cplx;
-    int n = p->size;
-    int forward = !p->inverse;
-
-    double *ri = new double[n];
-    double *ro = new double[n];
-    double *io = new double[n];
-
-    double *ii = 0;
-    if (!forward) ii = new double[n];
-
-    if (forward) {
-        for (int i = 0; i < n; ++i) {
-            ri[i] = real[i];
-        }
-    } else {
-        for (int i = 0; i < n/2+1; ++i) {
-            ri[i] = cplx[i][0];
-            ii[i] = cplx[i][1];
-            if (i > 0) {
-                ri[n-i] = ri[i];
-                ii[n-i] = -ii[i];
-            }
-        }
-    }
-
-    fft(n, !forward, ri, ii, ro, io);
-
-    if (forward) {
-        for (int i = 0; i < n/2+1; ++i) {
-            cplx[i][0] = ro[i];
-            cplx[i][1] = io[i];
-        }
-    } else {
-        for (int i = 0; i < n; ++i) {
-            real[i] = ro[i];
-        }
-    }
-
-    delete[] ri;
-    delete[] ro;
-    delete[] io;
-    if (ii) delete[] ii;
-}
-
-#endif
diff -r f50c0bbe9096 -r 604f369c247a data/fft/FFTapi.h
--- a/data/fft/FFTapi.h	Fri Nov 18 23:30:15 2016 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,103 +0,0 @@
-/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
-
-/*
-    Sonic Visualiser
-    An audio file viewer and annotation editor.
-    Centre for Digital Music, Queen Mary, University of London.
-    This file copyright 2006 Chris Cannam and QMUL.
-    
-    This program is free software; you can redistribute it and/or
-    modify it under the terms of the GNU General Public License as
-    published by the Free Software Foundation; either version 2 of the
-    License, or (at your option) any later version.  See the file
-    COPYING included with this distribution for more information.
-*/
-
-#ifndef _FFT_API_H_
-#define _FFT_API_H_
-
-#ifdef HAVE_FFTW3F
-
-#include <fftw3.h>
-
-#define fftf_complex fftwf_complex
-#define fftf_malloc fftwf_malloc
-#define fftf_free fftwf_free
-#define fftf_plan fftwf_plan
-#define fftf_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d
-#define fftf_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d
-#define fftf_execute fftwf_execute
-#define fftf_destroy_plan fftwf_destroy_plan
-
-#else
-
-// Provide a fallback FFT implementation if FFTW3f is not available.
-
-typedef float fftf_complex[2];
-#define fftf_malloc malloc
-#define fftf_free free
-
-struct fftf_plan_;
-typedef fftf_plan_ *fftf_plan;
-
-fftf_plan fftf_plan_dft_r2c_1d(int n, float *in, fftf_complex *out, unsigned);
-fftf_plan fftf_plan_dft_c2r_1d(int n, fftf_complex *in, float *out, unsigned);
-void fftf_execute(const fftf_plan p);
-void fftf_destroy_plan(fftf_plan p);
-
-#define FFTW_ESTIMATE 0
-#define FFTW_MEASURE 0
-
-#endif
-
-#include <vector>
-#include <complex>
-#include <mutex>
-
-class FFTForward // with fft shift but not window
-{
-    static std::mutex m_mutex;
-    
-public:
-    FFTForward(int size) :
-        m_size(size)
-    {
-        std::lock_guard<std::mutex> lock(m_mutex);
-        m_input = (float *)fftf_malloc(size * sizeof(float));
-        m_output = (fftf_complex *)fftf_malloc((size/2 + 1) * sizeof(fftf_complex));
-        m_plan = fftf_plan_dft_r2c_1d(size, m_input, m_output, FFTW_ESTIMATE);
-    }
-
-    ~FFTForward() {
-        std::lock_guard<std::mutex> lock(m_mutex);
-        fftf_destroy_plan(m_plan);
-        fftf_free(m_input);
-        fftf_free(m_output);
-    }
-
-    std::vector<std::complex<float> > process(const std::vector<float> &in) const {
-        const int hs = m_size/2;
-        for (int i = 0; i < hs; ++i) {
-            m_input[i] = in[i + hs];
-        }
-        for (int i = 0; i < hs; ++i) {
-            m_input[i + hs] = in[i];
-        }
-        fftf_execute(m_plan);
-        std::vector<std::complex<float> > result;
-        result.reserve(hs + 1);
-        for (int i = 0; i <= hs; ++i) {
-            result.push_back({ m_output[i][0], m_output[i][1] });
-        }
-        return result;
-    }
-
-private:
-    int m_size;
-    float *m_input;
-    fftf_complex *m_output;
-    fftf_plan m_plan;
-};
-
-#endif
-
diff -r f50c0bbe9096 -r 604f369c247a data/model/FFTModel.cpp
--- a/data/model/FFTModel.cpp	Fri Nov 18 23:30:15 2016 +0000
+++ b/data/model/FFTModel.cpp	Mon Nov 21 15:33:03 2016 +0000
@@ -52,6 +52,8 @@
         throw invalid_argument("FFTModel window size must be at least FFT size");
     }
 
+    m_fft.initFloat();
+
     connect(model, SIGNAL(modelChanged()), this, SIGNAL(modelChanged()));
     connect(model, SIGNAL(modelChangedWithin(sv_frame_t, sv_frame_t)),
             this, SIGNAL(modelChangedWithin(sv_frame_t, sv_frame_t)));
@@ -223,6 +225,8 @@
         return m_savedData.data;
     }
 
+    Profiler profiler("FFTModel::getSourceData (cache miss)");
+    
     if (range.first < m_savedData.range.second &&
         range.first >= m_savedData.range.first &&
         range.second > m_savedData.range.second) {
@@ -309,7 +313,12 @@
     
     auto samples = getSourceSamples(n);
     m_windower.cut(samples.data());
-    auto col = m_fft.process(samples);
+    breakfastquay::v_fftshift(samples.data(), m_fftSize);
+
+    vector<complex<float>> col(m_fftSize/2 + 1);
+    
+    m_fft.forwardInterleaved(samples.data(),
+                             reinterpret_cast<float *>(col.data()));
 
     SavedColumn sc { n, col };
     if (m_cached.size() >= m_cacheSize) {
diff -r f50c0bbe9096 -r 604f369c247a data/model/FFTModel.h
--- a/data/model/FFTModel.h	Fri Nov 18 23:30:15 2016 +0000
+++ b/data/model/FFTModel.h	Mon Nov 21 15:33:03 2016 +0000
@@ -21,7 +21,7 @@
 
 #include "base/Window.h"
 
-#include "data/fft/FFTapi.h"
+#include <bqfft/FFT.h>
 
 #include <set>
 #include <vector>
@@ -153,7 +153,7 @@
     int m_windowIncrement;
     int m_fftSize;
     Window<float> m_windower;
-    FFTForward m_fft;
+    mutable breakfastquay::FFT m_fft;
     
     int getPeakPickWindowSize(PeakPickType type, sv_samplerate_t sampleRate,
                               int bin, float &percentile) const;
diff -r f50c0bbe9096 -r 604f369c247a files.pri
--- a/files.pri	Fri Nov 18 23:30:15 2016 +0000
+++ b/files.pri	Mon Nov 21 15:33:03 2016 +0000
@@ -41,7 +41,6 @@
            base/Window.h \
            base/XmlExportable.h \
            base/ZoomConstraint.h \
-	   data/fft/FFTapi.h \
            data/fileio/AudioFileReader.h \
            data/fileio/AudioFileReaderFactory.h \
            data/fileio/AudioFileSizeEstimator.h \
@@ -174,7 +173,6 @@
            base/UnitDatabase.cpp \
            base/ViewManagerBase.cpp \
            base/XmlExportable.cpp \
-	   data/fft/FFTapi.cpp \
            data/fileio/AudioFileReader.cpp \
            data/fileio/AudioFileReaderFactory.cpp \
            data/fileio/AudioFileSizeEstimator.cpp \