diff --git a/common.mk b/common.mk index d4a63262..fe7341b7 100644 --- a/common.mk +++ b/common.mk @@ -87,6 +87,7 @@ SFIZZ_SOURCES = \ src/sfizz/FlexEGDescription.cpp \ src/sfizz/FlexEnvelope.cpp \ src/sfizz/FloatEnvelopes.cpp \ + src/sfizz/Interpolators.cpp \ src/sfizz/Logger.cpp \ src/sfizz/LFO.cpp \ src/sfizz/LFODescription.cpp \ @@ -121,7 +122,8 @@ SFIZZ_SOURCES = \ src/sfizz/Voice.cpp \ src/sfizz/VoiceManager.cpp \ src/sfizz/VoiceStealing.cpp \ - src/sfizz/Wavetables.cpp + src/sfizz/Wavetables.cpp \ + src/sfizz/WindowedSinc.cpp ### Other internal diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 47c24210..26de3632 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -112,6 +112,8 @@ set(SFIZZ_HEADERS sfizz/VoiceManager.h sfizz/VoiceStealing.h sfizz/Wavetables.h + sfizz/WindowedSinc.h + sfizz/WindowedSinc.hpp sfizz.h sfizz.hpp) @@ -151,6 +153,8 @@ set(SFIZZ_SOURCES sfizz/BeatClock.cpp sfizz/Metronome.cpp sfizz/SynthMessaging.cpp + sfizz/WindowedSinc.cpp + sfizz/Interpolators.cpp sfizz/modulations/ModId.cpp sfizz/modulations/ModKey.cpp sfizz/modulations/ModKeyHash.cpp diff --git a/src/sfizz/Interpolators.cpp b/src/sfizz/Interpolators.cpp new file mode 100644 index 00000000..5fc37a25 --- /dev/null +++ b/src/sfizz/Interpolators.cpp @@ -0,0 +1,23 @@ +// SPDX-License-Identifier: BSD-2-Clause + +// This code is part of the sfizz library and is licensed under a BSD 2-clause +// license. You should have receive a LICENSE.md file along with the code. +// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz + +#include "Interpolators.h" + +namespace sfz { + +void initializeInterpolators() +{ + SincInterpolatorTraits<8>::initialize(); + SincInterpolatorTraits<12>::initialize(); + SincInterpolatorTraits<16>::initialize(); + SincInterpolatorTraits<24>::initialize(); + SincInterpolatorTraits<36>::initialize(); + SincInterpolatorTraits<48>::initialize(); + SincInterpolatorTraits<60>::initialize(); + SincInterpolatorTraits<72>::initialize(); +} + +} // namespace sfz diff --git a/src/sfizz/Interpolators.h b/src/sfizz/Interpolators.h index c4fdbbd9..db17b720 100644 --- a/src/sfizz/Interpolators.h +++ b/src/sfizz/Interpolators.h @@ -17,8 +17,36 @@ enum InterpolatorModel : int { kInterpolatorHermite3, // a B-spline 3rd order interpolator kInterpolatorBspline3, + // a windowed-sinc 8-point interpolator + kInterpolatorSinc8, + // a windowed-sinc 12-point interpolator + kInterpolatorSinc12, + // a windowed-sinc 16-point interpolator + kInterpolatorSinc16, + // a windowed-sinc 24-point interpolator + kInterpolatorSinc24, + // a windowed-sinc 36-point interpolator + kInterpolatorSinc36, + // a windowed-sinc 48-point interpolator + kInterpolatorSinc48, + // a windowed-sinc 60-point interpolator + kInterpolatorSinc60, + // a windowed-sinc 72-point interpolator + kInterpolatorSinc72, }; +/** + * @brief Initialize interpolators + * + * This precomputes windowed-sinc tables globally. + * It needs to be called at least once, before using the windowed-sinc models. + * + * These are not computed at static initialization time, to prevent slowing down + * an audio plugin library scan (eg. VST). The static-local-variable method is + * avoided also, because we don't want this overhead on a frame-by-frame basis. + */ +void initializeInterpolators(); + /** * @brief Interpolate from a vector of values * diff --git a/src/sfizz/Interpolators.hpp b/src/sfizz/Interpolators.hpp index 5c302900..9cef3d7d 100644 --- a/src/sfizz/Interpolators.hpp +++ b/src/sfizz/Interpolators.hpp @@ -5,6 +5,7 @@ // If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz #include "Interpolators.h" +#include "WindowedSinc.h" #include "MathHelpers.h" #include "SIMDConfig.h" @@ -133,4 +134,135 @@ public: } }; +//------------------------------------------------------------------------------ +// Windowed sinc + +namespace SincInterpolatorDetail { + // See sfizz wiki page "Resampling". + constexpr size_t PointsMin = 8; + constexpr size_t PointsMax = 72; + + // Adjust Kaiser window Beta as necessary. + constexpr double BetaMin = 6.0; + constexpr double BetaMax = 10.0; + + constexpr double getBetaForNumPoints(size_t points) + { + return BetaMin + (BetaMax - BetaMin) * + (double(points - PointsMin) / double(PointsMax - PointsMin)); + } + + constexpr size_t getTableSizeForNumPoints(size_t /*points*/) + { + return 1u << 16; + } +} + +/// +template +struct SincInterpolatorTraits { + static_assert(Points == 8 || Points == 12 || Points == 16 || + Points == 24 || Points == 36 || Points == 48 || + Points == 60 || Points == 72, + "Windowed sinc size is not acceptable"); + + enum { + TableSize = SincInterpolatorDetail::getTableSizeForNumPoints(Points) + }; + + static void initialize() + { + static const FixedWindowedSinc globalInstance( + SincInterpolatorDetail::getBetaForNumPoints(Points)); + windowedSinc = &globalInstance; + } + + static const FixedWindowedSinc* windowedSinc; +}; + +template +const FixedWindowedSinc::TableSize>* +SincInterpolatorTraits::windowedSinc = nullptr; + +/// +template +class SincInterpolator; + +//------------------------------------------------------------------------------ +// Windowed sinc any order, SSE specialization +#if SFIZZ_HAVE_SSE +template +class SincInterpolator +{ +public: + static_assert(Points % 4 == 0, "Windowed sinc must be multiple of 4"); + + static inline float process(const float* values, float coeff) + { + const auto &ws = *SincInterpolatorTraits::windowedSinc; + + int j0 = 1 - int(Points) / 2; + + __m128 h[Points / 4]; + for (int i = 0; i < int(Points); ++i) + reinterpret_cast(h)[i] = ws.getUnchecked(j0 - coeff + i); + + __m128 y = _mm_mul_ps(h[0], _mm_loadu_ps(&values[j0])); + for (int i = 1; i < int(Points / 4); ++i) + y = _mm_add_ps(y, _mm_mul_ps(h[i], _mm_loadu_ps(&values[j0 + 4 * i]))); + + // sum 4 to 1 + __m128 xmm0 = y; + __m128 xmm1 = _mm_shuffle_ps(xmm0, xmm0, 0xe5); + __m128 xmm2 = _mm_movehl_ps(xmm0, xmm0); + xmm1 = _mm_add_ss(xmm1, xmm0); + xmm0 = _mm_shuffle_ps(xmm0, xmm0, 0xe7); + xmm2 = _mm_add_ss(xmm2, xmm1); + xmm0 = _mm_add_ss(xmm0, xmm2); + return _mm_cvtss_f32(xmm0); + } +}; +#endif + +//------------------------------------------------------------------------------ +// Windowed sinc any order, generic +template +class SincInterpolator +{ +public: + static inline R process(const R* values, R coeff) + { + const auto &ws = *SincInterpolatorTraits::windowedSinc; + + int j0 = 1 - int(Points) / 2; + + R h[Points]; + for (int i = 0; i < int(Points); ++i) + h[i] = R(ws.getUnchecked(j0 - coeff + i)); + + R y = h[0] * values[j0]; + for (int i = 1; i < int(Points); ++i) + y += h[i] * values[j0 + i]; + + return y; + } +}; + +template +class Interpolator : public SincInterpolator {}; +template +class Interpolator : public SincInterpolator {}; +template +class Interpolator : public SincInterpolator {}; +template +class Interpolator : public SincInterpolator {}; +template +class Interpolator : public SincInterpolator {}; +template +class Interpolator : public SincInterpolator {}; +template +class Interpolator : public SincInterpolator {}; +template +class Interpolator : public SincInterpolator {}; + } // namespace sfz diff --git a/src/sfizz/Synth.cpp b/src/sfizz/Synth.cpp index ef54fe98..ed7bf3b0 100644 --- a/src/sfizz/Synth.cpp +++ b/src/sfizz/Synth.cpp @@ -21,6 +21,7 @@ #include "utility/SpinMutex.h" #include "utility/XmlHelpers.h" #include "Voice.h" +#include "Interpolators.h" #include #include #include @@ -48,6 +49,7 @@ Synth::~Synth() Synth::Impl::Impl() { initializeSIMDDispatchers(); + initializeInterpolators(); const std::lock_guard disableCallback { callbackGuard_ }; parser_.setListener(this); diff --git a/src/sfizz/Voice.cpp b/src/sfizz/Voice.cpp index a759dee5..a93c925a 100644 --- a/src/sfizz/Voice.cpp +++ b/src/sfizz/Voice.cpp @@ -1149,18 +1149,20 @@ void Voice::Impl::fillInterpolatedWithQuality( absl::Span indices, absl::Span coeffs, absl::Span addingGains, int quality) { - switch (quality) { - default: - if (quality > 2) - goto high; // TODO sinc, not implemented - // fall through + switch (clamp(quality, 0, 10)) { + case 0: + { + constexpr auto itp = kInterpolatorNearest; + fillInterpolated(source, dest, indices, coeffs, addingGains); + } + break; case 1: { constexpr auto itp = kInterpolatorLinear; fillInterpolated(source, dest, indices, coeffs, addingGains); } break; - case 2: high: + case 2: { #if 0 // B-spline response has faster decay of aliasing, but not zero-crossings at integer positions @@ -1172,6 +1174,54 @@ void Voice::Impl::fillInterpolatedWithQuality( fillInterpolated(source, dest, indices, coeffs, addingGains); } break; + case 3: + { + constexpr auto itp = kInterpolatorSinc8; + fillInterpolated(source, dest, indices, coeffs, addingGains); + } + break; + case 4: + { + constexpr auto itp = kInterpolatorSinc12; + fillInterpolated(source, dest, indices, coeffs, addingGains); + } + break; + case 5: + { + constexpr auto itp = kInterpolatorSinc16; + fillInterpolated(source, dest, indices, coeffs, addingGains); + } + break; + case 6: + { + constexpr auto itp = kInterpolatorSinc24; + fillInterpolated(source, dest, indices, coeffs, addingGains); + } + break; + case 7: + { + constexpr auto itp = kInterpolatorSinc36; + fillInterpolated(source, dest, indices, coeffs, addingGains); + } + break; + case 8: + { + constexpr auto itp = kInterpolatorSinc48; + fillInterpolated(source, dest, indices, coeffs, addingGains); + } + break; + case 9: + { + constexpr auto itp = kInterpolatorSinc60; + fillInterpolated(source, dest, indices, coeffs, addingGains); + } + break; + case 10: + { + constexpr auto itp = kInterpolatorSinc72; + fillInterpolated(source, dest, indices, coeffs, addingGains); + } + break; } } diff --git a/src/sfizz/WindowedSinc.cpp b/src/sfizz/WindowedSinc.cpp new file mode 100644 index 00000000..94aec536 --- /dev/null +++ b/src/sfizz/WindowedSinc.cpp @@ -0,0 +1,39 @@ +// SPDX-License-Identifier: BSD-2-Clause + +// This code is part of the sfizz library and is licensed under a BSD 2-clause +// license. You should have receive a LICENSE.md file along with the code. +// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz + +#include "WindowedSinc.h" +#include "MathHelpers.h" +#include + +namespace sfz { + +void WindowedSincDetail::calculateTable(absl::Span table, size_t sincExtent, double beta, size_t extra) +{ + size_t tableSize = table.size(); + + auto window = absl::make_unique(tableSize); + kaiserWindow(beta, absl::MakeSpan(window.get(), tableSize)); + + // table domain [-N/2:+N/2] + double scale = sincExtent / static_cast(tableSize - 1); + double offset = sincExtent / -2.0; + + for (size_t i = 0; i < tableSize; ++i) { + double x = i * scale + offset; + table[i] = window[i] * normalizedSinc(x); + } + + for (size_t i = 0; i < extra; ++i) + table[extra + i] = table[tableSize - 1]; +} + +double WindowedSincDetail::calculateExact(double x, size_t sincExtent, double beta) +{ + return normalizedSinc(x) * + kaiserWindowSinglePoint(beta, (x + sincExtent / 2.0f) / sincExtent); +} + +} // namespace sfz diff --git a/src/sfizz/WindowedSinc.h b/src/sfizz/WindowedSinc.h new file mode 100644 index 00000000..94a8a9de --- /dev/null +++ b/src/sfizz/WindowedSinc.h @@ -0,0 +1,111 @@ +// SPDX-License-Identifier: BSD-2-Clause + +// This code is part of the sfizz library and is licensed under a BSD 2-clause +// license. You should have receive a LICENSE.md file along with the code. +// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz + +#pragma once +#include +#include + +namespace sfz { + +namespace WindowedSincDetail { + void calculateTable(absl::Span table, size_t sincExtent, double beta, size_t extra); + double calculateExact(double x, size_t sincExtent, double beta); +}; + +template +class AbstractWindowedSinc { +protected: + explicit AbstractWindowedSinc(double beta) noexcept : beta_(beta) {} + +public: + virtual ~AbstractWindowedSinc() noexcept {} + + // interpolate f(x), where x must be in domain [-Points/2:+Points/2] + float getUnchecked(float x) const noexcept; + + // calculate exact f(x), where x must be in domain [-Points/2:+Points/2] + double getExact(double x) const noexcept; + + // get the Kaiser window Beta parameter + double getBeta() const noexcept { return beta_; } + +protected: + void fillTable() noexcept; + + // allows interpolating f(Points/2), provided for safety + enum { TableExtra = 1 }; + +private: + double beta_ {}; +}; + +/** + * @brief Windowed-sinc using fixed compile-time parameters + * This can help to save some instructions in a resampler loop. + */ +template +class FixedWindowedSinc final : + public AbstractWindowedSinc> { +public: + using Self = FixedWindowedSinc; + using Super = AbstractWindowedSinc; + + explicit FixedWindowedSinc(double beta) : Super(beta) { Super::fillTable(); } + + // the number of points where this sinc will be evaluated (zero crossings + 1) + static constexpr size_t getNumPoints() noexcept { return Points; } + + // the size of the lookup table + static constexpr size_t getTableSize() noexcept { return TableSize; } + + // the lookup table + const float* getTablePointer() const noexcept { return table_; } + + // the lookup table + absl::Span getTableSpan() const noexcept { return absl::MakeConstSpan(table_, TableSize); } + +protected: + using Super::TableExtra; + +private: + float table_[TableSize + TableExtra]; +}; + +/** + * @brief Windowed-sinc using run-time parameters + */ +class WindowedSinc final : public AbstractWindowedSinc +{ +public: + using Self = WindowedSinc; + using Super = AbstractWindowedSinc; + + WindowedSinc(size_t points, size_t tableSize, double beta) + : Super(beta), points_(points), tableSize_(tableSize), + table_(new float[tableSize + TableExtra]) + { Super::fillTable(); } + + // the number of points where this sinc will be evaluated (zero crossings + 1) + size_t getNumPoints() const noexcept { return points_; } + + // the size of the lookup table + size_t getTableSize() const noexcept { return tableSize_; } + + // the lookup table + const float* getTablePointer() const noexcept { return table_.get(); } + + // the lookup table + absl::Span getTableSpan() const noexcept { return absl::MakeConstSpan(table_.get(), tableSize_); } + +private: + size_t points_ {}; + size_t tableSize_ {}; + std::unique_ptr table_; +}; + +} // namespace sfz + +#include "WindowedSinc.hpp" diff --git a/src/sfizz/WindowedSinc.hpp b/src/sfizz/WindowedSinc.hpp new file mode 100644 index 00000000..47d2c89d --- /dev/null +++ b/src/sfizz/WindowedSinc.hpp @@ -0,0 +1,45 @@ +// SPDX-License-Identifier: BSD-2-Clause + +// This code is part of the sfizz library and is licensed under a BSD 2-clause +// license. You should have receive a LICENSE.md file along with the code. +// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz + +#pragma once +#include "WindowedSinc.h" + +namespace sfz { + +template +inline void AbstractWindowedSinc::fillTable() noexcept +{ + float* table = const_cast(static_cast(this)->getTablePointer()); + size_t points = static_cast(this)->getNumPoints(); + size_t tableSize = static_cast(this)->getTableSize(); + + WindowedSincDetail::calculateTable( + absl::MakeSpan(table, tableSize), points, beta_, TableExtra); +} + +template +inline float AbstractWindowedSinc::getUnchecked(float x) const noexcept +{ + const float* table = static_cast(this)->getTablePointer(); + size_t points = static_cast(this)->getNumPoints(); + size_t tableSize = static_cast(this)->getTableSize(); + + float ix = (x + points / 2.0f) * ((tableSize - 1) / points); + int i0 = static_cast(ix); + float mu = ix - i0; + float y0 = table[i0]; + float dy = table[i0 + 1] - y0; + return y0 + mu * dy; +} + +template +inline double AbstractWindowedSinc::getExact(double x) const noexcept +{ + size_t points = static_cast(this)->getNumPoints(); + return WindowedSincDetail::calculateExact(x, points, beta_); +} + +} // namespace sfz diff --git a/tests/InterpolatorsT.cpp b/tests/InterpolatorsT.cpp index 93193a12..ed9c03d9 100644 --- a/tests/InterpolatorsT.cpp +++ b/tests/InterpolatorsT.cpp @@ -70,3 +70,60 @@ TEST_CASE("[Interpolators] Squares") == Approx(expected).margin(1e-2)); } } + +template +static std::pair windowedSincError(WS& ws, double step = 0.1, bool verbose = false) +{ + size_t points = ws.getNumPoints(); + double x1 = points / -2.0; + double x2 = points / +2.0; + double maxAbsErr = 0.0; + double meanAbsErr = 0.0; + //double meanSquareErr = 0.0; + + double x; + size_t n; + for (n = 0; (x = x1 + n * step) < x2; ++n) { + double val = ws.getUnchecked(x); + double ref = ws.getExact(x); + double absErr = std::fabs(val - ref); + maxAbsErr = std::max(maxAbsErr, absErr); + meanAbsErr += absErr; + //meanSquareErr += absErr * absErr; + } + meanAbsErr /= n; + //meanSquareErr /= n; + + if (verbose) { + std::cerr << "MaxAbsErr=" << maxAbsErr + << " MeanAbsErr=" << meanAbsErr + //<< " MeanSquareErr=" << meanSquareErr + << " with Points=" << points + << " TableSize=" << ws.getTableSize() + << "\n"; + } + + return { maxAbsErr, meanAbsErr }; +} + +TEST_CASE("[Interpolators] Windowed sinc precision") +{ + sfz::initializeInterpolators(); + + double maxAbsTolerance = 5e-2; + double meanAbsTolerance = 1e-3; + + auto Check = [=](std::pair maxAndMeanErr) { + REQUIRE(maxAndMeanErr.first < maxAbsTolerance); + REQUIRE(maxAndMeanErr.second < meanAbsTolerance); + }; + + Check(windowedSincError(*sfz::SincInterpolatorTraits<8>::windowedSinc)); + Check(windowedSincError(*sfz::SincInterpolatorTraits<12>::windowedSinc)); + Check(windowedSincError(*sfz::SincInterpolatorTraits<16>::windowedSinc)); + Check(windowedSincError(*sfz::SincInterpolatorTraits<24>::windowedSinc)); + Check(windowedSincError(*sfz::SincInterpolatorTraits<36>::windowedSinc)); + Check(windowedSincError(*sfz::SincInterpolatorTraits<48>::windowedSinc)); + Check(windowedSincError(*sfz::SincInterpolatorTraits<60>::windowedSinc)); + Check(windowedSincError(*sfz::SincInterpolatorTraits<72>::windowedSinc)); +}