Mercurial > hg > svcore
changeset 1270:bac86d3fc6c9 bqfft
Branch to use bqfft code and remove FFTapi
author | Chris Cannam |
---|---|
date | Mon, 21 Nov 2016 15:05:57 +0000 (2016-11-21) |
parents | f50c0bbe9096 |
children | 604f369c247a |
files | base/Window.h data/fft/FFTapi.cpp data/fft/FFTapi.h data/model/FFTModel.cpp data/model/FFTModel.h files.pri |
diffstat | 6 files changed, 47 insertions(+), 358 deletions(-) [+] |
line wrap: on
line diff
--- a/base/Window.h Fri Nov 18 23:30:15 2016 +0000 +++ b/base/Window.h Mon Nov 21 15:05:57 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];
--- 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
--- 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 -
--- a/data/model/FFTModel.cpp Fri Nov 18 23:30:15 2016 +0000 +++ b/data/model/FFTModel.cpp Mon Nov 21 15:05:57 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) {
--- a/data/model/FFTModel.h Fri Nov 18 23:30:15 2016 +0000 +++ b/data/model/FFTModel.h Mon Nov 21 15:05:57 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;
--- a/files.pri Fri Nov 18 23:30:15 2016 +0000 +++ b/files.pri Mon Nov 21 15:05:57 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 \