Reorganise into faster (R2) and finer (R3)

src/common/Allocators.cpp

@@ -0,0 +1,69 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#include "Allocators.h"
+
+#ifdef HAVE_IPP
+#include <ipps.h>
+#endif
+
+#include <iostream>
+using std::cerr;
+using std::endl;
+
+namespace RubberBand {
+
+#ifdef HAVE_IPP
+
+template <>
+float *allocate(size_t count)
+{
+    float *ptr = ippsMalloc_32f(count);
+    if (!ptr) throw (std::bad_alloc());
+    return ptr;
+}
+
+template <>
+double *allocate(size_t count)
+{
+    double *ptr = ippsMalloc_64f(count);
+    if (!ptr) throw (std::bad_alloc());
+    return ptr;
+}
+
+template <>
+void deallocate(float *ptr)
+{
+    if (ptr) ippsFree((void *)ptr);
+}
+
+template <>
+void deallocate(double *ptr)
+{
+    if (ptr) ippsFree((void *)ptr);
+}
+
+#endif
+
+}
+

src/common/Allocators.h

@@ -0,0 +1,407 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_ALLOCATORS_H
+#define RUBBERBAND_ALLOCATORS_H
+
+#include "VectorOps.h"
+
+#include <new> // for std::bad_alloc
+#include <stdlib.h>
+
+#include <stdexcept>
+
+#ifndef HAVE_POSIX_MEMALIGN
+#ifndef _WIN32
+#ifndef __APPLE__
+#ifndef LACK_POSIX_MEMALIGN
+#define HAVE_POSIX_MEMALIGN
+#endif
+#endif
+#endif
+#endif
+
+#ifndef MALLOC_IS_ALIGNED
+#ifndef MALLOC_IS_NOT_ALIGNED
+#ifdef __APPLE__
+#define MALLOC_IS_ALIGNED
+#endif
+#endif
+#endif
+
+#ifndef HAVE__ALIGNED_MALLOC
+#ifndef LACK__ALIGNED_MALLOC
+#ifdef _WIN32
+#define HAVE__ALIGNED_MALLOC
+#endif
+#endif
+#endif
+
+#ifdef HAVE_POSIX_MEMALIGN
+#include <sys/mman.h>
+#include <errno.h>
+#endif
+
+#ifdef LACK_BAD_ALLOC
+namespace std { struct bad_alloc { }; }
+#endif
+
+namespace RubberBand {
+
+template <typename T>
+T *allocate(size_t count)
+{
+    void *ptr = 0;
+
+    // We'd like to check HAVE_IPP first and, if it's defined, call
+    // ippsMalloc_8u(count * sizeof(T)). But that isn't a general
+    // replacement for malloc() because it only takes an int
+    // argument. So we save it for the specialisations of
+    // allocate<float> and allocate<double> below, where we're more
+    // likely to get away with it.
+
+#ifdef MALLOC_IS_ALIGNED
+    ptr = malloc(count * sizeof(T));
+#else /* !MALLOC_IS_ALIGNED */
+
+    // That's the "sufficiently aligned" functions dealt with, the
+    // rest need a specific alignment provided to the call. 32-byte
+    // alignment is required for at least OpenMAX
+    static const int alignment = 32;
+
+#ifdef HAVE__ALIGNED_MALLOC
+    ptr = _aligned_malloc(count * sizeof(T), alignment);
+#else /* !HAVE__ALIGNED_MALLOC */
+
+#ifdef HAVE_POSIX_MEMALIGN
+    int rv = posix_memalign(&ptr, alignment, count * sizeof(T));
+    if (rv) {
+#ifndef NO_EXCEPTIONS
+        if (rv == EINVAL) {
+            throw "Internal error: invalid alignment";
+        } else {
+            throw std::bad_alloc();
+        }
+#else
+        abort();
+#endif
+    }
+#else /* !HAVE_POSIX_MEMALIGN */
+    
+#ifdef USE_OWN_ALIGNED_MALLOC
+#pragma message("Rolling own aligned malloc: this is unlikely to perform as well as the alternatives")
+
+    // Alignment must be a power of two, bigger than the pointer
+    // size. Stuff the actual malloc'd pointer in just before the
+    // returned value.  This is the least desirable way to do this --
+    // the other options below are all better
+    size_t allocd = count * sizeof(T) + alignment;
+    void *buf = malloc(allocd);
+    if (buf) {
+        char *adj = (char *)buf;
+        while ((unsigned long long)adj & (alignment-1)) --adj;
+        ptr = ((char *)adj) + alignment;
+        new (((void **)ptr)[-1]) (void *);
+        ((void **)ptr)[-1] = buf;
+    }
+
+#else /* !USE_OWN_ALIGNED_MALLOC */
+
+#error "No aligned malloc available: define MALLOC_IS_ALIGNED to use system malloc, HAVE_POSIX_MEMALIGN if posix_memalign is available, HAVE__ALIGNED_MALLOC if _aligned_malloc is available, or USE_OWN_ALIGNED_MALLOC to roll our own"
+
+#endif /* !USE_OWN_ALIGNED_MALLOC */
+#endif /* !HAVE_POSIX_MEMALIGN */
+#endif /* !HAVE__ALIGNED_MALLOC */
+#endif /* !MALLOC_IS_ALIGNED */
+
+    if (!ptr) {
+#ifndef NO_EXCEPTIONS
+        throw std::bad_alloc();
+#else
+        abort();
+#endif
+    }
+
+    T *typed_ptr = static_cast<T *>(ptr);
+    for (size_t i = 0; i < count; ++i) {
+        new (typed_ptr + i) T;
+    }
+    return typed_ptr;
+}
+
+#ifdef HAVE_IPP
+
+template <>
+float *allocate(size_t count);
+
+template <>
+double *allocate(size_t count);
+
+#endif
+	
+template <typename T>
+T *allocate_and_zero(size_t count)
+{
+    T *ptr = allocate<T>(count);
+    v_zero(ptr, count);
+    return ptr;
+}
+
+template <typename T>
+void deallocate(T *ptr)
+{
+    if (!ptr) return;
+    
+#ifdef MALLOC_IS_ALIGNED
+    free((void *)ptr);
+#else /* !MALLOC_IS_ALIGNED */
+
+#ifdef HAVE__ALIGNED_MALLOC
+    _aligned_free((void *)ptr);
+#else /* !HAVE__ALIGNED_MALLOC */
+
+#ifdef HAVE_POSIX_MEMALIGN
+    free((void *)ptr);
+#else /* !HAVE_POSIX_MEMALIGN */
+    
+#ifdef USE_OWN_ALIGNED_MALLOC
+    free(((void **)ptr)[-1]);
+#else /* !USE_OWN_ALIGNED_MALLOC */
+
+#error "No aligned malloc available: define MALLOC_IS_ALIGNED to use system malloc, HAVE_POSIX_MEMALIGN if posix_memalign is available, or USE_OWN_ALIGNED_MALLOC to roll our own"
+
+#endif /* !USE_OWN_ALIGNED_MALLOC */
+#endif /* !HAVE_POSIX_MEMALIGN */
+#endif /* !HAVE__ALIGNED_MALLOC */
+#endif /* !MALLOC_IS_ALIGNED */
+}
+
+#ifdef HAVE_IPP
+
+template <>
+void deallocate(float *);
+
+template <>
+void deallocate(double *);
+
+#endif
+
+/// Reallocate preserving contents but leaving additional memory uninitialised	
+template <typename T>
+T *reallocate(T *ptr, size_t oldcount, size_t count)
+{
+    T *newptr = allocate<T>(count);
+    if (oldcount && ptr) {
+        v_copy(newptr, ptr, oldcount < count ? oldcount : count);
+    }
+    if (ptr) deallocate<T>(ptr);
+    return newptr;
+}
+
+/// Reallocate, zeroing all contents
+template <typename T>
+T *reallocate_and_zero(T *ptr, size_t oldcount, size_t count)
+{
+    ptr = reallocate(ptr, oldcount, count);
+    v_zero(ptr, count);
+    return ptr;
+}
+	
+/// Reallocate preserving contents and zeroing any additional memory	
+template <typename T>
+T *reallocate_and_zero_extension(T *ptr, size_t oldcount, size_t count)
+{
+    ptr = reallocate(ptr, oldcount, count);
+    if (count > oldcount) v_zero(ptr + oldcount, count - oldcount);
+    return ptr;
+}
+
+template <typename T>
+T **allocate_channels(size_t channels, size_t count)
+{
+    T **ptr = allocate<T *>(channels);
+    for (size_t c = 0; c < channels; ++c) {
+        ptr[c] = allocate<T>(count);
+    }
+    return ptr;
+}
+	
+template <typename T>
+T **allocate_and_zero_channels(size_t channels, size_t count)
+{
+    T **ptr = allocate<T *>(channels);
+    for (size_t c = 0; c < channels; ++c) {
+        ptr[c] = allocate_and_zero<T>(count);
+    }
+    return ptr;
+}
+
+template <typename T>
+void deallocate_channels(T **ptr, size_t channels)
+{
+    if (!ptr) return;
+    for (size_t c = 0; c < channels; ++c) {
+        deallocate<T>(ptr[c]);
+    }
+    deallocate<T *>(ptr);
+}
+	
+template <typename T>
+T **reallocate_channels(T **ptr,
+                        size_t oldchannels, size_t oldcount,
+                        size_t channels, size_t count)
+{
+    T **newptr = allocate_channels<T>(channels, count);
+    if (oldcount && ptr) {
+        for (size_t c = 0; c < oldchannels && c < channels; ++c) {
+            for (size_t i = 0; i < oldcount && i < count; ++i) {
+                newptr[c][i] = ptr[c][i];
+            }
+        }
+    } 
+    if (ptr) deallocate_channels<T>(ptr, oldchannels);
+    return newptr;
+}
+	
+template <typename T>
+T **reallocate_and_zero_extend_channels(T **ptr,
+                                        size_t oldchannels, size_t oldcount,
+                                        size_t channels, size_t count)
+{
+    T **newptr = allocate_and_zero_channels<T>(channels, count);
+    if (oldcount && ptr) {
+        for (size_t c = 0; c < oldchannels && c < channels; ++c) {
+            for (size_t i = 0; i < oldcount && i < count; ++i) {
+                newptr[c][i] = ptr[c][i];
+            }
+        }
+    } 
+    if (ptr) deallocate_channels<T>(ptr, oldchannels);
+    return newptr;
+}
+
+/// RAII class to call deallocate() on destruction
+template <typename T>
+class Deallocator
+{
+public:
+    Deallocator(T *t) : m_t(t) { }
+    ~Deallocator() { deallocate<T>(m_t); }
+private:
+    T *m_t;
+};
+
+/** Allocator for use with STL classes, e.g. vector, to ensure
+ *  alignment.  Based on example code by Stephan T. Lavavej.
+ *
+ *  e.g. std::vector<float, StlAllocator<float> > v;
+ */
+template <typename T>
+class StlAllocator
+{
+public:
+    typedef T *pointer;
+    typedef const T *const_pointer;
+    typedef T &reference;
+    typedef const T &const_reference;
+    typedef T value_type;
+    typedef size_t size_type;
+    typedef ptrdiff_t difference_type;
+
+    StlAllocator() { }
+    StlAllocator(const StlAllocator&) { }
+    template <typename U> StlAllocator(const StlAllocator<U>&) { }
+    ~StlAllocator() { }
+
+    T *
+    allocate(const size_t n) const {
+        if (n == 0) return 0;
+        if (n > max_size()) {
+#ifndef NO_EXCEPTIONS
+            throw std::length_error("Size overflow in StlAllocator::allocate()");
+#else
+            abort();
+#endif
+        }
+        return RubberBand::allocate<T>(n);
+    }
+
+    void
+    deallocate(T *const p, const size_t) const {
+        RubberBand::deallocate(p);
+    }
+
+    template <typename U>
+    T *
+    allocate(const size_t n, const U *) const {
+        return allocate(n);
+    }
+
+    T *
+    address(T &r) const {
+        return &r;
+    }
+
+    const T *
+    address(const T &s) const {
+        return &s;
+    }
+
+    size_t
+    max_size() const {
+        return (static_cast<size_t>(0) - static_cast<size_t>(1)) / sizeof(T);
+    }
+
+    template <typename U> struct rebind {
+        typedef StlAllocator<U> other;
+    };
+    
+    bool
+    operator==(const StlAllocator &) const {
+        return true;
+    }
+
+    bool
+    operator!=(const StlAllocator &) const {
+        return false;
+    }
+
+    void
+    construct(T *const p, const T &t) const {
+        void *const pv = static_cast<void *>(p);
+        new (pv) T(t);
+    }
+
+    void
+    destroy(T *const p) const {
+        p->~T();
+    }
+
+private:
+    StlAllocator& operator=(const StlAllocator&);
+};
+
+}
+
+#endif
+

src/common/BQResampler.cpp

@@ -0,0 +1,649 @@
+//* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#include "BQResampler.h"
+
+#include <cmath>
+
+#include <iostream>
+#include <algorithm>
+
+#include "Allocators.h"
+#include "VectorOps.h"
+
+#define BQ_R__ R__
+
+using std::vector;
+using std::cerr;
+using std::endl;
+using std::min;
+using std::max;
+
+namespace RubberBand {
+
+BQResampler::BQResampler(Parameters parameters, int channels) :
+    m_qparams(parameters.quality),
+    m_dynamism(parameters.dynamism),
+    m_ratio_change(parameters.ratioChange),
+    m_debug_level(parameters.debugLevel),
+    m_initial_rate(parameters.referenceSampleRate),
+    m_channels(channels),
+    m_fade_count(0),
+    m_initialised(false)
+{
+    if (m_debug_level > 0) {
+        cerr << "BQResampler::BQResampler: "
+             << (m_dynamism == RatioOftenChanging ? "often-changing" : "mostly-fixed")
+             << ", "
+             << (m_ratio_change == SmoothRatioChange ? "smooth" : "sudden")
+             << " ratio changes, ref " << m_initial_rate << " Hz" << endl;
+    }
+    
+    if (m_dynamism == RatioOftenChanging) {
+        m_proto_length = m_qparams.proto_p * m_qparams.p_multiple + 1;
+        if (m_debug_level > 0) {
+            cerr << "BQResampler: creating prototype filter of length "
+                 << m_proto_length << endl;
+        }
+        m_prototype = make_filter(m_proto_length, m_qparams.proto_p);
+        m_prototype.push_back(0.0); // interpolate without fear
+    }
+
+    int phase_reserve = 2 * int(round(m_initial_rate));
+    int buffer_reserve = 1000 * m_channels;
+    m_state_a.phase_info.reserve(phase_reserve);
+    m_state_a.buffer.reserve(buffer_reserve);
+
+    if (m_dynamism == RatioOftenChanging) {
+        m_state_b.phase_info.reserve(phase_reserve);
+        m_state_b.buffer.reserve(buffer_reserve);
+    }
+
+    m_s = &m_state_a;
+    m_fade = &m_state_b;
+}
+
+BQResampler::BQResampler(const BQResampler &other) :
+    m_qparams(other.m_qparams),
+    m_dynamism(other.m_dynamism),
+    m_ratio_change(other.m_ratio_change),
+    m_debug_level(other.m_debug_level),
+    m_initial_rate(other.m_initial_rate),
+    m_channels(other.m_channels),
+    m_state_a(other.m_state_a),
+    m_state_b(other.m_state_b),
+    m_fade_count(other.m_fade_count),
+    m_prototype(other.m_prototype),
+    m_proto_length(other.m_proto_length),
+    m_initialised(other.m_initialised)
+{
+    if (other.m_s == &(other.m_state_a)) {
+        m_s = &m_state_a;
+        m_fade = &m_state_b;
+    } else {
+        m_s = &m_state_b;
+        m_fade = &m_state_a;
+    }
+}
+
+void
+BQResampler::reset()
+{
+    m_initialised = false;
+    m_fade_count = 0;
+}
+
+BQResampler::QualityParams::QualityParams(Quality q)
+{
+    switch (q) {
+    case Fastest:
+        p_multiple = 12;
+        proto_p = 160;
+        k_snr = 70.0;
+        k_transition = 0.2;
+        cut = 0.9;
+        break;
+    case FastestTolerable:
+        p_multiple = 62;
+        proto_p = 160;
+        k_snr = 90.0;
+        k_transition = 0.05;
+        cut = 0.975;
+        break;
+    case Best:
+        p_multiple = 122;
+        proto_p = 800;
+        k_snr = 100.0;
+        k_transition = 0.01;
+        cut = 0.995;
+        break;
+    }
+}
+    
+int
+BQResampler::resampleInterleaved(float *const out,
+                                 int outspace,
+                                 const float *const in,
+                                 int incount,
+                                 double ratio,
+                                 bool final)
+{
+    int fade_length = int(round(m_initial_rate / 1000.0));
+    if (fade_length < 6) {
+        fade_length = 6;
+    }
+    int max_fade = min(outspace, int(floor(incount * ratio))) / 2;
+    if (fade_length > max_fade) {
+        fade_length = max_fade;
+    }
+        
+    if (!m_initialised) {
+        state_for_ratio(*m_s, ratio, *m_fade);
+        m_initialised = true;
+    } else if (ratio != m_s->parameters.ratio) {
+        state *tmp = m_fade;
+        m_fade = m_s;
+        m_s = tmp;
+        state_for_ratio(*m_s, ratio, *m_fade);
+        if (m_ratio_change == SmoothRatioChange) {
+            if (m_debug_level > 0) {
+                cerr << "BQResampler: ratio changed, beginning fade of length "
+                     << fade_length << endl;
+            }
+            m_fade_count = fade_length;
+        }
+    }
+
+    int i = 0, o = 0;
+    int bufsize = m_s->buffer.size();
+
+    int incount_samples = incount * m_channels;
+    int outspace_samples = outspace * m_channels;
+    
+    while (o < outspace_samples) {
+        while (i < incount_samples && m_s->fill < bufsize) {
+            m_s->buffer[m_s->fill++] = in[i++];
+        }
+        if (m_s->fill == bufsize) {
+            out[o++] = reconstruct_one(m_s);
+        } else if (final && m_s->fill > m_s->centre) {
+            out[o++] = reconstruct_one(m_s);
+        } else if (final && m_s->fill == m_s->centre &&
+                   m_s->current_phase != m_s->initial_phase) {
+            out[o++] = reconstruct_one(m_s);
+        } else {
+            break;
+        }
+    }
+
+    int fbufsize = m_fade->buffer.size();
+    int fi = 0, fo = 0;
+    while (fo < o && m_fade_count > 0) {
+        while (fi < incount_samples && m_fade->fill < fbufsize) {
+            m_fade->buffer[m_fade->fill++] = in[fi++];
+        }
+        if (m_fade->fill == fbufsize) {
+            double r = reconstruct_one(m_fade);
+            double fadeWith = out[fo];
+            double extent = double(m_fade_count - 1) / double(fade_length);
+            double mixture = 0.5 * (1.0 - cos(M_PI * extent));
+            double mixed = r * mixture + fadeWith * (1.0 - mixture);
+            out[fo] = mixed;
+            ++fo;
+            if (m_fade->current_channel == 0) {
+                --m_fade_count;
+            }
+        } else {
+            break;
+        }
+    }
+        
+    return o / m_channels;
+}
+
+double
+BQResampler::getEffectiveRatio(double ratio) const {
+    if (m_initialised && ratio == m_s->parameters.ratio) {
+        return m_s->parameters.effective;
+    } else {
+        return pick_params(ratio).effective;
+    }
+}
+    
+int
+BQResampler::gcd(int a, int b) const
+{
+    int c = a % b;
+    if (c == 0) return b;
+    else return gcd(b, c);
+}
+
+double
+BQResampler::bessel0(double x) const
+{
+    static double facsquared[] = {
+        0.0, 1.0, 4.0, 36.0,
+        576.0, 14400.0, 518400.0, 25401600.0,
+        1625702400.0, 131681894400.0, 1.316818944E13, 1.59335092224E15,
+        2.29442532803E17, 3.87757880436E19, 7.60005445655E21,
+        1.71001225272E24, 4.37763136697E26, 1.26513546506E29,
+        4.09903890678E31, 1.47975304535E34
+    };
+    static int nterms = sizeof(facsquared) / sizeof(facsquared[0]);
+    double b = 1.0;
+    for (int n = 1; n < nterms; ++n) {
+        double ff = facsquared[n];
+        double term = pow(x / 2.0, n * 2.0) / ff;
+        b += term;
+    }
+    return b;
+}
+
+vector<double>
+BQResampler::kaiser(double beta, int len) const
+{
+    double denominator = bessel0(beta);
+    int half = (len % 2 == 0 ? len/2 : (len+1)/2);
+    vector<double> v(len, 0.0);
+    for (int n = 0; n < half; ++n) {
+        double k = (2.0 * n) / (len-1) - 1.0;
+        v[n] = bessel0 (beta * sqrt(1.0 - k*k)) / denominator;
+    }
+    for (int n = half; n < len; ++n) {
+        v[n] = v[len-1 - n];
+    }
+    return v;
+}
+
+void
+BQResampler::kaiser_params(double attenuation,
+                           double transition,
+                           double &beta,
+                           int &len) const
+{
+    if (attenuation > 21.0) {
+        len = 1 + int(ceil((attenuation - 7.95) / (2.285 * transition)));
+    } else {
+        len = 1 + int(ceil(5.79 / transition));
+    }
+    beta = 0.0;
+    if (attenuation > 50.0) {
+        beta = 0.1102 * (attenuation - 8.7);
+    } else if (attenuation > 21.0) {
+        beta = 0.5842 * (pow (attenuation - 21.0, 0.4)) +
+            0.07886 * (attenuation - 21.0);
+    }
+}
+
+vector<double>
+BQResampler::kaiser_for(double attenuation,
+                        double transition,
+                        int minlen,
+                        int maxlen) const
+{
+    double beta;
+    int m;
+    kaiser_params(attenuation, transition, beta, m);
+    int mb = m;
+    if (maxlen > 0 && mb > maxlen - 1) {
+        mb = maxlen - 1;
+    } else if (minlen > 0 && mb < minlen) {
+        mb = minlen;
+    }
+    if (mb % 2 == 0) ++mb;
+    if (m_debug_level > 0) {
+        cerr << "BQResampler: window attenuation " << attenuation
+             << ", transition " << transition
+             << " -> length " << m << " adjusted to " << mb
+             << ", beta " << beta << endl;
+    }
+    return kaiser(beta, mb);
+}
+    
+void
+BQResampler::sinc_multiply(double peak_to_zero, vector<double> &buf) const
+{
+    int len = int(buf.size());
+    if (len < 2) return;
+
+    int left = len / 2;
+    int right = (len + 1) / 2;
+    double m = M_PI / peak_to_zero;
+
+    for (int i = 1; i <= right; ++i) {
+        double x = i * m;
+        double sinc = sin(x) / x;
+        if (i <= left) {
+            buf[left - i] *= sinc;
+        }
+        if (i < right) {
+            buf[i + left] *= sinc;
+        }
+    }
+}
+
+BQResampler::params
+BQResampler::fill_params(double ratio, double numd, double denomd) const
+{
+    int num = int(round(numd));
+    int denom = int(round(denomd));
+    params p;
+    int g = gcd (num, denom);
+    p.ratio = ratio;
+    p.numerator = num / g;
+    p.denominator = denom / g;
+    p.effective = double(p.numerator) / double(p.denominator);
+    p.peak_to_zero = max(p.denominator, p.numerator);
+    p.peak_to_zero /= m_qparams.cut;
+    p.scale = double(p.numerator) / double(p.peak_to_zero);
+
+    if (m_debug_level > 0) {
+        cerr << "BQResampler: ratio " << p.ratio
+             << " -> fraction " << p.numerator << "/" << p.denominator
+             << " with error " << p.effective - p.ratio
+             << endl;
+        cerr << "BQResampler: peak-to-zero " << p.peak_to_zero
+             << ", scale " << p.scale
+             << endl;
+    }
+    
+    return p;
+}
+    
+BQResampler::params
+BQResampler::pick_params(double ratio) const
+{
+    // Farey algorithm, see
+    // https://www.johndcook.com/blog/2010/10/20/best-rational-approximation/
+    int max_denom = 192000;
+    double a = 0.0, b = 1.0, c = 1.0, d = 0.0;
+    double pa = a, pb = b, pc = c, pd = d;
+    double eps = 1e-9;
+    while (b <= max_denom && d <= max_denom) {
+        double mediant = (a + c) / (b + d);
+        if (fabs(ratio - mediant) < eps) {
+            if (b + d <= max_denom) {
+                return fill_params(ratio, a + c, b + d);
+            } else if (d > b) {
+                return fill_params(ratio, c, d);
+            } else {
+                return fill_params(ratio, a, b);
+            }
+        }
+        if (ratio > mediant) {
+            pa = a; pb = b;
+            a += c; b += d;
+        } else {
+            pc = c; pd = d;
+            c += a; d += b;
+        }
+    }
+    if (fabs(ratio - (pc / pd)) < fabs(ratio - (pa / pb))) {
+        return fill_params(ratio, pc, pd);
+    } else {
+        return fill_params(ratio, pa, pb);
+    }
+}
+
+void
+BQResampler::phase_data_for(vector<BQResampler::phase_rec> &target_phase_data,
+                            floatbuf &target_phase_sorted_filter,
+                            int filter_length,
+                            const vector<double> *filter,
+                            int initial_phase,
+                            int input_spacing,
+                            int output_spacing) const
+{
+    target_phase_data.clear();
+    target_phase_data.reserve(input_spacing);
+        
+    for (int p = 0; p < input_spacing; ++p) {
+        int next_phase = p - output_spacing;
+        while (next_phase < 0) next_phase += input_spacing;
+        next_phase %= input_spacing;
+        double dspace = double(input_spacing);
+        int zip_length = int(ceil(double(filter_length - p) / dspace));
+        int drop = int(ceil(double(max(0, output_spacing - p)) / dspace));
+        phase_rec phase;
+        phase.next_phase = next_phase;
+        phase.drop = drop;
+        phase.length = zip_length;
+        phase.start_index = 0; // we fill this in below if needed
+        target_phase_data.push_back(phase);
+    }
+
+    if (m_dynamism == RatioMostlyFixed) {
+        if (!filter) throw std::logic_error("filter required at phase_data_for in RatioMostlyFixed mode");
+        target_phase_sorted_filter.clear();
+        target_phase_sorted_filter.reserve(filter_length);
+        for (int p = initial_phase; ; ) {
+            phase_rec &phase = target_phase_data[p];
+            phase.start_index = target_phase_sorted_filter.size();
+            for (int i = 0; i < phase.length; ++i) {
+                target_phase_sorted_filter.push_back
+                    ((*filter)[i * input_spacing + p]);
+            }
+            p = phase.next_phase;
+            if (p == initial_phase) {
+                break;
+            }
+        }
+    }
+}
+
+vector<double>
+BQResampler::make_filter(int filter_length, double peak_to_zero) const
+{
+    vector<double> filter;
+    filter.reserve(filter_length);
+
+    vector<double> kaiser = kaiser_for(m_qparams.k_snr, m_qparams.k_transition,
+                                       1, filter_length);
+    int k_length = kaiser.size();
+
+    if (k_length == filter_length) {
+        sinc_multiply(peak_to_zero, kaiser);
+        return kaiser;
+    } else {
+        kaiser.push_back(0.0);
+        double m = double(k_length - 1) / double(filter_length - 1);
+        for (int i = 0; i < filter_length; ++i) {
+            double ix = i * m;
+            int iix = int(floor(ix));
+            double remainder = ix - iix;
+            double value = 0.0;
+            value += kaiser[iix] * (1.0 - remainder);
+            value += kaiser[iix+1] * remainder;
+            filter.push_back(value);
+        }
+        sinc_multiply(peak_to_zero, filter);
+        return filter;
+    }
+}
+
+void
+BQResampler::state_for_ratio(BQResampler::state &target_state,
+                             double ratio,
+                             const BQResampler::state &BQ_R__ prev_state) const
+{
+    params parameters = pick_params(ratio);
+    target_state.parameters = parameters;
+    
+    target_state.filter_length =
+        int(parameters.peak_to_zero * m_qparams.p_multiple + 1);
+
+    if (target_state.filter_length % 2 == 0) {
+        ++target_state.filter_length;
+    }
+
+    int half_length = target_state.filter_length / 2; // nb length is odd
+    int input_spacing = parameters.numerator;
+    int initial_phase = half_length % input_spacing;
+
+    target_state.initial_phase = initial_phase;
+    target_state.current_phase = initial_phase;
+
+    if (m_dynamism == RatioMostlyFixed) {
+        
+        if (m_debug_level > 0) {
+            cerr << "BQResampler: creating filter of length "
+                 << target_state.filter_length << endl;
+        }
+
+        vector<double> filter =
+            make_filter(target_state.filter_length, parameters.peak_to_zero);
+
+        phase_data_for(target_state.phase_info,
+                       target_state.phase_sorted_filter,
+                       target_state.filter_length, &filter,
+                       target_state.initial_phase,
+                       input_spacing,
+                       parameters.denominator);
+    } else {
+        phase_data_for(target_state.phase_info,
+                       target_state.phase_sorted_filter,
+                       target_state.filter_length, 0,
+                       target_state.initial_phase,
+                       input_spacing,
+                       parameters.denominator);
+    }
+
+    int buffer_left = half_length / input_spacing;
+    int buffer_right = buffer_left + 1;
+
+    int buffer_length = buffer_left + buffer_right;
+    buffer_length = max(buffer_length,
+                        int(prev_state.buffer.size() / m_channels));
+
+    target_state.centre = buffer_length / 2;
+    target_state.left = target_state.centre - buffer_left;
+    target_state.fill = target_state.centre;
+
+    buffer_length *= m_channels;
+    target_state.centre *= m_channels;
+    target_state.left *= m_channels;
+    target_state.fill *= m_channels;
+    
+    int n_phases = int(target_state.phase_info.size());
+
+    if (m_debug_level > 0) {
+        cerr << "BQResampler: " << m_channels << " channel(s) interleaved"
+             << ", buffer left " << buffer_left
+             << ", right " << buffer_right
+             << ", total " << buffer_length << endl;
+    
+        cerr << "BQResampler: input spacing " << input_spacing
+             << ", output spacing " << parameters.denominator
+             << ", initial phase " << initial_phase
+             << " of " << n_phases << endl;
+    }
+
+    if (prev_state.buffer.size() > 0) {
+        if (int(prev_state.buffer.size()) == buffer_length) {
+            target_state.buffer = prev_state.buffer;
+            target_state.fill = prev_state.fill;
+        } else {
+            target_state.buffer = floatbuf(buffer_length, 0.0);
+            for (int i = 0; i < prev_state.fill; ++i) {
+                int offset = i - prev_state.centre;
+                int new_ix = offset + target_state.centre;
+                if (new_ix >= 0 && new_ix < buffer_length) {
+                    target_state.buffer[new_ix] = prev_state.buffer[i];
+                    target_state.fill = new_ix + 1;
+                }
+            }
+        }
+
+        int phases_then = int(prev_state.phase_info.size());
+        double distance_through =
+            double(prev_state.current_phase) / double(phases_then);
+        target_state.current_phase = int(round(n_phases * distance_through));
+        if (target_state.current_phase >= n_phases) {
+            target_state.current_phase = n_phases - 1;
+        }
+    } else {
+        target_state.buffer = floatbuf(buffer_length, 0.0);
+    }
+}
+
+double
+BQResampler::reconstruct_one(state *s) const
+{
+    const phase_rec &pr = s->phase_info[s->current_phase];
+    int phase_length = pr.length;
+    double result = 0.0;
+
+    int dot_length =
+        min(phase_length,
+            (int(s->buffer.size()) - s->left) / m_channels);
+
+    if (m_dynamism == RatioMostlyFixed) {
+        int phase_start = pr.start_index;
+        if (m_channels == 1) {
+            result = v_multiply_and_sum
+                (s->phase_sorted_filter.data() + phase_start,
+                 s->buffer.data() + s->left,
+                 dot_length);
+        } else {
+            for (int i = 0; i < dot_length; ++i) {
+                result +=
+                    s->phase_sorted_filter[phase_start + i] *
+                    s->buffer[s->left + i * m_channels + s->current_channel];
+            }
+        }
+    } else {
+        double m = double(m_proto_length - 1) / double(s->filter_length - 1);
+        for (int i = 0; i < dot_length; ++i) {
+            double sample =
+                s->buffer[s->left + i * m_channels + s->current_channel];
+            int filter_index = i * s->parameters.numerator + s->current_phase;
+            double proto_index = m * filter_index;
+            int iix = int(floor(proto_index));
+            double remainder = proto_index - iix;
+            double filter_value = m_prototype[iix] * (1.0 - remainder);
+            filter_value += m_prototype[iix+1] * remainder;
+            result += filter_value * sample;
+        }
+    }
+
+    s->current_channel = (s->current_channel + 1) % m_channels;
+    
+    if (s->current_channel == 0) {
+
+        if (pr.drop > 0) {
+            int drop = pr.drop * m_channels;
+            v_move(s->buffer.data(), s->buffer.data() + drop,
+                   int(s->buffer.size()) - drop);
+            for (int i = 1; i <= drop; ++i) {
+                s->buffer[s->buffer.size() - i] = 0.0;
+            }
+            s->fill -= drop;
+        }
+
+        s->current_phase = pr.next_phase;
+    }
+    
+    return result * s->parameters.scale;
+}
+
+}

src/common/BQResampler.h

@@ -0,0 +1,167 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef BQ_BQRESAMPLER_H
+#define BQ_BQRESAMPLER_H
+
+#include <vector>
+
+#include "Allocators.h"
+#include "VectorOps.h"
+
+namespace RubberBand {
+
+class BQResampler
+{
+public:
+    enum Quality { Best, FastestTolerable, Fastest };
+    enum Dynamism { RatioOftenChanging, RatioMostlyFixed };
+    enum RatioChange { SmoothRatioChange, SuddenRatioChange };
+    
+    struct Parameters {
+        Quality quality;
+        Dynamism dynamism; 
+        RatioChange ratioChange;
+        double referenceSampleRate;
+        int debugLevel;
+
+        Parameters() :
+            quality(FastestTolerable),
+            dynamism(RatioMostlyFixed),
+            ratioChange(SmoothRatioChange),
+            referenceSampleRate(44100),
+            debugLevel(0) { }
+    };
+
+    BQResampler(Parameters parameters, int channels);
+    BQResampler(const BQResampler &);
+
+    int resampleInterleaved(float *const out, int outspace,
+                            const float *const in, int incount,
+                            double ratio, bool final);
+
+    double getEffectiveRatio(double ratio) const;
+    
+    void reset();
+
+private:
+    struct QualityParams {
+        int p_multiple;
+        int proto_p;
+        double k_snr;
+        double k_transition;
+        double cut;
+        QualityParams(Quality);
+    };
+
+    const QualityParams m_qparams;
+    const Dynamism m_dynamism;
+    const RatioChange m_ratio_change;
+    const int m_debug_level;
+    const double m_initial_rate;
+    const int m_channels;
+
+    struct params {
+        double ratio;
+        int numerator;
+        int denominator;
+        double effective;
+        double peak_to_zero;
+        double scale;
+        params() : ratio(1.0), numerator(1), denominator(1),
+                   effective(1.0), peak_to_zero(0), scale(1.0) { }
+    };
+
+    struct phase_rec {
+        int next_phase;
+        int length;
+        int start_index;
+        int drop;
+        phase_rec() : next_phase(0), length(0), start_index(0), drop(0) { }
+    };
+
+    typedef std::vector<float, RubberBand::StlAllocator<float> > floatbuf;
+    
+    struct state {
+        params parameters;
+        int initial_phase;
+        int current_phase;
+        int current_channel;
+        int filter_length;
+        std::vector<phase_rec> phase_info;
+        floatbuf phase_sorted_filter;
+        floatbuf buffer;
+        int left;
+        int centre;
+        int fill;
+        state() : initial_phase(0), current_phase(0), current_channel(0),
+                  filter_length(0), left(0), centre(0), fill(0) { }
+    };
+
+    state m_state_a;
+    state m_state_b;
+
+    state *m_s;        // points at either m_state_a or m_state_b
+    state *m_fade;     // whichever one m_s does not point to
+    
+    int m_fade_count;
+    
+    std::vector<double> m_prototype;
+    int m_proto_length;
+    bool m_initialised;
+
+    int gcd(int a, int b) const;
+    double bessel0(double x) const;
+    std::vector<double> kaiser(double beta, int len) const;
+    void kaiser_params(double attenuation, double transition,
+                       double &beta, int &len) const;
+    std::vector<double> kaiser_for(double attenuation, double transition,
+                                   int minlen, int maxlen) const;
+    void sinc_multiply(double peak_to_zero, std::vector<double> &buf) const;
+
+    params fill_params(double ratio, double numd, double denomd) const;
+    params pick_params(double ratio) const;
+
+    std::vector<double> make_filter(int filter_length,
+                                    double peak_to_zero) const;
+    
+    void phase_data_for(std::vector<phase_rec> &target_phase_data,
+                        floatbuf &target_phase_sorted_filter,
+                        int filter_length,
+                        const std::vector<double> *filter,
+                        int initial_phase,
+                        int input_spacing,
+                        int output_spacing) const;
+    
+    void state_for_ratio(state &target_state,
+                         double ratio,
+                         const state &R__ prev_state) const;
+    
+    double reconstruct_one(state *s) const;
+
+    BQResampler &operator=(const BQResampler &); // not provided
+};
+
+}
+
+#endif

src/common/FFT.h

@@ -0,0 +1,135 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_FFT_H
+#define RUBBERBAND_FFT_H
+
+#include "sysutils.h"
+
+#include <string>
+#include <set>
+
+namespace RubberBand {
+
+class FFTImpl;
+
+/**
+ * Provide the basic FFT computations we need, using one of a set of
+ * candidate FFT implementations (depending on compile flags).
+ *
+ * Implements real->complex FFTs of power-of-two sizes only.  Note
+ * that only the first half of the output signal is returned (the
+ * complex conjugates half is omitted), so the "complex" arrays need
+ * room for size/2+1 elements.
+ *
+ * The "interleaved" functions use the format sometimes called CCS --
+ * size/2+1 real+imaginary pairs.  So, the array elements at indices 1
+ * and size+1 will always be zero (since the signal is real).
+ * 
+ * All pointer arguments must point to valid data. A NullArgument
+ * exception is thrown if any argument is NULL.
+ *
+ * Neither forward nor inverse transform is scaled.
+ *
+ * This class is reentrant but not thread safe: use a separate
+ * instance per thread, or use a mutex.
+ */
+class FFT
+{
+public:
+    enum Exception {
+        NullArgument, InvalidSize, InvalidImplementation, InternalError
+    };
+
+    FFT(int size, int debugLevel = 0); // may throw InvalidSize
+    ~FFT();
+
+    int getSize() const;
+    
+    void forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut);
+    void forwardInterleaved(const double *R__ realIn, double *R__ complexOut);
+    void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut);
+    void forwardMagnitude(const double *R__ realIn, double *R__ magOut);
+
+    void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut);
+    void forwardInterleaved(const float *R__ realIn, float *R__ complexOut);
+    void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut);
+    void forwardMagnitude(const float *R__ realIn, float *R__ magOut);
+
+    void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut);
+    void inverseInterleaved(const double *R__ complexIn, double *R__ realOut);
+    void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut);
+    void inverseCepstral(const double *R__ magIn, double *R__ cepOut);
+
+    void inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut);
+    void inverseInterleaved(const float *R__ complexIn, float *R__ realOut);
+    void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut);
+    void inverseCepstral(const float *R__ magIn, float *R__ cepOut);
+
+    // Calling one or both of these is optional -- if neither is
+    // called, the first call to a forward or inverse method will call
+    // init().  You only need call these if you don't want to risk
+    // expensive allocations etc happening in forward or inverse.
+    void initFloat();
+    void initDouble();
+
+    enum Precision {
+        SinglePrecision = 0x1,
+        DoublePrecision = 0x2
+    };
+    typedef int Precisions;
+
+    /**
+     * Return the OR of all precisions supported by this
+     * implementation. All of the functions (float and double) are
+     * available regardless of the supported implementations, but they
+     * will be calculated at the proper precision only if it is
+     * available. (So float functions will be calculated using doubles
+     * and then truncated if single-precision is unavailable, and
+     * double functions will use single-precision arithmetic if double
+     * is unavailable.)
+     */
+    Precisions getSupportedPrecisions() const;
+
+    static std::set<std::string> getImplementations();
+    static std::string getDefaultImplementation();
+    static void setDefaultImplementation(std::string);
+
+#ifdef FFT_MEASUREMENT
+    static std::string tune();
+#endif
+
+protected:
+    FFTImpl *d;
+    static std::string m_implementation;
+    static void pickDefaultImplementation();
+    
+private:
+    FFT(const FFT &); // not provided
+    FFT &operator=(const FFT &); // not provided
+};
+
+}
+
+#endif
+

src/common/MovingMedian.h

@@ -0,0 +1,132 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_MOVING_MEDIAN_H
+#define RUBBERBAND_MOVING_MEDIAN_H
+
+#include "SampleFilter.h"
+#include "Allocators.h"
+
+#include <algorithm>
+#include <iostream>
+
+namespace RubberBand
+{
+
+template <typename T>
+class MovingMedian : public SampleFilter<T>
+{
+    typedef SampleFilter<T> P;
+
+public:
+    MovingMedian(int size, float percentile = 50.f) :
+        SampleFilter<T>(size),
+	m_frame(allocate_and_zero<T>(size)),
+	m_sorted(allocate_and_zero<T>(size)),
+	m_sortend(m_sorted + P::m_size - 1) {
+        setPercentile(percentile);
+    }
+
+    ~MovingMedian() { 
+	deallocate(m_frame);
+	deallocate(m_sorted);
+    }
+
+    void setPercentile(float p) {
+        m_index = int((P::m_size * p) / 100.f);
+        if (m_index >= P::m_size) m_index = P::m_size-1;
+        if (m_index < 0) m_index = 0;
+    }
+
+    void push(T value) {
+        if (value != value) {
+            std::cerr << "WARNING: MovingMedian: NaN encountered" << std::endl;
+            value = T();
+        }
+	drop(m_frame[0]);
+	v_move(m_frame, m_frame+1, P::m_size-1);
+	m_frame[P::m_size-1] = value;
+	put(value);
+    }
+
+    T get() const {
+	return m_sorted[m_index];
+    }
+
+    void reset() {
+	v_zero(m_frame, P::m_size);
+	v_zero(m_sorted, P::m_size);
+    }
+
+    // Convenience function that applies a given filter to an array
+    // in-place. Array must have length equal to getSize(). Modifies
+    // both the filter and the array.
+    //
+    static void filter(MovingMedian<T> &mm, T *v) {
+        int n = mm.getSize();
+        int lag = n / 2;
+        mm.reset();
+        for (int i = 0; i < lag; ++i) {
+            mm.push(v[i]);
+        }
+        for (int i = lag; i < n; ++i) {
+            mm.push(v[i]);
+            v[i-lag] = mm.get();
+        }
+        for (int i = n; i < n + lag; ++i) {
+            mm.push(T());
+            v[i-lag] = mm.get();
+        }
+    }
+    
+private:
+    T *const m_frame;
+    T *const m_sorted;
+    T *const m_sortend;
+    int m_index;
+
+    void put(T value) {
+	// precondition: m_sorted contains m_size-1 values, packed at start
+	// postcondition: m_sorted contains m_size values, one of which is value
+	T *index = std::lower_bound(m_sorted, m_sortend, value);
+	v_move(index + 1, index, m_sortend - index);
+	*index = value;
+    }
+
+    void drop(T value) {
+	// precondition: m_sorted contains m_size values, one of which is value
+	// postcondition: m_sorted contains m_size-1 values, packed at start
+	T *index = std::lower_bound(m_sorted, m_sortend + 1, value);
+	assert(*index == value);
+	v_move(index, index + 1, m_sortend - index);
+	*m_sortend = T(0);
+    }
+
+    MovingMedian(const MovingMedian &) =delete;
+    MovingMedian &operator=(const MovingMedian &) =delete;
+};
+
+}
+
+#endif
+

src/common/Profiler.cpp

@@ -0,0 +1,239 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#include "Profiler.h"
+
+#include "Thread.h"
+
+#include <algorithm>
+#include <set>
+#include <string>
+#include <map>
+
+#include <stdio.h>
+
+#ifdef _MSC_VER
+// Ugh --cc
+#define snprintf sprintf_s
+#endif
+
+namespace RubberBand {
+
+#ifndef NO_TIMING
+
+Profiler::ProfileMap
+Profiler::m_profiles;
+
+Profiler::WorstCallMap
+Profiler::m_worstCalls;
+
+static Mutex profileMutex;
+
+void
+Profiler::add(const char *id, float ms)
+{
+    profileMutex.lock();
+    
+    ProfileMap::iterator pmi = m_profiles.find(id);
+    if (pmi != m_profiles.end()) {
+        ++pmi->second.first;
+        pmi->second.second += ms;
+    } else {
+        m_profiles[id] = TimePair(1, ms);
+    }
+
+    WorstCallMap::iterator wci = m_worstCalls.find(id);
+    if (wci != m_worstCalls.end()) {
+        if (ms > wci->second) wci->second = ms;
+    } else {
+        m_worstCalls[id] = ms;
+    }
+
+    profileMutex.unlock();
+}
+
+void
+Profiler::dump()
+{
+    std::string report = getReport();
+    fprintf(stderr, "%s", report.c_str());
+}
+
+std::string
+Profiler::getReport()
+{
+    profileMutex.lock();
+    
+    static const int buflen = 256;
+    char buffer[buflen];
+    std::string report;
+
+#ifdef PROFILE_CLOCKS
+    snprintf(buffer, buflen, "Profiling points [CPU time]:\n");
+#else
+    snprintf(buffer, buflen, "Profiling points [Wall time]:\n");
+#endif
+    report += buffer;
+
+    typedef std::multimap<float, const char *> TimeRMap;
+    typedef std::multimap<int, const char *> IntRMap;
+    TimeRMap totmap, avgmap, worstmap;
+    IntRMap ncallmap;
+
+    for (ProfileMap::const_iterator i = m_profiles.begin();
+         i != m_profiles.end(); ++i) {
+        totmap.insert(TimeRMap::value_type(i->second.second, i->first));
+        avgmap.insert(TimeRMap::value_type(i->second.second /
+                                           i->second.first, i->first));
+        ncallmap.insert(IntRMap::value_type(i->second.first, i->first));
+    }
+
+    for (WorstCallMap::const_iterator i = m_worstCalls.begin();
+         i != m_worstCalls.end(); ++i) {
+        worstmap.insert(TimeRMap::value_type(i->second, i->first));
+    }
+
+    snprintf(buffer, buflen, "\nBy total:\n");
+    report += buffer;
+    for (TimeRMap::const_iterator i = totmap.end(); i != totmap.begin(); ) {
+        --i;
+        snprintf(buffer, buflen, "%-40s  %f ms\n", i->second, i->first);
+        report += buffer;
+    }
+
+    snprintf(buffer, buflen, "\nBy average:\n");
+    report += buffer;
+    for (TimeRMap::const_iterator i = avgmap.end(); i != avgmap.begin(); ) {
+        --i;
+        snprintf(buffer, buflen, "%-40s  %f ms\n", i->second, i->first);
+        report += buffer;
+    }
+
+    snprintf(buffer, buflen, "\nBy worst case:\n");
+    report += buffer;
+    for (TimeRMap::const_iterator i = worstmap.end(); i != worstmap.begin(); ) {
+        --i;
+        snprintf(buffer, buflen, "%-40s  %f ms\n", i->second, i->first);
+        report += buffer;
+    }
+
+    snprintf(buffer, buflen, "\nBy number of calls:\n");
+    report += buffer;
+    for (IntRMap::const_iterator i = ncallmap.end(); i != ncallmap.begin(); ) {
+        --i;
+        snprintf(buffer, buflen, "%-40s  %d\n", i->second, i->first);
+        report += buffer;
+    }
+
+    snprintf(buffer, buflen, "\nBy name:\n");
+    report += buffer;
+
+    typedef std::set<const char *, std::less<std::string> > StringSet;
+
+    StringSet profileNames;
+    for (ProfileMap::const_iterator i = m_profiles.begin();
+         i != m_profiles.end(); ++i) {
+        profileNames.insert(i->first);
+    }
+
+    for (StringSet::const_iterator i = profileNames.begin();
+         i != profileNames.end(); ++i) {
+
+        ProfileMap::const_iterator j = m_profiles.find(*i);
+        if (j == m_profiles.end()) continue;
+
+        const TimePair &pp(j->second);
+        snprintf(buffer, buflen, "%s(%d):\n", *i, pp.first);
+        report += buffer;
+        snprintf(buffer, buflen, "\tReal: \t%f ms      \t[%f ms total]\n",
+                (pp.second / pp.first),
+                (pp.second));
+        report += buffer;
+
+        WorstCallMap::const_iterator k = m_worstCalls.find(*i);
+        if (k == m_worstCalls.end()) continue;
+        
+        snprintf(buffer, buflen, "\tWorst:\t%f ms/call\n", k->second);
+        report += buffer;
+    }
+
+    profileMutex.unlock();
+    
+    return report;
+}
+
+Profiler::Profiler(const char* c) :
+    m_c(c),
+    m_ended(false)
+{
+#ifdef PROFILE_CLOCKS
+    m_start = clock();
+#else
+    (void)gettimeofday(&m_start, 0);
+#endif
+}
+
+Profiler::~Profiler()
+{
+    if (!m_ended) end();
+}
+
+void
+Profiler::end()
+{
+#ifdef PROFILE_CLOCKS
+    clock_t end = clock();
+    clock_t elapsed = end - m_start;
+    float ms = float((double(elapsed) / double(CLOCKS_PER_SEC)) * 1000.0);
+#else
+    struct timeval tv;
+    (void)gettimeofday(&tv, 0);
+
+    tv.tv_sec -= m_start.tv_sec;
+    if (tv.tv_usec < m_start.tv_usec) {
+        tv.tv_usec += 1000000;
+        tv.tv_sec -= 1;
+    }
+    tv.tv_usec -= m_start.tv_usec;
+    float ms = float((double(tv.tv_sec) + (double(tv.tv_usec) / 1000000.0)) * 1000.0);
+#endif
+
+    add(m_c, ms);
+
+    m_ended = true;
+}
+ 
+#else /* NO_TIMING */
+
+#ifndef NO_TIMING_COMPLETE_NOOP
+
+Profiler::Profiler(const char *) { }
+Profiler::~Profiler() { }
+void Profiler::end() { }
+void Profiler::dump() { }
+
+#endif
+
+#endif
+
+}

src/common/Profiler.h

@@ -0,0 +1,130 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_PROFILER_H
+#define RUBBERBAND_PROFILER_H
+
+//#define NO_TIMING 1
+//#define WANT_TIMING 1
+//#define PROFILE_CLOCKS 1
+
+// Define NO_TIMING or NDEBUG to switch off profilers
+#ifdef NDEBUG
+#define NO_TIMING 1
+#endif
+
+// But we always allow WANT_TIMING to switch them back on again
+#ifdef WANT_TIMING
+#undef NO_TIMING
+#endif
+
+#ifndef NO_TIMING
+#ifdef PROFILE_CLOCKS
+#include <time.h>
+#else
+#include "sysutils.h"
+#ifndef _WIN32
+#include <sys/time.h>
+#endif
+#endif
+#endif
+
+#ifndef NO_TIMING
+#include <map>
+#include <string>
+#endif
+
+namespace RubberBand {
+
+#ifndef NO_TIMING
+
+class Profiler
+{
+public:
+    Profiler(const char *name);
+    ~Profiler();
+
+    void end(); // same action as dtor
+
+    static void dump();
+
+    // Unlike the other functions, this is only defined if NO_TIMING
+    // is not set (because it uses std::string which is otherwise
+    // unused here). So, treat this as a tricksy internal function
+    // rather than an API call and guard any call to it appropriately.
+    static std::string getReport();
+
+protected:
+    const char* m_c;
+#ifdef PROFILE_CLOCKS
+    clock_t m_start;
+#else
+    struct timeval m_start;
+#endif
+    bool m_showOnDestruct;
+    bool m_ended;
+
+    typedef std::pair<int, float> TimePair;
+    typedef std::map<const char *, TimePair> ProfileMap;
+    typedef std::map<const char *, float> WorstCallMap;
+    static ProfileMap m_profiles;
+    static WorstCallMap m_worstCalls;
+    static void add(const char *, float);
+};
+
+#else
+
+#ifdef NO_TIMING_COMPLETE_NOOP
+
+// Fastest for release builds, but annoying because it can't be linked
+// with code built in debug mode (expecting non-inline functions), so
+// not preferred during development
+
+class Profiler
+{
+public:
+    Profiler(const char *) { }
+    ~Profiler() { }
+
+    void end() { }
+    static void dump() { }
+};
+
+#else
+
+class Profiler
+{
+public:
+    Profiler(const char *);
+    ~Profiler();
+
+    void end();
+    static void dump();
+};
+
+#endif
+#endif
+
+}
+
+#endif

src/common/Resampler.h

@@ -0,0 +1,179 @@
+//* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_RESAMPLER_H
+#define RUBBERBAND_RESAMPLER_H
+
+#include "sysutils.h"
+
+namespace RubberBand {
+
+class Resampler
+{
+public:
+    enum Quality { Best, FastestTolerable, Fastest };
+    enum Dynamism { RatioOftenChanging, RatioMostlyFixed };
+    enum RatioChange { SmoothRatioChange, SuddenRatioChange };
+    
+    enum Exception { ImplementationError };
+
+    struct Parameters {
+
+        /**
+         * Resampler filter quality level.
+         */
+        Quality quality;
+
+        /**
+         * Performance hint indicating whether the ratio is expected
+         * to change regularly or not. If not, more work may happen on
+         * ratio changes to reduce work when ratio is unchanged.
+         */
+        Dynamism dynamism; 
+
+        /**
+         * Hint indicating whether to smooth transitions, via filter
+         * interpolation or some such method, at ratio change
+         * boundaries, or whether to make a precise switch to the new
+         * ratio without regard to audible artifacts. The actual
+         * effect of this depends on the implementation in use.
+         */
+        RatioChange ratioChange;
+        
+        /** 
+         * Rate of expected input prior to resampling: may be used to
+         * determine the filter bandwidth for the quality setting. If
+         * you don't know what this will be, you can provide an
+         * arbitrary rate (such as the default) and the resampler will
+         * work fine, but quality may not be as designed.
+         */
+        double initialSampleRate;
+
+        /** 
+         * Bound on the maximum incount size that may be passed to the
+         * resample function before the resampler needs to reallocate
+         * its internal buffers. Default is zero, so that buffer
+         * allocation will happen on the first call and any subsequent
+         * call with a greater incount.
+         */
+        int maxBufferSize;
+
+        /**
+         * Debug output level, from 0 to 3. Controls the amount of
+         * debug information printed to stderr.
+         */
+        int debugLevel;
+
+        Parameters() :
+            quality(FastestTolerable),
+            dynamism(RatioMostlyFixed),
+            ratioChange(SmoothRatioChange),
+            initialSampleRate(44100),
+            maxBufferSize(0),
+            debugLevel(0) { }
+    };
+    
+    /**
+     * Construct a resampler to process the given number of channels,
+     * with the given quality level, initial sample rate, and other
+     * parameters.
+     */
+    Resampler(Parameters parameters, int channels);
+    
+    ~Resampler();
+
+    /**
+     * Resample the given multi-channel buffers, where incount is the
+     * number of frames in the input buffers and outspace is the space
+     * available in the output buffers. Generally you want outspace to
+     * be at least ceil(incount * ratio).
+     *
+     * Returns the number of frames written to the output
+     * buffers. This may be smaller than outspace even where the ratio
+     * suggests otherwise, particularly at the start of processing
+     * where there may be a filter tail to allow for.
+     */
+#ifdef __GNUC__
+    __attribute__((warn_unused_result))
+#endif
+    int resample(float *const R__ *const R__ out,
+                 int outspace,
+                 const float *const R__ *const R__ in,
+                 int incount,
+                 double ratio,
+                 bool final = false);
+
+    /**
+     * Resample the given interleaved buffer, where incount is the
+     * number of frames in the input buffer (i.e. it has incount *
+     * getChannelCount() samples) and outspace is the space available
+     * in frames in the output buffer (i.e. it has space for at least
+     * outspace * getChannelCount() samples). Generally you want
+     * outspace to be at least ceil(incount * ratio).
+     *
+     * Returns the number of frames written to the output buffer. This
+     * may be smaller than outspace even where the ratio suggests
+     * otherwise, particularly at the start of processing where there
+     * may be a filter tail to allow for.
+     */
+#ifdef __GNUC__
+    __attribute__((warn_unused_result))
+#endif
+    int resampleInterleaved(float *const R__ out,
+                            int outspace,
+                            const float *const R__ in,
+                            int incount,
+                            double ratio,
+                            bool final = false);
+
+    /**
+     * Return the channel count provided on construction.
+     */
+    int getChannelCount() const;
+
+    /**
+     * Return the ratio that will be actually used when the given
+     * ratio is requested. For example, if the resampler internally
+     * uses a rational approximation of the given ratio, this will
+     * return the closest double to that approximation. Not all
+     * implementations support this; an implementation that does not
+     * will just return the given ratio.
+     */
+    double getEffectiveRatio(double ratio) const;
+    
+    /**
+     * Reset the internal processing state so that the next call
+     * behaves as if the resampler had just been constructed.
+     */
+    void reset();
+
+    class Impl;
+
+protected:
+    Impl *d;
+    int m_method;
+};
+
+}
+
+#endif

src/common/RingBuffer.h

@@ -0,0 +1,526 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_RINGBUFFER_H
+#define RUBBERBAND_RINGBUFFER_H
+
+#include <sys/types.h>
+
+//#define DEBUG_RINGBUFFER 1
+
+#include "sysutils.h"
+#include "Allocators.h"
+
+#include <iostream>
+
+#include <atomic>
+
+namespace RubberBand {
+
+/**
+ * RingBuffer implements a lock-free ring buffer for one writer and
+ * one reader, that is to be used to store a sample type T.
+ *
+ * RingBuffer is thread-safe provided only one thread writes and only
+ * one thread reads.
+ */
+
+template <typename T>
+class RingBuffer
+{
+public:
+    /**
+     * Create a ring buffer with room to write n samples.
+     *
+     * Note that the internal storage size will actually be n+1
+     * samples, as one element is unavailable for administrative
+     * reasons.  Since the ring buffer performs best if its size is a
+     * power of two, this means n should ideally be some power of two
+     * minus one.
+     */
+    RingBuffer(int n);
+
+    virtual ~RingBuffer();
+
+    /**
+     * Return the total capacity of the ring buffer in samples.
+     * (This is the argument n passed to the constructor.)
+     */
+    int getSize() const;
+
+    /**
+     * Return a new ring buffer (allocated with "new" -- caller must
+     * delete when no longer needed) of the given size, containing the
+     * same data as this one.  If another thread reads from or writes
+     * to this buffer during the call, the results may be incomplete
+     * or inconsistent.  If this buffer's data will not fit in the new
+     * size, the contents are undefined.
+     */
+    RingBuffer<T> *resized(int newSize) const;
+
+    /**
+     * Lock the ring buffer into physical memory.  Returns true
+     * for success.
+     */
+    bool mlock();
+
+    /**
+     * Reset read and write pointers, thus emptying the buffer.
+     * Should be called from the write thread.
+     */
+    void reset();
+
+    /**
+     * Return the amount of data available for reading, in samples.
+     */
+    int getReadSpace() const;
+
+    /**
+     * Return the amount of space available for writing, in samples.
+     */
+    int getWriteSpace() const;
+
+    /**
+     * Read n samples from the buffer.  If fewer than n are available,
+     * the remainder will be zeroed out.  Returns the number of
+     * samples actually read.
+     *
+     * This is a template function, taking an argument S for the target
+     * sample type, which is permitted to differ from T if the two
+     * types are compatible for arithmetic operations.
+     */
+    template <typename S>
+    int read(S *const R__ destination, int n);
+
+    /**
+     * Read n samples from the buffer, adding them to the destination.
+     * If fewer than n are available, the remainder will be left
+     * alone.  Returns the number of samples actually read.
+     *
+     * This is a template function, taking an argument S for the target
+     * sample type, which is permitted to differ from T if the two
+     * types are compatible for arithmetic operations.
+     */
+    template <typename S>
+    int readAdding(S *const R__ destination, int n);
+
+    /**
+     * Read one sample from the buffer.  If no sample is available,
+     * this will silently return zero.  Calling this repeatedly is
+     * obviously slower than calling read once, but it may be good
+     * enough if you don't want to allocate a buffer to read into.
+     */
+    T readOne();
+
+    /**
+     * Read n samples from the buffer, if available, without advancing
+     * the read pointer -- i.e. a subsequent read() or skip() will be
+     * necessary to empty the buffer.  If fewer than n are available,
+     * the remainder will be zeroed out.  Returns the number of
+     * samples actually read.
+     */
+    int peek(T *const R__ destination, int n) const;
+
+    /**
+     * Read one sample from the buffer, if available, without
+     * advancing the read pointer -- i.e. a subsequent read() or
+     * skip() will be necessary to empty the buffer.  Returns zero if
+     * no sample was available.
+     */
+    T peekOne() const;
+
+    /**
+     * Pretend to read n samples from the buffer, without actually
+     * returning them (i.e. discard the next n samples).  Returns the
+     * number of samples actually available for discarding.
+     */
+    int skip(int n);
+
+    /**
+     * Write n samples to the buffer.  If insufficient space is
+     * available, not all samples may actually be written.  Returns
+     * the number of samples actually written.
+     *
+     * This is a template function, taking an argument S for the source
+     * sample type, which is permitted to differ from T if the two
+     * types are compatible for assignment.
+     */
+    template <typename S>
+    int write(const S *const R__ source, int n);
+
+    /**
+     * Write n zero-value samples to the buffer.  If insufficient
+     * space is available, not all zeros may actually be written.
+     * Returns the number of zeroes actually written.
+     */
+    int zero(int n);
+
+protected:
+    T *const R__      m_buffer;
+    std::atomic<int>  m_writer;
+    std::atomic<int>  m_reader;
+    const int         m_size;
+    bool              m_mlocked;
+
+    int readSpaceFor(int w, int r) const {
+        int space;
+        if (w > r) space = w - r;
+        else if (w < r) space = (w + m_size) - r;
+        else space = 0;
+        return space;
+    }
+
+    int writeSpaceFor(int w, int r) const {
+        int space = (r + m_size - w - 1);
+        if (space >= m_size) space -= m_size;
+        return space;
+    }
+
+private:
+    RingBuffer(const RingBuffer &); // not provided
+    RingBuffer &operator=(const RingBuffer &); // not provided
+};
+
+template <typename T>
+RingBuffer<T>::RingBuffer(int n) :
+    m_buffer(allocate<T>(n + 1)),
+    m_writer(0),
+    m_size(n + 1),
+    m_mlocked(false)
+{
+#ifdef DEBUG_RINGBUFFER
+    std::cerr << "RingBuffer<T>[" << this << "]::RingBuffer(" << n << ")" << std::endl;
+#endif
+
+    m_reader = 0;
+}
+
+template <typename T>
+RingBuffer<T>::~RingBuffer()
+{
+#ifdef DEBUG_RINGBUFFER
+    std::cerr << "RingBuffer<T>[" << this << "]::~RingBuffer" << std::endl;
+#endif
+
+    if (m_mlocked) {
+	MUNLOCK((void *)m_buffer, m_size * sizeof(T));
+    }
+
+    deallocate(m_buffer);
+}
+
+template <typename T>
+int
+RingBuffer<T>::getSize() const
+{
+#ifdef DEBUG_RINGBUFFER
+    std::cerr << "RingBuffer<T>[" << this << "]::getSize(): " << m_size-1 << std::endl;
+#endif
+
+    return m_size - 1;
+}
+
+template <typename T>
+RingBuffer<T> *
+RingBuffer<T>::resized(int newSize) const
+{
+    RingBuffer<T> *newBuffer = new RingBuffer<T>(newSize);
+    
+    MBARRIER();
+    int w = m_writer;
+    int r = m_reader;
+
+    while (r != w) {
+        T value = m_buffer[r];
+        newBuffer->write(&value, 1);
+        if (++r == m_size) r = 0;
+    }
+
+    return newBuffer;
+}
+
+template <typename T>
+bool
+RingBuffer<T>::mlock()
+{
+    if (MLOCK((void *)m_buffer, m_size * sizeof(T))) return false;
+    m_mlocked = true;
+    return true;
+}
+
+template <typename T>
+void
+RingBuffer<T>::reset()
+{
+#ifdef DEBUG_RINGBUFFER
+    std::cerr << "RingBuffer<T>[" << this << "]::reset" << std::endl;
+#endif
+
+    int r = m_reader;
+    m_writer = r;
+}
+
+template <typename T>
+int
+RingBuffer<T>::getReadSpace() const
+{
+    return readSpaceFor(m_writer, m_reader);
+}
+
+template <typename T>
+int
+RingBuffer<T>::getWriteSpace() const
+{
+    return writeSpaceFor(m_writer, m_reader);
+}
+
+template <typename T>
+template <typename S>
+int
+RingBuffer<T>::read(S *const R__ destination, int n)
+{
+    int w = m_writer;
+    int r = m_reader;
+
+    int available = readSpaceFor(w, r);
+    if (n > available) {
+	std::cerr << "WARNING: RingBuffer::read: " << n << " requested, only "
+                  << available << " available" << std::endl;
+	n = available;
+    }
+    if (n == 0) return n;
+
+    int here = m_size - r;
+    T *const R__ bufbase = m_buffer + r;
+
+    if (here >= n) {
+        v_convert(destination, bufbase, n);
+    } else {
+        v_convert(destination, bufbase, here);
+        v_convert(destination + here, m_buffer, n - here);
+    }
+
+    r += n;
+    while (r >= m_size) r -= m_size;
+
+    m_reader = r;
+
+    return n;
+}
+
+template <typename T>
+template <typename S>
+int
+RingBuffer<T>::readAdding(S *const R__ destination, int n)
+{
+    int w = m_writer;
+    int r = m_reader;
+
+    int available = readSpaceFor(w, r);
+    if (n > available) {
+	std::cerr << "WARNING: RingBuffer::read: " << n << " requested, only "
+                  << available << " available" << std::endl;
+	n = available;
+    }
+    if (n == 0) return n;
+
+    int here = m_size - r;
+    T *const R__ bufbase = m_buffer + r;
+
+    if (here >= n) {
+        v_add(destination, bufbase, n);
+    } else {
+        v_add(destination, bufbase, here);
+        v_add(destination + here, m_buffer, n - here);
+    }
+
+    r += n;
+    while (r >= m_size) r -= m_size;
+
+    m_reader = r;
+
+    return n;
+}
+
+template <typename T>
+T
+RingBuffer<T>::readOne()
+{
+    int w = m_writer;
+    int r = m_reader;
+
+    if (w == r) {
+	std::cerr << "WARNING: RingBuffer::readOne: no sample available"
+		  << std::endl;
+	return T();
+    }
+
+    T value = m_buffer[r];
+    if (++r == m_size) r = 0;
+
+    m_reader = r;
+
+    return value;
+}
+
+template <typename T>
+int
+RingBuffer<T>::peek(T *const R__ destination, int n) const
+{
+    int w = m_writer;
+    int r = m_reader;
+
+    int available = readSpaceFor(w, r);
+    if (n > available) {
+	std::cerr << "WARNING: RingBuffer::peek: " << n << " requested, only "
+                  << available << " available" << std::endl;
+	memset(destination + available, 0, (n - available) * sizeof(T));
+	n = available;
+    }
+    if (n == 0) return n;
+
+    int here = m_size - r;
+    const T *const R__ bufbase = m_buffer + r;
+
+    if (here >= n) {
+        v_copy(destination, bufbase, n);
+    } else {
+        v_copy(destination, bufbase, here);
+        v_copy(destination + here, m_buffer, n - here);
+    }
+
+    return n;
+}
+
+template <typename T>
+T
+RingBuffer<T>::peekOne() const
+{
+    int w = m_writer;
+    int r = m_reader;
+
+    if (w == r) {
+	std::cerr << "WARNING: RingBuffer::peekOne: no sample available"
+		  << std::endl;
+	return 0;
+    }
+
+    T value = m_buffer[r];
+    return value;
+}
+
+template <typename T>
+int
+RingBuffer<T>::skip(int n)
+{
+    int w = m_writer;
+    int r = m_reader;
+
+    int available = readSpaceFor(w, r);
+    if (n > available) {
+	std::cerr << "WARNING: RingBuffer::skip: " << n << " requested, only "
+                  << available << " available" << std::endl;
+	n = available;
+    }
+    if (n == 0) return n;
+
+    r += n;
+    while (r >= m_size) r -= m_size;
+
+    m_reader = r;
+
+    return n;
+}
+
+template <typename T>
+template <typename S>
+int
+RingBuffer<T>::write(const S *const R__ source, int n)
+{
+    int w = m_writer;
+    int r = m_reader;
+
+    int available = writeSpaceFor(w, r);
+    if (n > available) {
+	std::cerr << "WARNING: RingBuffer::write: " << n
+                  << " requested, only room for " << available << std::endl;
+	n = available;
+    }
+    if (n == 0) return n;
+
+    int here = m_size - w;
+    T *const R__ bufbase = m_buffer + w;
+
+    if (here >= n) {
+        v_convert<S, T>(bufbase, source, n);
+    } else {
+        v_convert<S, T>(bufbase, source, here);
+        v_convert<S, T>(m_buffer, source + here, n - here);
+    }
+
+    w += n;
+    while (w >= m_size) w -= m_size;
+
+    MBARRIER();
+    m_writer = w;
+
+    return n;
+}
+
+template <typename T>
+int
+RingBuffer<T>::zero(int n)
+{
+    int w = m_writer;
+    int r = m_reader;
+
+    int available = writeSpaceFor(w, r);
+    if (n > available) {
+	std::cerr << "WARNING: RingBuffer::zero: " << n
+                  << " requested, only room for " << available << std::endl;
+	n = available;
+    }
+    if (n == 0) return n;
+
+    int here = m_size - w;
+    T *const R__ bufbase = m_buffer + w;
+
+    if (here >= n) {
+        v_zero(bufbase, n);
+    } else {
+        v_zero(bufbase, here);
+        v_zero(m_buffer, n - here);
+    }
+
+    w += n;
+    while (w >= m_size) w -= m_size;
+
+    MBARRIER();
+    m_writer = w;
+
+    return n;
+}
+
+}
+
+#endif // _RINGBUFFER_H

src/common/SampleFilter.h

@@ -0,0 +1,59 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_SAMPLE_FILTER_H
+#define RUBBERBAND_SAMPLE_FILTER_H
+
+#include <cassert>
+
+namespace RubberBand
+{
+
+template <typename T>
+class SampleFilter
+{
+public:
+    SampleFilter(int size) : m_size(size) {
+	assert(m_size > 0);
+    }
+
+    virtual ~SampleFilter() { }
+
+    int getSize() const { return m_size; }
+
+    virtual void push(T) = 0;
+    virtual T get() const = 0;
+    virtual void reset() = 0;
+
+protected:
+    const int m_size;
+
+private:
+    SampleFilter(const SampleFilter &);
+    SampleFilter &operator=(const SampleFilter &);
+};
+
+}
+
+#endif
+

src/common/Scavenger.h

@@ -0,0 +1,250 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_SCAVENGER_H
+#define RUBBERBAND_SCAVENGER_H
+
+#include <vector>
+#include <list>
+#include <utility>
+#include <iostream>
+
+#ifndef _MSC_VER
+#include <sys/time.h>
+#endif
+
+#include "Thread.h"
+#include "sysutils.h"
+#include "Allocators.h"
+
+//#define DEBUG_SCAVENGER 1
+
+namespace RubberBand {
+
+/**
+ * A very simple class that facilitates running things like plugins
+ * without locking, by collecting unwanted objects and deleting them
+ * after a delay so as to be sure nobody's in the middle of using
+ * them.  Requires scavenge() to be called regularly from a non-RT
+ * thread.
+ *
+ * This is currently not at all suitable for large numbers of objects
+ * -- it's just a quick hack for use with things like plugins.
+ */
+
+//!!! should review this, it's not really thread safe owing to lack of
+//!!! atomic updates
+
+template <typename T>
+class Scavenger
+{
+public:
+    Scavenger(int sec = 2, int defaultObjectListSize = 200);
+    ~Scavenger();
+
+    /**
+     * Call from an RT thread etc., to pass ownership of t to us.
+     * Only one thread should be calling this on any given scavenger.
+     */
+    void claim(T *t);
+
+    /**
+     * Call from a non-RT thread.
+     * Only one thread should be calling this on any given scavenger.
+     */
+    void scavenge(bool clearNow = false);
+
+protected:
+    typedef std::pair<T *, int> ObjectTimePair;
+    typedef std::vector<ObjectTimePair> ObjectTimeList;
+    ObjectTimeList m_objects;
+    int m_sec;
+
+    typedef std::list<T *> ObjectList;
+    ObjectList m_excess;
+    int m_lastExcess;
+    Mutex m_excessMutex;
+    void pushExcess(T *);
+    void clearExcess(int);
+
+    unsigned int m_claimed;
+    unsigned int m_scavenged;
+    unsigned int m_asExcess;
+};
+
+
+/**
+ * A wrapper to permit arrays allocated with new[] to be scavenged.
+ */
+
+template <typename T>
+class ScavengerArrayWrapper
+{
+public:
+    ScavengerArrayWrapper(T *array) : m_array(array) { }
+    ~ScavengerArrayWrapper() { delete[] m_array; }
+
+private:
+    T *m_array;
+};
+
+
+/**
+ * A wrapper to permit arrays allocated with the Allocators functions
+ * to be scavenged.
+ */
+
+template <typename T>
+class ScavengerAllocArrayWrapper
+{
+public:
+    ScavengerAllocArrayWrapper(T *array) : m_array(array) { }
+    ~ScavengerAllocArrayWrapper() { deallocate<T>(m_array); }
+
+private:
+    T *m_array;
+};
+
+
+template <typename T>
+Scavenger<T>::Scavenger(int sec, int defaultObjectListSize) :
+    m_objects(ObjectTimeList(defaultObjectListSize)),
+    m_sec(sec),
+    m_lastExcess(0),
+    m_claimed(0),
+    m_scavenged(0),
+    m_asExcess(0)
+{
+}
+
+template <typename T>
+Scavenger<T>::~Scavenger()
+{
+    if (m_scavenged < m_claimed) {
+	for (size_t i = 0; i < m_objects.size(); ++i) {
+	    ObjectTimePair &pair = m_objects[i];
+	    if (pair.first != 0) {
+		T *ot = pair.first;
+		pair.first = 0;
+		delete ot;
+		++m_scavenged;
+	    }
+	}
+    }
+
+    clearExcess(0);
+}
+
+template <typename T>
+void
+Scavenger<T>::claim(T *t)
+{
+//    std::cerr << "Scavenger::claim(" << t << ")" << std::endl;
+
+    struct timeval tv;
+    (void)gettimeofday(&tv, 0);
+    int sec = tv.tv_sec;
+
+    for (size_t i = 0; i < m_objects.size(); ++i) {
+	ObjectTimePair &pair = m_objects[i];
+	if (pair.first == 0) {
+	    pair.second = sec;
+	    pair.first = t;
+	    ++m_claimed;
+	    return;
+	}
+    }
+
+#ifdef DEBUG_SCAVENGER
+    std::cerr << "WARNING: Scavenger::claim(" << t << "): run out of slots (at "
+              << m_objects.size() << "), using non-RT-safe method" << std::endl;
+#endif
+    pushExcess(t);
+}
+
+template <typename T>
+void
+Scavenger<T>::scavenge(bool clearNow)
+{
+#ifdef DEBUG_SCAVENGER
+    std::cerr << "Scavenger::scavenge: claimed " << m_claimed << ", scavenged " << m_scavenged << ", cleared as excess " << m_asExcess << std::endl;
+#endif
+
+    if (m_scavenged >= m_claimed) return;
+    
+    struct timeval tv;
+    (void)gettimeofday(&tv, 0);
+    int sec = tv.tv_sec;
+    bool anything = false;
+
+    for (size_t i = 0; i < m_objects.size(); ++i) {
+	ObjectTimePair &pair = m_objects[i];
+        if (!pair.first) continue;
+	if (clearNow || pair.second + m_sec < sec) {
+	    T *ot = pair.first;
+	    pair.first = 0;
+	    delete ot;
+	    ++m_scavenged;
+            anything = true;
+	}
+    }
+
+    if (clearNow || anything || (sec > m_lastExcess + m_sec)) {
+        clearExcess(sec);
+    }
+}
+
+template <typename T>
+void
+Scavenger<T>::pushExcess(T *t)
+{
+    m_excessMutex.lock();
+    m_excess.push_back(t);
+    struct timeval tv;
+    (void)gettimeofday(&tv, 0);
+    m_lastExcess = tv.tv_sec;
+    m_excessMutex.unlock();
+}
+
+template <typename T>
+void
+Scavenger<T>::clearExcess(int sec)
+{
+#ifdef DEBUG_SCAVENGER
+    std::cerr << "Scavenger::clearExcess: Excess now " << m_excess.size() << std::endl;
+#endif
+
+    m_excessMutex.lock();
+    for (typename ObjectList::iterator i = m_excess.begin();
+	 i != m_excess.end(); ++i) {
+	delete *i;
+        ++m_asExcess;
+    }
+    m_excess.clear();
+    m_lastExcess = sec;
+    m_excessMutex.unlock();
+}
+
+}
+
+#endif

src/common/Thread.cpp

@@ -0,0 +1,671 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef NO_THREADING
+
+#include "Thread.h"
+
+#include <iostream>
+#include <cstdlib>
+
+#ifdef USE_PTHREADS
+#include <sys/time.h>
+#include <time.h>
+#endif
+
+using std::cerr;
+using std::endl;
+using std::string;
+
+namespace RubberBand
+{
+
+#ifdef _WIN32
+
+Thread::Thread() :
+    m_id(0),
+    m_extant(false)
+{
+#ifdef DEBUG_THREAD
+    cerr << "THREAD DEBUG: Created thread object " << this << endl;
+#endif
+}
+
+Thread::~Thread()
+{
+#ifdef DEBUG_THREAD
+    cerr << "THREAD DEBUG: Destroying thread object " << this << ", id " << m_id << endl;
+#endif
+    if (m_extant) {
+        WaitForSingleObject(m_id, INFINITE);
+    }
+#ifdef DEBUG_THREAD
+    cerr << "THREAD DEBUG: Destroyed thread object " << this << endl;
+#endif
+}
+
+void
+Thread::start()
+{
+    m_id = CreateThread(NULL, 0, staticRun, this, 0, 0);
+    if (!m_id) {
+        cerr << "ERROR: thread creation failed" << endl;
+        exit(1);
+    } else {
+#ifdef DEBUG_THREAD
+        cerr << "THREAD DEBUG: Created thread " << m_id << " for thread object " << this << endl;
+#endif
+        m_extant = true;
+    }
+}    
+
+void 
+Thread::wait()
+{
+    if (m_extant) {
+#ifdef DEBUG_THREAD
+        cerr << "THREAD DEBUG: Waiting on thread " << m_id << " for thread object " << this << endl;
+#endif
+        WaitForSingleObject(m_id, INFINITE);
+#ifdef DEBUG_THREAD
+        cerr << "THREAD DEBUG: Waited on thread " << m_id << " for thread object " << this << endl;
+#endif
+        m_extant = false;
+    }
+}
+
+Thread::Id
+Thread::id()
+{
+    return m_id;
+}
+
+bool
+Thread::threadingAvailable()
+{
+    return true;
+}
+
+DWORD
+Thread::staticRun(LPVOID arg)
+{
+    Thread *thread = static_cast<Thread *>(arg);
+#ifdef DEBUG_THREAD
+    cerr << "THREAD DEBUG: " << (void *)GetCurrentThreadId() << ": Running thread " << thread->m_id << " for thread object " << thread << endl;
+#endif
+    thread->run();
+    return 0;
+}
+
+Mutex::Mutex()
+#ifndef NO_THREAD_CHECKS
+    :
+    m_lockedBy(-1)
+#endif
+{
+    m_mutex = CreateMutex(NULL, FALSE, NULL);
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)GetCurrentThreadId() << ": Initialised mutex " << &m_mutex << endl;
+#endif
+}
+
+Mutex::~Mutex()
+{
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)GetCurrentThreadId() << ": Destroying mutex " << &m_mutex << endl;
+#endif
+    CloseHandle(m_mutex);
+}
+
+void
+Mutex::lock()
+{
+#ifndef NO_THREAD_CHECKS
+    DWORD tid = GetCurrentThreadId();
+    if (m_lockedBy == tid) {
+        cerr << "ERROR: Deadlock on mutex " << &m_mutex << endl;
+    }
+#endif
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)tid << ": Want to lock mutex " << &m_mutex << endl;
+#endif
+    WaitForSingleObject(m_mutex, INFINITE);
+#ifndef NO_THREAD_CHECKS
+    m_lockedBy = tid;
+#endif
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)tid << ": Locked mutex " << &m_mutex << endl;
+#endif
+}
+
+void
+Mutex::unlock()
+{
+#ifndef NO_THREAD_CHECKS
+    DWORD tid = GetCurrentThreadId();
+    if (m_lockedBy != tid) {
+        cerr << "ERROR: Mutex " << &m_mutex << " not owned by unlocking thread" << endl;
+        return;
+    }
+#endif
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)tid << ": Unlocking mutex " << &m_mutex << endl;
+#endif
+#ifndef NO_THREAD_CHECKS
+    m_lockedBy = -1;
+#endif
+    ReleaseMutex(m_mutex);
+}
+
+bool
+Mutex::trylock()
+{
+#ifndef NO_THREAD_CHECKS
+    DWORD tid = GetCurrentThreadId();
+#endif
+    DWORD result = WaitForSingleObject(m_mutex, 0);
+    if (result == WAIT_TIMEOUT || result == WAIT_FAILED) {
+#ifdef DEBUG_MUTEX
+        cerr << "MUTEX DEBUG: " << (void *)tid << ": Mutex " << &m_mutex << " unavailable" << endl;
+#endif
+        return false;
+    } else {
+#ifndef NO_THREAD_CHECKS
+        m_lockedBy = tid;
+#endif
+#ifdef DEBUG_MUTEX
+        cerr << "MUTEX DEBUG: " << (void *)tid << ": Locked mutex " << &m_mutex << " (from trylock)" << endl;
+#endif
+        return true;
+    }
+}
+
+Condition::Condition(string
+#ifdef DEBUG_CONDITION
+                     name
+#endif
+    ) :
+    m_locked(false)
+#ifdef DEBUG_CONDITION
+    , m_name(name)
+#endif
+{
+    m_mutex = CreateMutex(NULL, FALSE, NULL);
+    m_condition = CreateEvent(NULL, FALSE, FALSE, NULL);
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Initialised condition " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+}
+
+Condition::~Condition()
+{
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Destroying condition " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    if (m_locked) ReleaseMutex(m_mutex);
+    CloseHandle(m_condition);
+    CloseHandle(m_mutex);
+}
+
+void
+Condition::lock()
+{
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Want to lock " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    WaitForSingleObject(m_mutex, INFINITE);
+    m_locked = true;
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Locked " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+}
+
+void
+Condition::unlock()
+{
+    if (!m_locked) {
+#ifdef DEBUG_CONDITION
+        cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Not locked " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+        return;
+    }
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Unlocking " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    m_locked = false;
+    ReleaseMutex(m_mutex);
+}
+
+void 
+Condition::wait(int us)
+{
+    if (us == 0) {
+
+#ifdef DEBUG_CONDITION
+        cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Waiting on " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+        SignalObjectAndWait(m_mutex, m_condition, INFINITE, FALSE);
+        WaitForSingleObject(m_mutex, INFINITE);
+
+    } else {
+
+        DWORD ms = us / 1000;
+        if (us > 0 && ms == 0) ms = 1;
+    
+#ifdef DEBUG_CONDITION
+        cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Timed waiting on " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+        SignalObjectAndWait(m_mutex, m_condition, ms, FALSE);
+        WaitForSingleObject(m_mutex, INFINITE);
+    }
+
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Wait done on " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    m_locked = true;
+}
+
+void
+Condition::signal()
+{
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)GetCurrentThreadId() << ": Signalling " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    SetEvent(m_condition);
+}
+
+#else /* !_WIN32 */
+
+#ifdef USE_PTHREADS
+
+Thread::Thread() :
+    m_id(0),
+    m_extant(false)
+{
+#ifdef DEBUG_THREAD
+    cerr << "THREAD DEBUG: Created thread object " << this << endl;
+#endif
+}
+
+Thread::~Thread()
+{
+#ifdef DEBUG_THREAD
+    cerr << "THREAD DEBUG: Destroying thread object " << this << ", id " << m_id << endl;
+#endif
+    if (m_extant) {
+        pthread_join(m_id, 0);
+    }
+#ifdef DEBUG_THREAD
+    cerr << "THREAD DEBUG: Destroyed thread object " << this << endl;
+#endif
+}
+
+void
+Thread::start()
+{
+    if (pthread_create(&m_id, 0, staticRun, this)) {
+        cerr << "ERROR: thread creation failed" << endl;
+        exit(1);
+    } else {
+#ifdef DEBUG_THREAD
+        cerr << "THREAD DEBUG: Created thread " << m_id << " for thread object " << this << endl;
+#endif
+        m_extant = true;
+    }
+}    
+
+void 
+Thread::wait()
+{
+    if (m_extant) {
+#ifdef DEBUG_THREAD
+        cerr << "THREAD DEBUG: Waiting on thread " << m_id << " for thread object " << this << endl;
+#endif
+        pthread_join(m_id, 0);
+#ifdef DEBUG_THREAD
+        cerr << "THREAD DEBUG: Waited on thread " << m_id << " for thread object " << this << endl;
+#endif
+        m_extant = false;
+    }
+}
+
+Thread::Id
+Thread::id()
+{
+    return m_id;
+}
+
+bool
+Thread::threadingAvailable()
+{
+    return true;
+}
+
+void *
+Thread::staticRun(void *arg)
+{
+    Thread *thread = static_cast<Thread *>(arg);
+#ifdef DEBUG_THREAD
+    cerr << "THREAD DEBUG: " << (void *)pthread_self() << ": Running thread " << thread->m_id << " for thread object " << thread << endl;
+#endif
+    thread->run();
+    return 0;
+}
+
+Mutex::Mutex()
+#ifndef NO_THREAD_CHECKS
+    :
+    m_lockedBy(0),
+    m_locked(false)
+#endif
+{
+    pthread_mutex_init(&m_mutex, 0);
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)pthread_self() << ": Initialised mutex " << &m_mutex << endl;
+#endif
+}
+
+Mutex::~Mutex()
+{
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)pthread_self() << ": Destroying mutex " << &m_mutex << endl;
+#endif
+    pthread_mutex_destroy(&m_mutex);
+}
+
+void
+Mutex::lock()
+{
+#ifndef NO_THREAD_CHECKS
+    pthread_t tid = pthread_self();
+    if (m_locked && m_lockedBy == tid) {
+        cerr << "ERROR: Deadlock on mutex " << &m_mutex << endl;
+    }
+#endif
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)tid << ": Want to lock mutex " << &m_mutex << endl;
+#endif
+    pthread_mutex_lock(&m_mutex);
+#ifndef NO_THREAD_CHECKS
+    m_lockedBy = tid;
+    m_locked = true;
+#endif
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)tid << ": Locked mutex " << &m_mutex << endl;
+#endif
+}
+
+void
+Mutex::unlock()
+{
+#ifndef NO_THREAD_CHECKS
+    pthread_t tid = pthread_self();
+    if (!m_locked) {
+        cerr << "ERROR: Mutex " << &m_mutex << " not locked in unlock" << endl;
+        return;
+    } else if (m_lockedBy != tid) {
+        cerr << "ERROR: Mutex " << &m_mutex << " not owned by unlocking thread" << endl;
+        return;
+    }
+#endif
+#ifdef DEBUG_MUTEX
+    cerr << "MUTEX DEBUG: " << (void *)tid << ": Unlocking mutex " << &m_mutex << endl;
+#endif
+#ifndef NO_THREAD_CHECKS
+    m_locked = false;
+#endif
+    pthread_mutex_unlock(&m_mutex);
+}
+
+bool
+Mutex::trylock()
+{
+#ifndef NO_THREAD_CHECKS
+    pthread_t tid = pthread_self();
+#endif
+    if (pthread_mutex_trylock(&m_mutex)) {
+#ifdef DEBUG_MUTEX
+        cerr << "MUTEX DEBUG: " << (void *)tid << ": Mutex " << &m_mutex << " unavailable" << endl;
+#endif
+        return false;
+    } else {
+#ifndef NO_THREAD_CHECKS
+        m_lockedBy = tid;
+        m_locked = true;
+#endif
+#ifdef DEBUG_MUTEX
+        cerr << "MUTEX DEBUG: " << (void *)tid << ": Locked mutex " << &m_mutex << " (from trylock)" << endl;
+#endif
+        return true;
+    }
+}
+
+Condition::Condition(string
+#ifdef DEBUG_CONDITION
+                     name
+#endif
+    ) :
+    m_locked(false)
+#ifdef DEBUG_CONDITION
+    , m_name(name)
+#endif
+{
+    pthread_mutex_init(&m_mutex, 0);
+    pthread_cond_init(&m_condition, 0);
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Initialised condition " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+}
+
+Condition::~Condition()
+{
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Destroying condition " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    if (m_locked) pthread_mutex_unlock(&m_mutex);
+    pthread_cond_destroy(&m_condition);
+    pthread_mutex_destroy(&m_mutex);
+}
+
+void
+Condition::lock()
+{
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Want to lock " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    pthread_mutex_lock(&m_mutex);
+    m_locked = true;
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Locked " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+}
+
+void
+Condition::unlock()
+{
+    if (!m_locked) {
+#ifdef DEBUG_CONDITION
+        cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Not locked " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+        return;
+    }
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Unlocking " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    m_locked = false;
+    pthread_mutex_unlock(&m_mutex);
+}
+
+void 
+Condition::wait(int us)
+{
+    if (us == 0) {
+
+#ifdef DEBUG_CONDITION
+        cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Waiting on " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+        pthread_cond_wait(&m_condition, &m_mutex);
+
+    } else {
+
+        struct timeval now;
+        gettimeofday(&now, 0);
+
+        now.tv_usec += us;
+        while (now.tv_usec > 1000000) {
+            now.tv_usec -= 1000000;
+            ++now.tv_sec;
+        }
+
+        struct timespec timeout;
+        timeout.tv_sec = now.tv_sec;
+        timeout.tv_nsec = now.tv_usec * 1000;
+    
+#ifdef DEBUG_CONDITION
+        cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Timed waiting on " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+        pthread_cond_timedwait(&m_condition, &m_mutex, &timeout);
+    }
+
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Wait done on " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    m_locked = true;
+}
+
+void
+Condition::signal()
+{
+#ifdef DEBUG_CONDITION
+    cerr << "CONDITION DEBUG: " << (void *)pthread_self() << ": Signalling " << &m_condition << " \"" << m_name << "\"" << endl;
+#endif
+    pthread_cond_signal(&m_condition);
+}
+
+#else /* !USE_PTHREADS */
+
+Thread::Thread()
+{
+}
+
+Thread::~Thread()
+{
+}
+
+void
+Thread::start()
+{
+    abort();
+}    
+
+void 
+Thread::wait()
+{
+    abort();
+}
+
+Thread::Id
+Thread::id()
+{
+    abort();
+}
+
+bool
+Thread::threadingAvailable()
+{
+    return false;
+}
+
+Mutex::Mutex()
+{
+}
+
+Mutex::~Mutex()
+{
+}
+
+void
+Mutex::lock()
+{
+    abort();
+}
+
+void
+Mutex::unlock()
+{
+    abort();
+}
+
+bool
+Mutex::trylock()
+{
+    abort();
+}
+
+Condition::Condition(const char *)
+{
+}
+
+Condition::~Condition()
+{
+}
+
+void
+Condition::lock()
+{
+    abort();
+}
+
+void 
+Condition::wait(int us)
+{
+    abort();
+}
+
+void
+Condition::signal()
+{
+    abort();
+}
+
+#endif /* !USE_PTHREADS */
+#endif /* !_WIN32 */
+
+MutexLocker::MutexLocker(Mutex *mutex) :
+    m_mutex(mutex)
+{
+    if (m_mutex) {
+        m_mutex->lock();
+    }
+}
+
+MutexLocker::~MutexLocker()
+{
+    if (m_mutex) {
+        m_mutex->unlock();
+    }
+}
+
+}
+
+#endif

src/common/Thread.h

@@ -0,0 +1,232 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_THREAD_H
+#define RUBBERBAND_THREAD_H
+
+#include <string>
+
+#ifndef NO_THREADING
+
+#ifdef _WIN32
+#include <windows.h>
+#else /* !_WIN32 */
+#ifdef USE_PTHREADS
+#include <pthread.h>
+#else /* !USE_PTHREADS */
+#error No thread implementation selected
+#endif /* !USE_PTHREADS */
+#endif /* !_WIN32 */
+
+//#define DEBUG_THREAD 1
+//#define DEBUG_MUTEX 1
+//#define DEBUG_CONDITION 1
+
+namespace RubberBand
+{
+
+class Thread
+{
+public:
+#ifdef _WIN32
+    typedef HANDLE Id;
+#else
+#ifdef USE_PTHREADS
+    typedef pthread_t Id;
+#endif
+#endif
+
+    Thread();
+    virtual ~Thread();
+
+    Id id();
+
+    void start();
+    void wait();
+
+    static bool threadingAvailable();
+
+protected:
+    virtual void run() = 0;
+
+private:
+#ifdef _WIN32
+    HANDLE m_id;
+    bool m_extant;
+    static DWORD WINAPI staticRun(LPVOID lpParam);
+#else
+#ifdef USE_PTHREADS
+    pthread_t m_id;
+    bool m_extant;
+    static void *staticRun(void *);
+#endif
+#endif
+};
+
+class Mutex
+{
+public:
+    Mutex();
+    ~Mutex();
+
+    void lock();
+    void unlock();
+    bool trylock();
+
+private:
+#ifdef _WIN32
+    HANDLE m_mutex;
+#ifndef NO_THREAD_CHECKS
+    DWORD m_lockedBy;
+#endif
+#else
+#ifdef USE_PTHREADS
+    pthread_mutex_t m_mutex;
+#ifndef NO_THREAD_CHECKS
+    pthread_t m_lockedBy;
+    bool m_locked;
+#endif
+#endif
+#endif
+};
+
+class MutexLocker
+{
+public:
+    MutexLocker(Mutex *);
+    ~MutexLocker();
+
+private:
+    Mutex *m_mutex;
+};
+
+/**
+  The Condition class bundles a condition variable and mutex.
+
+  To wait on a condition, call lock(), test the termination condition
+  if desired, then wait().  The condition will be unlocked during the
+  wait and re-locked when wait() returns (which will happen when the
+  condition is signalled or the timer times out).
+
+  To signal a condition, call signal().  If the condition is signalled
+  between lock() and wait(), the signal may be missed by the waiting
+  thread.  To avoid this, the signalling thread should also lock the
+  condition before calling signal() and unlock it afterwards.
+*/
+
+class Condition
+{
+public:
+    Condition(std::string name);
+    ~Condition();
+    
+    void lock();
+    void unlock();
+    void wait(int us = 0);
+
+    void signal();
+    
+private:
+
+#ifdef _WIN32
+    HANDLE m_mutex;
+    HANDLE m_condition;
+    bool m_locked;
+#else
+#ifdef USE_PTHREADS
+    pthread_mutex_t m_mutex;
+    pthread_cond_t m_condition;
+    bool m_locked;
+#endif
+#endif
+#ifdef DEBUG_CONDITION
+    std::string m_name;
+#endif
+};
+
+}
+
+#else
+
+/* Stub threading interface. We do not have threading support in this code. */
+
+namespace RubberBand
+{
+
+class Thread
+{
+public:
+    typedef unsigned int Id;
+
+    Thread() { }
+    virtual ~Thread() { }
+
+    Id id() { return 0; }
+
+    void start() { } 
+    void wait() { }
+
+    static bool threadingAvailable() { return false; }
+
+protected:
+    virtual void run() = 0;
+
+private:
+};
+
+class Mutex
+{
+public:
+    Mutex() { }
+    ~Mutex() { }
+
+    void lock() { }
+    void unlock() { }
+    bool trylock() { return false; }
+};
+
+class MutexLocker
+{
+public:
+    MutexLocker(Mutex *) { }
+    ~MutexLocker() { }
+};
+
+class Condition
+{
+public:
+    Condition(std::string name) { }
+    ~Condition() { }
+    
+    void lock() { }
+    void unlock() { }
+    void wait(int us = 0) { }
+
+    void signal() { }
+};
+
+}
+
+#endif /* NO_THREADING */
+
+#endif

src/common/VectorOps.h

@@ -0,0 +1,869 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_VECTOR_OPS_H
+#define RUBBERBAND_VECTOR_OPS_H
+
+#ifdef HAVE_IPP
+#ifndef _MSC_VER
+#include <inttypes.h>
+#endif
+#include <ippversion.h>
+#include <ipps.h>
+#if (IPP_VERSION_MAJOR <= 7)
+// Deprecated in v8, removed in v9
+#include <ippac.h>
+#endif
+#endif
+
+#ifdef HAVE_VDSP
+#include <Accelerate/Accelerate.h>
+#include <TargetConditionals.h>
+#include <alloca.h>
+#endif
+
+#include <cstring>
+#include "sysutils.h"
+
+namespace RubberBand {
+
+// Note that all functions with a "target" vector have their arguments
+// in the same order as memcpy and friends, i.e. target vector first.
+// This is the reverse order from the IPP functions.
+
+// The ideal here is to write the basic loops in such a way as to be
+// auto-vectorizable by a sensible compiler (definitely gcc-4.3 on
+// Linux, ideally also gcc-4.0 on OS/X).
+
+template<typename T>
+inline void v_zero(T *const R__ ptr,
+                   const int count)
+{
+    const T value = T(0);
+    for (int i = 0; i < count; ++i) {
+        ptr[i] = value;
+    }
+}
+
+#if defined HAVE_IPP
+template<> 
+inline void v_zero(float *const R__ ptr, 
+                   const int count)
+{
+    ippsZero_32f(ptr, count);
+}
+template<> 
+inline void v_zero(double *const R__ ptr,
+                   const int count)
+{
+    ippsZero_64f(ptr, count);
+}
+#elif defined HAVE_VDSP
+template<> 
+inline void v_zero(float *const R__ ptr, 
+                   const int count)
+{
+    vDSP_vclr(ptr, 1, count);
+}
+template<> 
+inline void v_zero(double *const R__ ptr,
+                   const int count)
+{
+    vDSP_vclrD(ptr, 1, count);
+}
+#endif
+
+template<typename T>
+inline void v_zero_channels(T *const R__ *const R__ ptr,
+                            const int channels,
+                            const int count)
+{
+    for (int c = 0; c < channels; ++c) {
+        v_zero(ptr[c], count);
+    }
+}
+
+template<typename T>
+inline void v_set(T *const R__ ptr,
+                  const T value,
+                  const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        ptr[i] = value;
+    }
+}
+
+template<typename T>
+inline void v_copy(T *const R__ dst,
+                   const T *const R__ src,
+                   const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] = src[i];
+    }
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_copy(float *const R__ dst,
+                   const float *const R__ src,
+                   const int count)
+{
+    ippsCopy_32f(src, dst, count);
+}
+template<>
+inline void v_copy(double *const R__ dst,
+                   const double *const R__ src,
+                   const int count)
+{
+    ippsCopy_64f(src, dst, count);
+}
+#endif
+
+template<typename T>
+inline void v_copy_channels(T *const R__ *const R__ dst,
+                            const T *const R__ *const R__ src,
+                            const int channels,
+                            const int count)
+{
+    for (int c = 0; c < channels; ++c) {
+        v_copy(dst[c], src[c], count);
+    }
+}
+
+// src and dst alias by definition, so not restricted
+template<typename T>
+inline void v_move(T *const dst,
+                   const T *const src,
+                   const int count)
+{
+    memmove(dst, src, count * sizeof(T));
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_move(float *const dst,
+                   const float *const src,
+                   const int count)
+{
+    ippsMove_32f(src, dst, count);
+}
+template<>
+inline void v_move(double *const dst,
+                   const double *const src,
+                   const int count)
+{
+    ippsMove_64f(src, dst, count);
+}
+#endif
+
+template<typename T, typename U>
+inline void v_convert(U *const R__ dst,
+                      const T *const R__ src,
+                      const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] = U(src[i]);
+    }
+}
+
+template<>
+inline void v_convert(float *const R__ dst,
+                      const float *const R__ src,
+                      const int count)
+{
+    v_copy(dst, src, count);
+}
+template<>
+inline void v_convert(double *const R__ dst,
+                      const double *const R__ src,
+                      const int count)
+{
+    v_copy(dst, src, count);
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_convert(double *const R__ dst,
+                      const float *const R__ src,
+                      const int count)
+{
+    ippsConvert_32f64f(src, dst, count);
+}
+template<>
+inline void v_convert(float *const R__ dst,
+                      const double *const R__ src,
+                      const int count)
+{
+    ippsConvert_64f32f(src, dst, count);
+}
+#elif defined HAVE_VDSP
+template<>
+inline void v_convert(double *const R__ dst,
+                      const float *const R__ src,
+                      const int count)
+{
+    vDSP_vspdp((float *)src, 1, dst, 1, count);
+}
+template<>
+inline void v_convert(float *const R__ dst,
+                      const double *const R__ src,
+                      const int count)
+{
+    vDSP_vdpsp((double *)src, 1, dst, 1, count);
+}
+#endif
+
+template<typename T, typename U>
+inline void v_convert_channels(U *const R__ *const R__ dst,
+                               const T *const R__ *const R__ src,
+                               const int channels,
+                               const int count)
+{
+    for (int c = 0; c < channels; ++c) {
+        v_convert(dst[c], src[c], count);
+    }
+}
+
+template<typename T>
+inline void v_add(T *const R__ dst,
+                  const T *const R__ src,
+                  const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] += src[i];
+    }
+}
+
+template<typename T>
+inline void v_add(T *const R__ dst,
+                  const T value,
+                  const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] += value;
+    }
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_add(float *const R__ dst,
+                  const float *const R__ src,
+                  const int count)
+{
+    ippsAdd_32f_I(src, dst, count);
+}    
+inline void v_add(double *const R__ dst,
+                  const double *const R__ src,
+                  const int count)
+{
+    ippsAdd_64f_I(src, dst, count);
+}    
+#endif
+
+template<typename T>
+inline void v_add_channels(T *const R__ *const R__ dst,
+                           const T *const R__ *const R__ src,
+                           const int channels, const int count)
+{
+    for (int c = 0; c < channels; ++c) {
+        v_add(dst[c], src[c], count);
+    }
+}
+
+template<typename T, typename G>
+inline void v_add_with_gain(T *const R__ dst,
+                            const T *const R__ src,
+                            const G gain,
+                            const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] += src[i] * gain;
+    }
+}
+
+template<typename T, typename G>
+inline void v_add_channels_with_gain(T *const R__ *const R__ dst,
+                                     const T *const R__ *const R__ src,
+                                     const G gain,
+                                     const int channels,
+                                     const int count)
+{
+    for (int c = 0; c < channels; ++c) {
+        v_add_with_gain(dst[c], src[c], gain, count);
+    }
+}
+
+template<typename T>
+inline void v_subtract(T *const R__ dst,
+                       const T *const R__ src,
+                       const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] -= src[i];
+    }
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_subtract(float *const R__ dst,
+                       const float *const R__ src,
+                       const int count)
+{
+    ippsSub_32f_I(src, dst, count);
+}    
+inline void v_subtract(double *const R__ dst,
+                       const double *const R__ src,
+                       const int count)
+{
+    ippsSub_64f_I(src, dst, count);
+}    
+#endif
+
+template<typename T, typename G>
+inline void v_scale(T *const R__ dst,
+                    const G gain,
+                    const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] *= gain;
+    }
+}
+
+#if defined HAVE_IPP 
+template<>
+inline void v_scale(float *const R__ dst,
+                    const float gain,
+                    const int count)
+{
+    ippsMulC_32f_I(gain, dst, count);
+}
+template<>
+inline void v_scale(double *const R__ dst,
+                    const double gain,
+                    const int count)
+{
+    ippsMulC_64f_I(gain, dst, count);
+}
+#endif
+
+template<typename T, typename S>
+inline void v_multiply(T *const R__ srcdst,
+                       const S *const R__ src,
+                       const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        srcdst[i] *= src[i];
+    }
+}
+
+#if defined HAVE_IPP 
+template<>
+inline void v_multiply(float *const R__ srcdst,
+                       const float *const R__ src,
+                       const int count)
+{
+    ippsMul_32f_I(src, srcdst, count);
+}
+template<>
+inline void v_multiply(double *const R__ srcdst,
+                       const double *const R__ src,
+                       const int count)
+{
+    ippsMul_64f_I(src, srcdst, count);
+}
+#endif // HAVE_IPP
+
+template<typename T>
+inline void v_multiply(T *const R__ dst,
+                       const T *const R__ src1,
+                       const T *const R__ src2,
+                       const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] = src1[i] * src2[i];
+    }
+}
+
+template<typename T>
+inline void v_divide(T *const R__ dst,
+                     const T *const R__ src,
+                     const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] /= src[i];
+    }
+}
+
+#if defined HAVE_IPP 
+template<>
+inline void v_divide(float *const R__ dst,
+                     const float *const R__ src,
+                     const int count)
+{
+    ippsDiv_32f_I(src, dst, count);
+}
+template<>
+inline void v_divide(double *const R__ dst,
+                     const double *const R__ src,
+                     const int count)
+{
+    ippsDiv_64f_I(src, dst, count);
+}
+#endif
+
+#if defined HAVE_IPP 
+template<>
+inline void v_multiply(float *const R__ dst,
+                       const float *const R__ src1,
+                       const float *const R__ src2,
+                       const int count)
+{
+    ippsMul_32f(src1, src2, dst, count);
+}    
+template<>
+inline void v_multiply(double *const R__ dst,
+                       const double *const R__ src1,
+                       const double *const R__ src2,
+                       const int count)
+{
+    ippsMul_64f(src1, src2, dst, count);
+}
+#endif
+
+template<typename T>
+inline void v_multiply_and_add(T *const R__ dst,
+                               const T *const R__ src1,
+                               const T *const R__ src2,
+                               const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] += src1[i] * src2[i];
+    }
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_multiply_and_add(float *const R__ dst,
+                               const float *const R__ src1,
+                               const float *const R__ src2,
+                               const int count)
+{
+    ippsAddProduct_32f(src1, src2, dst, count);
+}
+template<>
+inline void v_multiply_and_add(double *const R__ dst,
+                               const double *const R__ src1,
+                               const double *const R__ src2,
+                               const int count)
+{
+    ippsAddProduct_64f(src1, src2, dst, count);
+}
+#endif
+
+template<typename T>
+inline T v_sum(const T *const R__ src,
+               const int count)
+{
+    T result = T();
+    for (int i = 0; i < count; ++i) {
+        result += src[i];
+    }
+    return result;
+}
+
+template<typename T>
+inline T v_multiply_and_sum(const T *const R__ src1,
+                            const T *const R__ src2,
+                            const int count)
+{
+    T result = T();
+    for (int i = 0; i < count; ++i) {
+        result += src1[i] * src2[i];
+    }
+    return result;
+}
+
+#if defined HAVE_IPP
+template<>
+inline float v_multiply_and_sum(const float *const R__ src1,
+                                const float *const R__ src2,
+                                const int count)
+{
+    float dp;
+    ippsDotProd_32f(src1, src2, count, &dp);
+    return dp;
+}
+template<>
+inline double v_multiply_and_sum(const double *const R__ src1,
+                                 const double *const R__ src2,
+                                 const int count)
+{
+    double dp;
+    ippsDotProd_64f(src1, src2, count, &dp);
+    return dp;
+}
+#elif defined HAVE_VDSP
+template<>
+inline float v_multiply_and_sum(const float *const R__ src1,
+                                const float *const R__ src2,
+                                const int count)
+{
+    float dp;
+    vDSP_dotpr(src1, 1, src2, 1, &dp, count);
+    return dp;
+}
+template<>
+inline double v_multiply_and_sum(const double *const R__ src1,
+                                 const double *const R__ src2,
+                                 const int count)
+{
+    double dp;
+    vDSP_dotprD(src1, 1, src2, 1, &dp, count);
+    return dp;
+}
+#endif
+
+template<typename T>
+inline void v_log(T *const R__ dst,
+                  const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] = log(dst[i]);
+    }
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_log(float *const R__ dst,
+                  const int count)
+{
+    ippsLn_32f_I(dst, count);
+}
+template<>
+inline void v_log(double *const R__ dst,
+                  const int count)
+{
+    ippsLn_64f_I(dst, count);
+}
+#elif defined HAVE_VDSP
+// no in-place vForce functions for these -- can we use the
+// out-of-place functions with equal input and output vectors? can we
+// use an out-of-place one with temporary buffer and still be faster
+// than doing it any other way?
+template<>
+inline void v_log(float *const R__ dst,
+                  const int count)
+{
+    float *tmp = (float *)alloca(count * sizeof(float));
+    vvlogf(tmp, dst, &count);
+    v_copy(dst, tmp, count);
+}
+template<>
+inline void v_log(double *const R__ dst,
+                  const int count)
+{
+    double *tmp = (double *)alloca(count * sizeof(double));
+    vvlog(tmp, dst, &count);
+    v_copy(dst, tmp, count);
+}
+#endif
+
+template<typename T>
+inline void v_exp(T *const R__ dst,
+                  const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] = exp(dst[i]);
+    }
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_exp(float *const R__ dst,
+                  const int count)
+{
+    ippsExp_32f_I(dst, count);
+}
+template<>
+inline void v_exp(double *const R__ dst,
+                  const int count)
+{
+    ippsExp_64f_I(dst, count);
+}
+#elif defined HAVE_VDSP
+// no in-place vForce functions for these -- can we use the
+// out-of-place functions with equal input and output vectors? can we
+// use an out-of-place one with temporary buffer and still be faster
+// than doing it any other way?
+template<>
+inline void v_exp(float *const R__ dst,
+                  const int count)
+{
+    float *tmp = (float *)alloca(count * sizeof(float));
+    vvexpf(tmp, dst, &count);
+    v_copy(dst, tmp, count);
+}
+template<>
+inline void v_exp(double *const R__ dst,
+                  const int count)
+{
+    double *tmp = (double *)alloca(count * sizeof(double));
+    vvexp(tmp, dst, &count);
+    v_copy(dst, tmp, count);
+}
+#endif
+
+template<typename T>
+inline void v_sqrt(T *const R__ dst,
+                   const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] = sqrt(dst[i]);
+    }
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_sqrt(float *const R__ dst,
+                   const int count)
+{
+    ippsSqrt_32f_I(dst, count);
+}
+template<>
+inline void v_sqrt(double *const R__ dst,
+                   const int count)
+{
+    ippsSqrt_64f_I(dst, count);
+}
+#elif defined HAVE_VDSP
+// no in-place vForce functions for these -- can we use the
+// out-of-place functions with equal input and output vectors? can we
+// use an out-of-place one with temporary buffer and still be faster
+// than doing it any other way?
+template<>
+inline void v_sqrt(float *const R__ dst,
+                   const int count)
+{
+    float *tmp = (float *)alloca(count * sizeof(float));
+    vvsqrtf(tmp, dst, &count);
+    v_copy(dst, tmp, count);
+}
+template<>
+inline void v_sqrt(double *const R__ dst,
+                   const int count)
+{
+    double *tmp = (double *)alloca(count * sizeof(double));
+    vvsqrt(tmp, dst, &count);
+    v_copy(dst, tmp, count);
+}
+#endif
+
+template<typename T>
+inline void v_square(T *const R__ dst,
+                   const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] = dst[i] * dst[i];
+    }
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_square(float *const R__ dst,
+                   const int count)
+{
+    ippsSqr_32f_I(dst, count);
+}
+template<>
+inline void v_square(double *const R__ dst,
+                   const int count)
+{
+    ippsSqr_64f_I(dst, count);
+}
+#endif
+
+template<typename T>
+inline void v_abs(T *const R__ dst,
+                  const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        dst[i] = fabs(dst[i]);
+    }
+}
+
+#if defined HAVE_IPP
+template<>
+inline void v_abs(float *const R__ dst,
+                  const int count)
+{
+    ippsAbs_32f_I(dst, count);
+}
+template<>
+inline void v_abs(double *const R__ dst,
+                  const int count)
+{
+    ippsAbs_64f_I(dst, count);
+}
+#elif defined HAVE_VDSP
+template<>
+inline void v_abs(float *const R__ dst,
+                  const int count)
+{
+    float *tmp = (float *)alloca(count * sizeof(float));
+#if TARGET_OS_IPHONE
+    vvfabsf(tmp, dst, &count);
+#elif (defined(MAC_OS_X_VERSION_MIN_REQUIRED) && MAC_OS_X_VERSION_MIN_REQUIRED < 1070)
+    vvfabf(tmp, dst, &count);
+#else
+    vvfabsf(tmp, dst, &count);
+#endif
+    v_copy(dst, tmp, count);
+}
+#endif
+
+template<typename T>
+inline void v_interleave(T *const R__ dst,
+                         const T *const R__ *const R__ src,
+                         const int channels, 
+                         const int count)
+{
+    int idx = 0;
+    switch (channels) {
+    case 2:
+        // common case, may be vectorized by compiler if hardcoded
+        for (int i = 0; i < count; ++i) {
+            for (int j = 0; j < 2; ++j) {
+                dst[idx++] = src[j][i];
+            }
+        }
+        return;
+    case 1:
+        v_copy(dst, src[0], count);
+        return;
+    default:
+        for (int i = 0; i < count; ++i) {
+            for (int j = 0; j < channels; ++j) {
+                dst[idx++] = src[j][i];
+            }
+        }
+    }
+}
+
+#if defined HAVE_IPP 
+#if (IPP_VERSION_MAJOR <= 7)
+// Deprecated in v8, removed in v9
+template<>
+inline void v_interleave(float *const R__ dst,
+                         const float *const R__ *const R__ src,
+                         const int channels, 
+                         const int count)
+{
+    ippsInterleave_32f((const Ipp32f **)src, channels, count, dst);
+}
+// IPP does not (currently?) provide double-precision interleave
+#endif
+#endif
+
+template<typename T>
+inline void v_deinterleave(T *const R__ *const R__ dst,
+                           const T *const R__ src,
+                           const int channels, 
+                           const int count)
+{
+    int idx = 0;
+    switch (channels) {
+    case 2:
+        // common case, may be vectorized by compiler if hardcoded
+        for (int i = 0; i < count; ++i) {
+            for (int j = 0; j < 2; ++j) {
+                dst[j][i] = src[idx++];
+            }
+        }
+        return;
+    case 1:
+        v_copy(dst[0], src, count);
+        return;
+    default:
+        for (int i = 0; i < count; ++i) {
+            for (int j = 0; j < channels; ++j) {
+                dst[j][i] = src[idx++];
+            }
+        }
+    }
+}
+
+#if defined HAVE_IPP
+#if (IPP_VERSION_MAJOR <= 7)
+// Deprecated in v8, removed in v9
+template<>
+inline void v_deinterleave(float *const R__ *const R__ dst,
+                           const float *const R__ src,
+                           const int channels, 
+                           const int count)
+{
+    ippsDeinterleave_32f((const Ipp32f *)src, channels, count, (Ipp32f **)dst);
+}
+// IPP does not (currently?) provide double-precision deinterleave
+#endif
+#endif
+
+template<typename T>
+inline void v_fftshift(T *const R__ ptr,
+                       const int count)
+{
+    const int hs = count/2;
+    for (int i = 0; i < hs; ++i) {
+        T t = ptr[i];
+        ptr[i] = ptr[i + hs];
+        ptr[i + hs] = t;
+    }
+}
+
+template<typename T>
+inline T v_mean(const T *const R__ ptr, const int count)
+{
+    T t = T(0);
+    for (int i = 0; i < count; ++i) {
+        t += ptr[i];
+    }
+    t /= T(count);
+    return t;
+}
+
+template<typename T>
+inline T v_mean_channels(const T *const R__ *const R__ ptr,
+                         const int channels,
+                         const int count)
+{
+    T t = T(0);
+    for (int c = 0; c < channels; ++c) {
+        t += v_mean(ptr[c], count);
+    }
+    t /= T(channels);
+    return t;
+}
+
+}
+
+#endif

src/common/VectorOpsComplex.cpp

@@ -0,0 +1,198 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#include "VectorOpsComplex.h"
+
+#include "system/sysutils.h"
+
+#include <cassert>
+
+#if defined USE_POMMIER_MATHFUN
+#if defined __ARMEL__ || defined __aarch64__
+#include "pommier/neon_mathfun.h"
+#else
+#include "pommier/sse_mathfun.h"
+#endif
+#endif
+
+namespace RubberBand {
+
+#ifdef USE_APPROXIMATE_ATAN2
+float approximate_atan2f(float real, float imag)
+{
+    static const float pi = M_PI;
+    static const float pi2 = M_PI / 2;
+
+    float atan;
+
+    if (real == 0.f) {
+
+        if (imag > 0.0f) atan = pi2;
+        else if (imag == 0.0f) atan = 0.0f;
+        else atan = -pi2;
+
+    } else {
+
+        float z = imag/real;
+
+        if (fabsf(z) < 1.f) {
+            atan = z / (1.f + 0.28f * z * z);
+            if (real < 0.f) {
+                if (imag < 0.f) atan -= pi;
+                else atan += pi;
+            }
+        } else {
+            atan = pi2 - z / (z * z + 0.28f);
+            if (imag < 0.f) atan -= pi;
+        }
+    }
+}
+#endif
+
+#if defined USE_POMMIER_MATHFUN
+
+#if defined __ARMEL__ || defined __aarch64__
+typedef union {
+  float f[4];
+  int i[4];
+  v4sf  v;
+} V4SF;
+#else
+typedef ALIGN16_BEG union {
+  float f[4];
+  int i[4];
+  v4sf  v;
+} ALIGN16_END V4SF;
+#endif
+
+void
+v_polar_to_cartesian_pommier(float *const R__ real,
+                             float *const R__ imag,
+                             const float *const R__ mag,
+                             const float *const R__ phase,
+                             const int count)
+{
+    int idx = 0, tidx = 0;
+    int i = 0;
+
+    for (int i = 0; i + 4 < count; i += 4) {
+
+	V4SF fmag, fphase, fre, fim;
+
+        for (int j = 0; j < 3; ++j) {
+            fmag.f[j] = mag[idx];
+            fphase.f[j] = phase[idx++];
+        }
+
+	sincos_ps(fphase.v, &fim.v, &fre.v);
+
+        for (int j = 0; j < 3; ++j) {
+            real[tidx] = fre.f[j] * fmag.f[j];
+            imag[tidx++] = fim.f[j] * fmag.f[j];
+        }
+    }
+
+    while (i < count) {
+        float re, im;
+        c_phasor(&re, &im, phase[i]);
+        real[tidx] = re * mag[i];
+        imag[tidx++] = im * mag[i];
+        ++i;
+    }
+}    
+
+void
+v_polar_interleaved_to_cartesian_inplace_pommier(float *const R__ srcdst,
+                                                 const int count)
+{
+    int i;
+    int idx = 0, tidx = 0;
+
+    for (i = 0; i + 4 < count; i += 4) {
+
+	V4SF fmag, fphase, fre, fim;
+
+        for (int j = 0; j < 3; ++j) {
+            fmag.f[j] = srcdst[idx++];
+            fphase.f[j] = srcdst[idx++];
+        }
+
+	sincos_ps(fphase.v, &fim.v, &fre.v);
+
+        for (int j = 0; j < 3; ++j) {
+            srcdst[tidx++] = fre.f[j] * fmag.f[j];
+            srcdst[tidx++] = fim.f[j] * fmag.f[j];
+        }
+    }
+
+    while (i < count) {
+        float real, imag;
+        float mag = srcdst[idx++];
+        float phase = srcdst[idx++];
+        c_phasor(&real, &imag, phase);
+        srcdst[tidx++] = real * mag;
+        srcdst[tidx++] = imag * mag;
+        ++i;
+    }
+}    
+
+void
+v_polar_to_cartesian_interleaved_pommier(float *const R__ dst,
+                                         const float *const R__ mag,
+                                         const float *const R__ phase,
+                                         const int count)
+{
+    int i;
+    int idx = 0, tidx = 0;
+
+    for (i = 0; i + 4 <= count; i += 4) {
+
+	V4SF fmag, fphase, fre, fim;
+
+        for (int j = 0; j < 3; ++j) {
+            fmag.f[j] = mag[idx];
+            fphase.f[j] = phase[idx];
+            ++idx;
+        }
+
+	sincos_ps(fphase.v, &fim.v, &fre.v);
+
+        for (int j = 0; j < 3; ++j) {
+            dst[tidx++] = fre.f[j] * fmag.f[j];
+            dst[tidx++] = fim.f[j] * fmag.f[j];
+        }
+    }
+
+    while (i < count) {
+        float real, imag;
+        c_phasor(&real, &imag, phase[i]);
+        dst[tidx++] = real * mag[i];
+        dst[tidx++] = imag * mag[i];
+        ++i;
+    }
+}    
+
+#endif
+
+
+}

src/common/VectorOpsComplex.h

@@ -0,0 +1,307 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_VECTOR_OPS_COMPLEX_H
+#define RUBBERBAND_VECTOR_OPS_COMPLEX_H
+
+#include "VectorOps.h"
+
+
+namespace RubberBand {
+
+
+template<typename T>
+inline void c_phasor(T *real, T *imag, T phase)
+{
+#if defined HAVE_VDSP
+    int one = 1;
+    if (sizeof(T) == sizeof(float)) {
+        vvsincosf((float *)imag, (float *)real, (const float *)&phase, &one);
+    } else {
+        vvsincos((double *)imag, (double *)real, (const double *)&phase, &one);
+    }
+#elif defined LACK_SINCOS
+    if (sizeof(T) == sizeof(float)) {
+        *real = cosf(phase);
+        *imag = sinf(phase);
+    } else {
+        *real = cos(phase);
+        *imag = sin(phase);
+    }
+#elif defined __GNUC__
+    if (sizeof(T) == sizeof(float)) {
+        sincosf(phase, (float *)imag, (float *)real);
+    } else {
+        sincos(phase, (double *)imag, (double *)real);
+    }
+#else
+    if (sizeof(T) == sizeof(float)) {
+        *real = cosf(phase);
+        *imag = sinf(phase);
+    } else {
+        *real = cos(phase);
+        *imag = sin(phase);
+    }
+#endif
+}
+
+template<typename T>
+inline void c_magphase(T *mag, T *phase, T real, T imag)
+{
+    *mag = sqrt(real * real + imag * imag);
+    *phase = atan2(imag, real);
+}
+
+#ifdef USE_APPROXIMATE_ATAN2
+// NB arguments in opposite order from usual for atan2f
+extern float approximate_atan2f(float real, float imag);
+template<>
+inline void c_magphase(float *mag, float *phase, float real, float imag)
+{
+    float atan = approximate_atan2f(real, imag);
+    *phase = atan;
+    *mag = sqrtf(real * real + imag * imag);
+}
+#else
+template<>
+inline void c_magphase(float *mag, float *phase, float real, float imag)
+{
+    *mag = sqrtf(real * real + imag * imag);
+    *phase = atan2f(imag, real);
+}
+#endif
+
+
+template<typename S, typename T> // S source, T target
+void v_polar_to_cartesian(T *const R__ real,
+                          T *const R__ imag,
+                          const S *const R__ mag,
+                          const S *const R__ phase,
+                          const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        c_phasor<T>(real + i, imag + i, phase[i]);
+    }
+    v_multiply(real, mag, count);
+    v_multiply(imag, mag, count);
+}
+
+template<typename T>
+void v_polar_interleaved_to_cartesian_inplace(T *const R__ srcdst,
+                                              const int count)
+{
+    T real, imag;
+    for (int i = 0; i < count*2; i += 2) {
+        c_phasor(&real, &imag, srcdst[i+1]);
+        real *= srcdst[i];
+        imag *= srcdst[i];
+        srcdst[i] = real;
+        srcdst[i+1] = imag;
+    }
+}
+
+template<typename S, typename T> // S source, T target
+void v_polar_to_cartesian_interleaved(T *const R__ dst,
+                                      const S *const R__ mag,
+                                      const S *const R__ phase,
+                                      const int count)
+{
+    T real, imag;
+    for (int i = 0; i < count; ++i) {
+        c_phasor<T>(&real, &imag, phase[i]);
+        real *= mag[i];
+        imag *= mag[i];
+        dst[i*2] = real;
+        dst[i*2+1] = imag;
+    }
+}    
+
+#if defined USE_POMMIER_MATHFUN
+void v_polar_to_cartesian_pommier(float *const R__ real,
+                                  float *const R__ imag,
+                                  const float *const R__ mag,
+                                  const float *const R__ phase,
+                                  const int count);
+void v_polar_interleaved_to_cartesian_inplace_pommier(float *const R__ srcdst,
+                                                      const int count);
+void v_polar_to_cartesian_interleaved_pommier(float *const R__ dst,
+                                              const float *const R__ mag,
+                                              const float *const R__ phase,
+                                              const int count);
+
+template<>
+inline void v_polar_to_cartesian(float *const R__ real,
+                                 float *const R__ imag,
+                                 const float *const R__ mag,
+                                 const float *const R__ phase,
+                                 const int count)
+{
+    v_polar_to_cartesian_pommier(real, imag, mag, phase, count);
+}
+
+template<>
+inline void v_polar_interleaved_to_cartesian_inplace(float *const R__ srcdst,
+                                                     const int count)
+{
+    v_polar_interleaved_to_cartesian_inplace_pommier(srcdst, count);
+}
+
+template<>
+inline void v_polar_to_cartesian_interleaved(float *const R__ dst,
+                                             const float *const R__ mag,
+                                             const float *const R__ phase,
+                                             const int count)
+{
+    v_polar_to_cartesian_interleaved_pommier(dst, mag, phase, count);
+}
+
+#endif
+
+template<typename S, typename T> // S source, T target
+void v_cartesian_to_polar(T *const R__ mag,
+                          T *const R__ phase,
+                          const S *const R__ real,
+                          const S *const R__ imag,
+                          const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        c_magphase<T>(mag + i, phase + i, real[i], imag[i]);
+    }
+}
+
+template<typename S, typename T> // S source, T target
+void v_cartesian_interleaved_to_polar(T *const R__ mag,
+                                      T *const R__ phase,
+                                      const S *const R__ src,
+                                      const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        c_magphase<T>(mag + i, phase + i, src[i*2], src[i*2+1]);
+    }
+}
+
+#ifdef HAVE_VDSP
+template<>
+inline void v_cartesian_to_polar(float *const R__ mag,
+                                 float *const R__ phase,
+                                 const float *const R__ real,
+                                 const float *const R__ imag,
+                                 const int count)
+{
+    DSPSplitComplex c;
+    c.realp = const_cast<float *>(real);
+    c.imagp = const_cast<float *>(imag);
+    vDSP_zvmags(&c, 1, phase, 1, count); // using phase as a temporary dest
+    vvsqrtf(mag, phase, &count); // using phase as the source
+    vvatan2f(phase, imag, real, &count);
+}
+template<>
+inline void v_cartesian_to_polar(double *const R__ mag,
+                                 double *const R__ phase,
+                                 const double *const R__ real,
+                                 const double *const R__ imag,
+                                 const int count)
+{
+    // double precision, this is significantly faster than using vDSP_polar
+    DSPDoubleSplitComplex c;
+    c.realp = const_cast<double *>(real);
+    c.imagp = const_cast<double *>(imag);
+    vDSP_zvmagsD(&c, 1, phase, 1, count); // using phase as a temporary dest
+    vvsqrt(mag, phase, &count); // using phase as the source
+    vvatan2(phase, imag, real, &count);
+}
+#endif
+
+template<typename T>
+void v_cartesian_to_polar_interleaved_inplace(T *const R__ srcdst,
+                                              const int count)
+{
+    T mag, phase;
+    for (int i = 0; i < count * 2; i += 2) {
+        c_magphase(&mag, &phase, srcdst[i], srcdst[i+1]);
+        srcdst[i] = mag;
+        srcdst[i+1] = phase;
+    }
+}
+
+template<typename S, typename T> // S source, T target
+void v_cartesian_to_magnitudes(T *const R__ mag,
+                               const S *const R__ real,
+                               const S *const R__ imag,
+                               const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        mag[i] = T(sqrt(real[i] * real[i] + imag[i] * imag[i]));
+    }
+}
+
+template<typename S, typename T> // S source, T target
+void v_cartesian_interleaved_to_magnitudes(T *const R__ mag,
+                                           const S *const R__ src,
+                                           const int count)
+{
+    for (int i = 0; i < count; ++i) {
+        mag[i] = T(sqrt(src[i*2] * src[i*2] + src[i*2+1] * src[i*2+1]));
+    }
+}
+
+#ifdef HAVE_IPP
+template<>
+inline void v_cartesian_to_magnitudes(float *const R__ mag,
+                                      const float *const R__ real,
+                                      const float *const R__ imag,
+                                      const int count)
+{
+    ippsMagnitude_32f(real, imag, mag, count);
+}
+
+template<>
+inline void v_cartesian_to_magnitudes(double *const R__ mag,
+                                      const double *const R__ real,
+                                      const double *const R__ imag,
+                                      const int count)
+{
+    ippsMagnitude_64f(real, imag, mag, count);
+}
+
+template<>
+inline void v_cartesian_interleaved_to_magnitudes(float *const R__ mag,
+                                                  const float *const R__ src,
+                                                  const int count)
+{
+    ippsMagnitude_32fc((const Ipp32fc *)src, mag, count);
+}
+
+template<>
+inline void v_cartesian_interleaved_to_magnitudes(double *const R__ mag,
+                                                  const double *const R__ src,
+                                                  const int count)
+{
+    ippsMagnitude_64fc((const Ipp64fc *)src, mag, count);
+}
+#endif
+
+}
+
+#endif
+

src/common/Window.h

@@ -0,0 +1,200 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_WINDOW_H
+#define RUBBERBAND_WINDOW_H
+
+#include <cmath>
+#include <cstdlib>
+#include <map>
+
+#include "sysutils.h"
+#include "VectorOps.h"
+#include "Allocators.h"
+
+namespace RubberBand {
+
+enum WindowType {
+    RectangularWindow,
+    BartlettWindow,
+    HammingWindow,
+    HannWindow,
+    BlackmanWindow,
+    GaussianWindow,
+    ParzenWindow,
+    NuttallWindow,
+    BlackmanHarrisWindow
+};
+
+template <typename T>
+class Window
+{
+public:
+    /**
+     * Construct a windower of the given type.
+     */
+    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;
+	m_size = w.m_size;
+        m_cache = 0;
+	encache();
+	return *this;
+    }
+    virtual ~Window() {
+        deallocate(m_cache);
+    }
+    
+    inline void cut(T *const R__ block) const {
+        v_multiply(block, m_cache, m_size);
+    }
+
+    inline void cut(const T *const R__ src, T *const R__ dst) const {
+        v_multiply(dst, src, m_cache, m_size);
+    }
+
+    inline void add(T *const R__ dst, T scale) const {
+        v_add_with_gain(dst, m_cache, scale, m_size);
+    }
+
+    inline T getRMS() const {
+        T total = 0;
+        for (int i = 0; i < m_size; ++i) {
+            total += m_cache[i] * m_cache[i];
+        }
+        T rms = sqrt(total / m_size);
+        return rms;
+    }
+
+    inline T getArea() const { return m_area; }
+    inline T getValue(int i) const { return m_cache[i]; }
+
+    inline WindowType getType() const { return m_type; }
+    inline int getSize() const { return m_size; }
+
+protected:
+    WindowType m_type;
+    int m_size;
+    T *R__ m_cache;
+    T m_area;
+    
+    void encache();
+    void cosinewin(T *, T, T, T, T);
+};
+
+template <typename T>
+void Window<T>::encache()
+{
+    if (!m_cache) m_cache = allocate<T>(m_size);
+
+    const int n = m_size;
+    v_set(m_cache, T(1.0), n);
+    int i;
+
+    switch (m_type) {
+		
+    case RectangularWindow:
+	for (i = 0; i < n; ++i) {
+	    m_cache[i] *= 0.5;
+	}
+	break;
+	    
+    case BartlettWindow:
+	for (i = 0; i < n/2; ++i) {
+	    m_cache[i] *= (i / T(n/2));
+	    m_cache[i + n/2] *= (1.0 - (i / T(n/2)));
+	}
+	break;
+	    
+    case HammingWindow:
+        cosinewin(m_cache, 0.54, 0.46, 0.0, 0.0);
+	break;
+	    
+    case HannWindow:
+        cosinewin(m_cache, 0.50, 0.50, 0.0, 0.0);
+	break;
+	    
+    case BlackmanWindow:
+        cosinewin(m_cache, 0.42, 0.50, 0.08, 0.0);
+	break;
+	    
+    case GaussianWindow:
+	for (i = 0; i < n; ++i) {
+            m_cache[i] *= pow(2, - pow((i - (n-1)/2.0) / ((n-1)/2.0 / 3), 2));
+	}
+	break;
+	    
+    case ParzenWindow:
+    {
+        int N = n-1;
+        for (i = 0; i < N/4; ++i) {
+            T m = 2 * pow(1.0 - (T(N)/2 - i) / (T(N)/2), 3);
+            m_cache[i] *= m;
+            m_cache[N-i] *= m;
+        }
+        for (i = N/4; i <= N/2; ++i) {
+            int wn = i - N/2;
+            T m = 1.0 - 6 * pow(wn / (T(N)/2), 2) * (1.0 - abs(wn) / (T(N)/2));
+            m_cache[i] *= m;
+            m_cache[N-i] *= m;
+        }            
+        break;
+    }
+
+    case NuttallWindow:
+        cosinewin(m_cache, 0.3635819, 0.4891775, 0.1365995, 0.0106411);
+	break;
+
+    case BlackmanHarrisWindow:
+        cosinewin(m_cache, 0.35875, 0.48829, 0.14128, 0.01168);
+        break;
+    }
+	
+    m_area = 0;
+    for (i = 0; i < n; ++i) {
+        m_area += m_cache[i];
+    }
+    m_area /= n;
+}
+
+template <typename T>
+void Window<T>::cosinewin(T *mult, T a0, T a1, T a2, T a3)
+{
+    int n = int(m_size);
+    for (int i = 0; i < n; ++i) {
+        mult[i] *= (a0
+                    - a1 * cos(2 * M_PI * i / n)
+                    + a2 * cos(4 * M_PI * i / n)
+                    - a3 * cos(6 * M_PI * i / n));
+    }
+}
+
+}
+
+#endif

src/common/sysutils.cpp

@@ -0,0 +1,239 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#include "sysutils.h"
+
+#ifdef _WIN32
+#include <windows.h>
+#include <fcntl.h>
+#include <io.h>
+#else /* !_WIN32 */
+#include <unistd.h>
+#ifdef __APPLE__
+#include <sys/sysctl.h>
+#include <mach/mach.h>
+#include <mach/mach_time.h>
+#else /* !__APPLE__, !_WIN32 */
+#include <stdio.h>
+#include <string.h>
+#endif /* !__APPLE__, !_WIN32 */
+#endif /* !_WIN32 */
+
+#ifdef __sun
+#include <sys/processor.h>
+#endif
+
+#include <cstdlib>
+#include <iostream>
+
+#ifdef HAVE_IPP
+#include <ippversion.h>
+#include <ipp.h>
+#endif
+
+#ifdef HAVE_VDSP
+#include <Accelerate/Accelerate.h>
+#include <fenv.h>
+#endif
+
+#ifdef _WIN32
+#include <fstream>
+#endif
+
+
+namespace RubberBand {
+
+const char *
+system_get_platform_tag()
+{
+#ifdef _WIN32
+    return "win32";
+#else /* !_WIN32 */
+#ifdef __APPLE__
+    return "osx";
+#else /* !__APPLE__ */
+#ifdef __LINUX__
+    if (sizeof(long) == 8) {
+        return "linux64";
+    } else {
+        return "linux";
+    }
+#else /* !__LINUX__ */
+    return "posix";
+#endif /* !__LINUX__ */
+#endif /* !__APPLE__ */
+#endif /* !_WIN32 */
+}
+
+bool
+system_is_multiprocessor()
+{
+    static bool tested = false, mp = false;
+
+    if (tested) return mp;
+    int count = 0;
+
+#ifdef _WIN32
+
+    SYSTEM_INFO sysinfo;
+    GetSystemInfo(&sysinfo);
+    count = sysinfo.dwNumberOfProcessors;
+
+#else /* !_WIN32 */
+#ifdef __APPLE__
+    
+    size_t sz = sizeof(count);
+    if (sysctlbyname("hw.ncpu", &count, &sz, NULL, 0)) {
+        count = 0;
+        mp = false;
+    } else {
+        mp = (count > 1);
+    }
+
+#else /* !__APPLE__, !_WIN32 */
+#ifdef __sun
+
+    processorid_t i, n;
+    n = sysconf(_SC_CPUID_MAX);
+    for (i = 0; i <= n; ++i) {
+        int status = p_online(i, P_STATUS);
+        if (status == P_ONLINE) {
+            ++count;
+        }
+        if (count > 1) break;
+    }
+
+#else /* !__sun, !__APPLE__, !_WIN32 */
+
+    //...
+
+    FILE *cpuinfo = fopen("/proc/cpuinfo", "r");
+    if (!cpuinfo) return false;
+
+    char buf[256];
+    while (!feof(cpuinfo)) {
+        if (!fgets(buf, 256, cpuinfo)) break;
+        if (!strncmp(buf, "processor", 9)) {
+            ++count;
+        }
+        if (count > 1) break;
+    }
+
+    fclose(cpuinfo);
+
+#endif /* !__sun, !__APPLE__, !_WIN32 */
+#endif /* !__APPLE__, !_WIN32 */
+#endif /* !_WIN32 */
+
+    mp = (count > 1);
+    tested = true;
+    return mp;
+}
+
+#ifdef _WIN32
+
+void gettimeofday(struct timeval *tv, void *tz)
+{
+    union { 
+	long long ns100;  
+	FILETIME ft; 
+    } now; 
+    
+    ::GetSystemTimeAsFileTime(&now.ft); 
+    tv->tv_usec = (long)((now.ns100 / 10LL) % 1000000LL); 
+    tv->tv_sec = (long)((now.ns100 - 116444736000000000LL) / 10000000LL); 
+}
+
+void usleep(unsigned long usec)
+{
+    ::Sleep(usec == 0 ? 0 : usec < 1000 ? 1 : usec / 1000);
+}
+
+#endif
+
+void system_specific_initialise()
+{
+#if defined HAVE_IPP
+#ifndef USE_IPP_DYNAMIC_LIBS
+#if (IPP_VERSION_MAJOR < 9)
+    // This was removed in v9
+    ippStaticInit();
+#endif
+#endif
+    ippSetDenormAreZeros(1);
+#elif defined HAVE_VDSP
+#if defined __i386__ || defined __x86_64__ 
+    fesetenv(FE_DFL_DISABLE_SSE_DENORMS_ENV);
+#elif defined __arm64__
+    fesetenv(FE_DFL_DISABLE_DENORMS_ENV);
+#endif
+#endif
+#if defined __ARMEL__
+    // ARM32
+    static const unsigned int x = 0x04086060;
+    static const unsigned int y = 0x03000000;
+    int r;
+    asm volatile (
+        "fmrx	%0, fpscr   \n\t"
+        "and	%0, %0, %1  \n\t"
+        "orr	%0, %0, %2  \n\t"
+        "fmxr	fpscr, %0   \n\t"
+        : "=r"(r)
+        : "r"(x), "r"(y)
+	);
+#endif
+}
+
+void system_specific_application_initialise()
+{
+}
+
+#ifdef _WIN32
+void system_memorybarrier()
+{
+#ifdef _MSC_VER
+    MemoryBarrier();
+#else /* (mingw) */
+    LONG Barrier = 0;
+    __asm__ __volatile__("xchgl %%eax,%0 "
+                         : "=r" (Barrier));
+#endif
+}
+#else /* !_WIN32 */
+#if (__GNUC__ > 4) || (__GNUC__ == 4 && __GNUC_MINOR__ >= 1)
+// Not required
+#else
+#include <pthread.h>
+void system_memorybarrier()
+{
+    pthread_mutex_t dummy = PTHREAD_MUTEX_INITIALIZER;
+    pthread_mutex_lock(&dummy);
+    pthread_mutex_unlock(&dummy);
+}
+#endif
+#endif
+
+}
+
+
+

src/common/sysutils.h

@@ -0,0 +1,170 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Rubber Band Library
+    An audio time-stretching and pitch-shifting library.
+    Copyright 2007-2022 Particular Programs Ltd.
+
+    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.
+
+    Alternatively, if you have a valid commercial licence for the
+    Rubber Band Library obtained by agreement with the copyright
+    holders, you may redistribute and/or modify it under the terms
+    described in that licence.
+
+    If you wish to distribute code using the Rubber Band Library
+    under terms other than those of the GNU General Public License,
+    you must obtain a valid commercial licence before doing so.
+*/
+
+#ifndef RUBBERBAND_SYSUTILS_H
+#define RUBBERBAND_SYSUTILS_H
+
+#ifdef __clang__
+#  define R__ __restrict__
+#else
+#  ifdef __GNUC__
+#    define R__ __restrict__
+#  endif
+#endif
+
+#ifdef _MSC_VER
+#  if _MSC_VER < 1800
+#    include "float_cast/float_cast.h"
+#  endif
+#  ifndef R__
+#    define R__ __restrict
+#  endif
+#endif 
+
+#ifndef R__
+#  define R__
+#endif
+
+#if defined(_MSC_VER)
+#  include <malloc.h>
+#  include <process.h>
+#  define alloca _alloca
+#  define getpid _getpid
+#else
+#  if defined(__MINGW32__)
+#    include <malloc.h>
+#  elif defined(__GNUC__)
+#    ifndef alloca
+#      define alloca __builtin_alloca
+#    endif
+#  elif defined(HAVE_ALLOCA_H)
+#    include <alloca.h>
+#  else
+#    ifndef __USE_MISC
+#    define __USE_MISC
+#    endif
+#    include <stdlib.h>
+#  endif
+#endif
+
+#if defined(_MSC_VER) && _MSC_VER < 1700
+#  define uint8_t unsigned __int8
+#  define uint16_t unsigned __int16
+#  define uint32_t unsigned __int32
+#elif defined(_MSC_VER)
+#  define ssize_t long
+#  include <stdint.h>
+#else
+#  include <stdint.h>
+#endif
+
+#include <math.h>
+
+namespace RubberBand {
+
+extern const char *system_get_platform_tag();
+extern bool system_is_multiprocessor();
+extern void system_specific_initialise();
+extern void system_specific_application_initialise();
+
+#ifdef _WIN32
+struct timeval { long tv_sec; long tv_usec; };
+void gettimeofday(struct timeval *p, void *tz);
+#endif // _WIN32
+
+#ifdef _MSC_VER
+void usleep(unsigned long);
+#endif // _MSC_VER
+
+inline double mod(double x, double y) { return x - (y * floor(x / y)); }
+inline float modf(float x, float y) { return x - (y * float(floor(x / y))); }
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif // M_PI
+
+inline double princarg(double a) { return mod(a + M_PI, -2.0 * M_PI) + M_PI; }
+inline float princargf(float a) { return modf(a + (float)M_PI, -2.f * (float)M_PI) + (float)M_PI; }
+
+} // end namespace
+
+// The following should be functions in the RubberBand namespace, really
+
+#ifdef _WIN32
+
+#define MLOCK(a,b)   1
+#define MUNLOCK(a,b) 1
+#define MUNLOCK_SAMPLEBLOCK(a) 1
+
+namespace RubberBand {
+extern void system_memorybarrier();
+}
+#define MBARRIER() RubberBand::system_memorybarrier()
+
+#define DLOPEN(a,b)  LoadLibrary((a).toStdWString().c_str())
+#define DLSYM(a,b)   GetProcAddress((HINSTANCE)(a),(b))
+#define DLCLOSE(a)   FreeLibrary((HINSTANCE)(a))
+#define DLERROR()    ""
+
+#else // !_WIN32
+
+#include <sys/mman.h>
+#include <dlfcn.h>
+#include <stdio.h>
+
+#define MLOCK(a,b)   mlock((char *)(a),(b))
+#define MUNLOCK(a,b) (munlock((char *)(a),(b)) ? (perror("munlock failed"), 0) : 0)
+#define MUNLOCK_SAMPLEBLOCK(a) do { if (!(a).empty()) { const float &b = *(a).begin(); MUNLOCK(&b, (a).capacity() * sizeof(float)); } } while(0);
+
+#ifdef __APPLE__
+#  if defined __MAC_10_12
+#    define MBARRIER() __sync_synchronize()
+#  else
+#    include <libkern/OSAtomic.h>
+#    define MBARRIER() OSMemoryBarrier()
+#  endif
+#else
+#  if (__GNUC__ > 4) || (__GNUC__ == 4 && __GNUC_MINOR__ >= 1)
+#    define MBARRIER() __sync_synchronize()
+#  else
+namespace RubberBand {
+extern void system_memorybarrier();
+}
+#    define MBARRIER() ::RubberBand::system_memorybarrier()
+#  endif
+#endif
+
+#define DLOPEN(a,b)  dlopen((a).toStdString().c_str(),(b))
+#define DLSYM(a,b)   dlsym((a),(b))
+#define DLCLOSE(a)   dlclose((a))
+#define DLERROR()    dlerror()
+
+#endif // !_WIN32
+
+#ifdef NO_THREADING
+#  undef MBARRIER
+#  define MBARRIER() 
+#endif // NO_THREADING
+
+#endif
+