Merge pull request #287 from jpcima/wavetables-hq

Wavetables hq

src/sfizz/Defaults.h

@@ -67,6 +67,8 @@ namespace Default
     constexpr Range<int> oscillatorMultiRange { 1, config::oscillatorsPerVoice };
     constexpr float oscillatorDetune { 0 };
     constexpr Range<float> oscillatorDetuneRange { -9600, 9600 };
+    constexpr int oscillatorQuality { 1 };
+    constexpr Range<int> oscillatorQualityRange { 0, 3 };
 
     // Instrument setting: voice lifecycle
 	constexpr uint32_t group { 0 };

src/sfizz/Interpolators.h

@@ -9,6 +9,8 @@
 namespace sfz {
 
 enum InterpolatorModel : int {
+    // a nearest interpolator
+    kInterpolatorNearest,
     // a linear interpolator
     kInterpolatorLinear,
     // a Hermite 3rd order interpolator

src/sfizz/Interpolators.hpp

@@ -19,6 +19,19 @@ inline R interpolate(const R* values, R coeff)
     return Interpolator<M, R>::process(values, coeff);
 }
 
+//------------------------------------------------------------------------------
+// Nearest
+
+template <class R>
+class Interpolator<kInterpolatorNearest, R>
+{
+public:
+    static inline R process(const R* values, R coeff)
+    {
+        return values[coeff > static_cast<R>(0.5)];
+    }
+};
+
 //------------------------------------------------------------------------------
 // Linear
 

src/sfizz/Region.cpp

@@ -135,6 +135,12 @@ bool sfz::Region::parseOpcode(const Opcode& rawOpcode)
     case hash("oscillator_detune"):
         setValueFromOpcode(opcode, oscillatorDetune, Default::oscillatorDetuneRange);
         break;
+    case hash("oscillator_quality"):
+        if (opcode.value == "-1")
+            oscillatorQuality.reset();
+        else
+            setValueFromOpcode(opcode, oscillatorQuality, Default::oscillatorQualityRange);
+        break;
 
     // Instrument settings: voice lifecycle
     case hash("group"): // also polyphony_group

src/sfizz/Region.h

@@ -276,6 +276,7 @@ struct Region {
     bool oscillator = false;
     int oscillatorMulti = Default::oscillatorMulti;
     float oscillatorDetune = Default::oscillatorDetune;
+    absl::optional<int> oscillatorQuality;
 
     // Instrument settings: voice lifecycle
     uint32_t group { Default::group }; // group

src/sfizz/Voice.cpp

@@ -65,16 +65,22 @@ void sfz::Voice::startVoice(Region* region, int delay, int number, float value,
             wave = resources.wavePool.getWaveSaw();
             break;
         }
+        const float phase = region->getPhase();
+        const int quality = region->oscillatorQuality.value_or(Default::oscillatorQuality);
         for (WavetableOscillator& osc : waveOscillators) {
             osc.setWavetable(wave);
-            osc.setPhase(region->getPhase());
+            osc.setPhase(phase);
+            osc.setQuality(quality);
         }
         setupOscillatorUnison();
     } else if (region->oscillator) {
         const WavetableMulti* wave = resources.wavePool.getFileWave(region->sampleId.filename());
+        const float phase = region->getPhase();
+        const int quality = region->oscillatorQuality.value_or(Default::oscillatorQuality);
         for (WavetableOscillator& osc : waveOscillators) {
             osc.setWavetable(wave);
-            osc.setPhase(region->getPhase());
+            osc.setPhase(phase);
+            osc.setQuality(quality);
         }
         setupOscillatorUnison();
     } else {

src/sfizz/Wavetables.cpp

@@ -6,6 +6,7 @@
 
 #include "Wavetables.h"
 #include "FilePool.h"
+#include "Interpolators.h"
 #include "MathHelpers.h"
 #include "absl/meta/type_traits.h"
 #include <kiss_fftr.h>
@@ -35,7 +36,8 @@ void WavetableOscillator::setPhase(float phase)
     _phase = phase;
 }
 
-void WavetableOscillator::process(float frequency, float detuneRatio, float* output, unsigned nframes)
+template <InterpolatorModel M>
+void WavetableOscillator::processSingle(float frequency, float detuneRatio, float* output, unsigned nframes)
 {
     float phase = _phase;
     float phaseInc = frequency * (detuneRatio * _sampleInterval);
@@ -48,7 +50,7 @@ void WavetableOscillator::process(float frequency, float detuneRatio, float* out
         float position = phase * tableSize;
         unsigned index = static_cast<unsigned>(position);
         float frac = position - index;
-        output[i] = interpolate(&table[index], frac);
+        output[i] = interpolate<M>(&table[index], frac);
 
         phase += phaseInc;
         phase -= static_cast<int>(phase);
@@ -57,7 +59,8 @@ void WavetableOscillator::process(float frequency, float detuneRatio, float* out
     _phase = phase;
 }
 
-void WavetableOscillator::processModulated(const float* frequencies, float detuneRatio, float* output, unsigned nframes)
+template <InterpolatorModel M>
+void WavetableOscillator::processModulatedSingle(const float* frequencies, float detuneRatio, float* output, unsigned nframes)
 {
     float phase = _phase;
     float sampleInterval = _sampleInterval;
@@ -73,7 +76,7 @@ void WavetableOscillator::processModulated(const float* frequencies, float detun
         float position = phase * tableSize;
         unsigned index = static_cast<unsigned>(position);
         float frac = position - index;
-        output[i] = interpolate(&table[index], frac);
+        output[i] = interpolate<M>(&table[index], frac);
 
         phase += phaseInc;
         phase -= static_cast<int>(phase);
@@ -82,9 +85,98 @@ void WavetableOscillator::processModulated(const float* frequencies, float detun
     _phase = phase;
 }
 
-float WavetableOscillator::interpolate(const float* x, float delta)
+template <InterpolatorModel M>
+void WavetableOscillator::processDual(float frequency, float detuneRatio, float* output, unsigned nframes)
 {
-    return x[0] + delta * (x[1] - x[0]);
+    float phase = _phase;
+    float phaseInc = frequency * (detuneRatio * _sampleInterval);
+
+    const WavetableMulti& multi = *_multi;
+    unsigned tableSize = multi.tableSize();
+    WavetableMulti::DualTable dt = multi.getInterpolationPairForFrequency(frequency);
+
+    for (unsigned i = 0; i < nframes; ++i) {
+        float position = phase * tableSize;
+        unsigned index = static_cast<unsigned>(position);
+        float frac = position - index;
+        output[i] =
+            (1 - dt.delta) * interpolate<M>(&dt.table1[index], frac) +
+            dt.delta * interpolate<M>(&dt.table2[index], frac);
+
+        phase += phaseInc;
+        phase -= static_cast<int>(phase);
+    }
+
+    _phase = phase;
+}
+
+template <InterpolatorModel M>
+void WavetableOscillator::processModulatedDual(const float* frequencies, float detuneRatio, float* output, unsigned nframes)
+{
+    float phase = _phase;
+    float sampleInterval = _sampleInterval;
+
+    const WavetableMulti& multi = *_multi;
+    unsigned tableSize = multi.tableSize();
+
+    for (unsigned i = 0; i < nframes; ++i) {
+        float frequency = frequencies[i];
+        float phaseInc = frequency * (detuneRatio * sampleInterval);
+
+        WavetableMulti::DualTable dt = multi.getInterpolationPairForFrequency(frequency);
+
+        float position = phase * tableSize;
+        unsigned index = static_cast<unsigned>(position);
+        float frac = position - index;
+        output[i] =
+            (1 - dt.delta) * interpolate<M>(&dt.table1[index], frac) +
+            dt.delta * interpolate<M>(&dt.table2[index], frac);
+
+        phase += phaseInc;
+        phase -= static_cast<int>(phase);
+    }
+
+    _phase = phase;
+}
+
+void WavetableOscillator::process(float frequency, float detuneRatio, float* output, unsigned nframes)
+{
+    int quality = clamp(_quality, 0, 3);
+
+    switch (quality) {
+    case 0:
+        processSingle<kInterpolatorNearest>(frequency, detuneRatio, output, nframes);
+        break;
+    case 1:
+        processSingle<kInterpolatorLinear>(frequency, detuneRatio, output, nframes);
+        break;
+    case 2:
+        processSingle<kInterpolatorBspline3>(frequency, detuneRatio, output, nframes);
+        break;
+    case 3:
+        processDual<kInterpolatorBspline3>(frequency, detuneRatio, output, nframes);
+        break;
+    }
+}
+
+void WavetableOscillator::processModulated(const float* frequencies, float detuneRatio, float* output, unsigned nframes)
+{
+    int quality = clamp(_quality, 0, 3);
+
+    switch (quality) {
+    case 0:
+        processModulatedSingle<kInterpolatorNearest>(frequencies, detuneRatio, output, nframes);
+        break;
+    case 1:
+        processModulatedSingle<kInterpolatorLinear>(frequencies, detuneRatio, output, nframes);
+        break;
+    case 2:
+        processModulatedSingle<kInterpolatorBspline3>(frequencies, detuneRatio, output, nframes);
+        break;
+    case 3:
+        processModulatedDual<kInterpolatorBspline3>(frequencies, detuneRatio, output, nframes);
+        break;
+    }
 }
 
 //------------------------------------------------------------------------------
@@ -200,6 +292,39 @@ unsigned WavetableRange::getOctaveForFrequency(float f)
     return clamp<int>(oct, 0, countOctaves - 1);
 }
 
+static const auto octaveForFrequencyTable = []()
+{
+    static constexpr unsigned N = 1024;
+    std::array<float, N> table;
+
+    constexpr double fmin = 1 / WavetableRange::frequencyScaleFactor;
+    constexpr double fmax = (1 << (WavetableRange::countOctaves - 1)) / WavetableRange::frequencyScaleFactor;
+
+    for (unsigned i = 0; i < N; ++i) {
+        double f = fmin + (i * (1.0 / (N - 1))) * (fmax - fmin);
+        table[i] = std::log2(f * WavetableRange::frequencyScaleFactor);
+    }
+
+    return table;
+}();
+
+float WavetableRange::getFractionalOctaveForFrequency(float f)
+{
+    static constexpr unsigned N = octaveForFrequencyTable.size();
+
+    constexpr double fmin = 1 / WavetableRange::frequencyScaleFactor;
+    constexpr double fmax = (1 << (WavetableRange::countOctaves - 1)) / WavetableRange::frequencyScaleFactor;
+
+    float pos = (f - fmin) * ((N - 1) / static_cast<float>(fmax - fmin));
+    int index1 = static_cast<int>(pos);
+    float frac = pos - index1;
+    index1 = clamp<int>(index1, 0, N - 1);
+    int index2 = std::min<int>(index1 + 1, N - 1);
+
+    return (1.0f - frac) * octaveForFrequencyTable[index1] +
+        frac * octaveForFrequencyTable[index2];
+}
+
 WavetableRange WavetableRange::getRangeForOctave(int o)
 {
     WavetableRange range;
@@ -271,7 +396,7 @@ const WavetableMulti* WavetableMulti::getSilenceWavetable()
 
 void WavetableMulti::allocateStorage(unsigned tableSize)
 {
-    _multiData.resize((tableSize + _tableExtra) * numTables());
+    _multiData.resize((tableSize + 2 * _tableExtra) * numTables());
     _tableSize = tableSize;
 }
 
@@ -282,9 +407,22 @@ void WavetableMulti::fillExtra()
     constexpr unsigned numTables = WavetableMulti::numTables();
 
     for (unsigned m = 0; m < numTables; ++m) {
-        float* ptr = const_cast<float*>(getTablePointer(m));
+        float* beg = const_cast<float*>(getTablePointer(m));
-        for (unsigned i = 0; i < tableExtra; ++i)
+        float* end = beg + tableSize;
-            ptr[tableSize + i] = ptr[i % tableSize];
+        // fill right
+        float* src = beg;
+        float* dst = end;
+        for (unsigned i = 0; i < tableExtra; ++i) {
+            *dst++ = *src;
+            src = (src + 1 != end) ? (src + 1) : beg;
+        }
+        // fill left
+        src = end - 1;
+        dst = beg - 1;
+        for (unsigned i = 0; i < tableExtra; ++i) {
+            *dst-- = *src;
+            src = (src != beg) ? (src - 1) : (end - 1);
+        }
     }
 }
 

src/sfizz/Wavetables.h

@@ -8,6 +8,7 @@
 #include "Config.h"
 #include "LeakDetector.h"
 #include "Buffer.h"
+#include "MathHelpers.h"
 #include <absl/types/span.h>
 #include <absl/container/flat_hash_map.h>
 #include <memory>
@@ -18,6 +19,8 @@ class FilePool;
 
 class WavetableMulti;
 
+enum InterpolatorModel : int;
+
 /**
    An oscillator based on wavetables
  */
@@ -44,6 +47,16 @@ public:
      */
     void setPhase(float phase);
 
+    /**
+       Set the quality of this oscillator. (cf. `oscillator_quality`)
+
+       0: nearest
+       1: linear
+       2: high
+       3: dual-high
+     */
+    void setQuality(int q) { _quality = q; }
+
     /**
        Compute a cycle of the oscillator, with constant frequency.
      */
@@ -55,17 +68,23 @@ public:
     void processModulated(const float* frequencies, float detuneRatio, float* output, unsigned nframes);
 
 private:
-    /**
+    // single-table interpolation
-       Interpolate a value from a part of table, with delta in 0 to 1 excluded.
+    template <InterpolatorModel M>
-       There are `TableExtra` elements available for reading.
+    void processSingle(float frequency, float detuneRatio, float* output, unsigned nframes);
-       (cf. WavetableMulti)
+    template <InterpolatorModel M>
-     */
+    void processModulatedSingle(const float* frequencies, float detuneRatio, float* output, unsigned nframes);
-    static float interpolate(const float* x, float delta);
+
+    // dual-table interpolation
+    template <InterpolatorModel M>
+    void processDual(float frequency, float detuneRatio, float* output, unsigned nframes);
+    template <InterpolatorModel M>
+    void processModulatedDual(const float* frequencies, float detuneRatio, float* output, unsigned nframes);
 
 private:
     float _phase = 0.0f;
     float _sampleInterval = 0.0f;
     const WavetableMulti* _multi = nullptr;
+    int _quality = 1;
     LEAK_DETECTOR(WavetableOscillator);
 };
 
@@ -113,6 +132,7 @@ public:
     static constexpr float frequencyScaleFactor = 0.05;
 
     static unsigned getOctaveForFrequency(float f);
+    static float getFractionalOctaveForFrequency(float f);
     static WavetableRange getRangeForOctave(int o);
     static WavetableRange getRangeForFrequency(float f);
 
@@ -153,6 +173,31 @@ public:
         return getTable(WavetableRange::getOctaveForFrequency(freq));
     }
 
+    // adjacent tables with interpolation factor between them
+    struct DualTable {
+        const float* table1;
+        const float* table2;
+        float delta;
+    };
+
+    // get the pair of tables at the fractional multisample position (range checked)
+    DualTable getInterpolationPair(float position) const
+    {
+        DualTable dt;
+        int index = static_cast<int>(position);
+        dt.delta = position - index;
+        dt.table1 = getTablePointer(clamp<int>(index, 0, WavetableRange::countOctaves - 1));
+        dt.table2 = getTablePointer(clamp<int>(index + 1, 0, WavetableRange::countOctaves - 1));
+        return dt;
+    }
+
+    // get the pair of tables for the given playback frequency (range checked)
+    DualTable getInterpolationPairForFrequency(float freq) const
+    {
+        float position = WavetableRange::getFractionalOctaveForFrequency(freq);
+        return getInterpolationPair(position);
+    }
+
     // create a multisample according to a given harmonic profile
     // the reference sample rate is the minimum value accepted by the DSP
     // system (most defavorable wrt. aliasing)
@@ -166,7 +211,7 @@ private:
     // get a pointer to the beginning of the N-th table
     const float* getTablePointer(unsigned index) const
     {
-        return _multiData.data() + index * (_tableSize + _tableExtra);
+        return _multiData.data() + index * (_tableSize + 2 * _tableExtra) + _tableExtra;
     }
 
     // allocate the internal data for tables of the given size

tests/WavetablesT.cpp

@@ -8,6 +8,7 @@
 #include "sfizz/MathHelpers.h"
 #include "catch2/catch.hpp"
 #include <algorithm>
+#include <cmath>
 
 TEST_CASE("[Wavetables] Frequency ranges")
 {
@@ -38,3 +39,18 @@ TEST_CASE("[Wavetables] Frequency ranges")
     REQUIRE(min_oct == 0);
     REQUIRE(max_oct == sfz::WavetableRange::countOctaves - 1);
 }
+
+TEST_CASE("[Wavetables] Octave number lookup")
+{
+    for (int note = 0; note < 128; ++note) {
+        double f = midiNoteFrequency(note);
+
+        float ref = std::log2(f * sfz::WavetableRange::frequencyScaleFactor);
+        float oct = sfz::WavetableRange::getFractionalOctaveForFrequency(f);
+
+        ref = clamp<float>(ref, 0, sfz::WavetableRange::countOctaves - 1);
+        oct = clamp<float>(oct, 0, sfz::WavetableRange::countOctaves - 1);
+
+        REQUIRE(oct == Approx(ref).margin(0.03f));
+    }
+}