Merge pull request #609 from jpcima/sinc

Windowed sinc interpolation
2021-02-01 09:53:10 +01:00 · 2021-02-01 09:53:10 +01:00 · 9629fcad63
commit 9629fcad63
parent b0b0806d9e c7a75ee2fc
18 changed files with 1028 additions and 38 deletions
--- a/benchmarks/BM_interpolators.cpp
+++ b/benchmarks/BM_interpolators.cpp
@ -13,6 +13,11 @@

 class Interpolators : public benchmark::Fixture {
 public:
+    Interpolators()
+    {
+        sfz::initializeInterpolators();
+    }
+
    void SetUp(const ::benchmark::State& state)
    {
        std::random_device rd { };
@ -55,33 +60,27 @@ static void doInterpolation(absl::Span<const float> input, absl::Span<float> out
    }
 }

-BENCHMARK_DEFINE_F(Interpolators, Linear)(benchmark::State& state)
-{
-    ScopedFTZ ftz;
+#define ADD_INTERPOLATOR_BENCHMARK(Type)                                \
+    BENCHMARK_DEFINE_F(Interpolators, Type)(benchmark::State& state)    \
+    {                                                                   \
+        ScopedFTZ ftz;                                                  \
+        for (auto _ : state) {                                          \
+            absl::Span<float> span = absl::MakeSpan(output);            \
+            doInterpolation<sfz::kInterpolator##Type>(input, span);     \
+        }                                                               \
+    }                                                                   \
+    BENCHMARK_REGISTER_F(Interpolators, Type)                           \
+        ->RangeMultiplier(4)->Range(1 << 4, 1 << 12);

-    for (auto _ : state) {
-        doInterpolation<sfz::kInterpolatorLinear>(input, absl::MakeSpan(output));
-    }
-}
-
-BENCHMARK_DEFINE_F(Interpolators, Hermite3)(benchmark::State& state)
-{
-    ScopedFTZ ftz;
-
-    for (auto _ : state) {
-        doInterpolation<sfz::kInterpolatorHermite3>(input, absl::MakeSpan(output));
-    }
-}
-
-BENCHMARK_DEFINE_F(Interpolators, Bspline3)(benchmark::State& state)
-{
-    ScopedFTZ ftz;
-
-    for (auto _ : state) {
-        doInterpolation<sfz::kInterpolatorBspline3>(input, absl::MakeSpan(output));
-    }
-}
-
-BENCHMARK_REGISTER_F(Interpolators, Linear)->RangeMultiplier(4)->Range(1 << 4, 1 << 12);
-BENCHMARK_REGISTER_F(Interpolators, Hermite3)->RangeMultiplier(4)->Range(1 << 4, 1 << 12);
-BENCHMARK_REGISTER_F(Interpolators, Bspline3)->RangeMultiplier(4)->Range(1 << 4, 1 << 12);
+ADD_INTERPOLATOR_BENCHMARK(Nearest)
+ADD_INTERPOLATOR_BENCHMARK(Linear)
+ADD_INTERPOLATOR_BENCHMARK(Hermite3)
+ADD_INTERPOLATOR_BENCHMARK(Bspline3)
+ADD_INTERPOLATOR_BENCHMARK(Sinc8)
+ADD_INTERPOLATOR_BENCHMARK(Sinc12)
+ADD_INTERPOLATOR_BENCHMARK(Sinc16)
+ADD_INTERPOLATOR_BENCHMARK(Sinc24)
+ADD_INTERPOLATOR_BENCHMARK(Sinc36)
+ADD_INTERPOLATOR_BENCHMARK(Sinc48)
+ADD_INTERPOLATOR_BENCHMARK(Sinc60)
+ADD_INTERPOLATOR_BENCHMARK(Sinc72)
--- a/cmake/SfizzDeps.cmake
+++ b/cmake/SfizzDeps.cmake
@ -96,6 +96,12 @@ target_include_directories(sfizz_kissfft
    PUBLIC "src/external/kiss_fft"
    PUBLIC "src/external/kiss_fft/tools")

+# The cephes library
+add_library(sfizz_cephes STATIC
+    "external/cephes/src/chbevl.c"
+    "external/cephes/src/i0.c")
+add_library(sfizz::cephes ALIAS sfizz_cephes)
+
 # The cpuid library
 add_library(sfizz_cpuid STATIC
    "src/external/cpuid/src/cpuid/cpuinfo.cpp"
--- 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

--- a/external/cephes/LICENSE.txt
+++ b/external/cephes/LICENSE.txt
@ -0,0 +1,119 @@
+==== NOTE ====
+The actual cephes library, shipped with and wrapped by this package, is available on The Netlib at http://www.netlib.org/cephes/ . It does not have any license specified. However, its original authors, Stephen Moshier, has kindly granted permission for inclusion in a BSD-licensed package. See email snippet below for reference.
+
+Return-Path: <steve@moshier.net>
+X-Original-To: julien@cornebise.com
+Delivered-To: julien@cornebise.com
+Received: from atl4mhob11.myregisteredsite.com (atl4mhob11.myregisteredsite.com [209.17.115.49])
+    by cornebise.com (Postfix) with ESMTP id D47B139FC0
+    for <julien@cornebise.com>; Fri, 25 Oct 2013 16:32:40 +0200 (CEST)
+Received: from mailpod1.hostingplatform.com ([10.30.71.116])
+    by atl4mhob11.myregisteredsite.com (8.14.4/8.14.4) with ESMTP id r9PEWcwQ003543
+    for <julien@cornebise.com>; Fri, 25 Oct 2013 10:32:38 -0400
+Received: (qmail 11948 invoked by uid 0); 25 Oct 2013 12:36:20 -0000
+X-TCPREMOTEIP: 76.24.25.74
+X-Authenticated-UID: steve@moshier.net
+Received: from unknown (HELO d510.local) (steve@moshier.net@76.24.25.74)
+  by 0 with ESMTPA; 25 Oct 2013 12:36:20 -0000
+Date: Fri, 25 Oct 2013 08:36:19 -0400 (EDT)
+From: Stephen Moshier <steve@moshier.net>
+X-X-Sender: steve@d510
+To: Julien Cornebise <julien@cornebise.com>
+Subject: Re: Cephes: permission to wrap+distribute for Lua
+In-Reply-To: <52653AD3.1010004@cornebise.com>
+Message-ID: <alpine.DEB.2.02.1310250827040.17646@d510>
+References: <52653AD3.1010004@cornebise.com>
+User-Agent: Alpine 2.02 (DEB 1266 2009-07-14)
+MIME-Version: 1.0
+Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed
+
+
+Julien, thank you for writing.
+BSD license is fine, modification is OK.
+There are more build scripts available in the web site distributions than 
+there are on the Netlib.  I think there is an update to Planck's radiation
+function that I haven't sent to Netlib yet.  But Netlib is a more stable
+site, so it is better to cite that as a reference.
+
+
+On Mon, 21 Oct 2013, Julien Cornebise wrote:
+
+> -----BEGIN PGP SIGNED MESSAGE-----
+> Hash: SHA1
+>
+> Dear Mr Moshier
+>
+> I am a researcher in mathematics and machine learning in London, and
+> am writing about your awesome Cephes library, whom I found at the
+> heart of Scipy.
+>
+> It is so useful that, with your permission, I would like to wrap it
+> for Lua and Torch (a machine learning overlay to Lua, specialized in
+> neural nets, see http://www.torch.ch). I would like to distribute it
+> as a package for Torch, including your source code along the wrapping
+> code.
+> This wouldbe a public package, distributed under BSD License. I have
+> put a first draft on github:
+> https://github.com/jucor/torch-cephes
+>
+> Hence my three questions, please:
+>
+> 1/ How would you like to be acknowledged, beyond the comments that are
+> already in your code? Do you have any standard header/disclaimer that
+> I could add to the documentation?
+>
+> 2/ At the moment, your code is left untouched. However, if I ever need
+> to modify bits of the code, what are the conditions/restrictions?
+> Nothing huge -- I definitely do not want to mess with it: I was
+> planning to use the natural completion of some functions on the
+> completed real line (e.g. CDF returing 1 when called with "infinity",
+> or quantiles returning -Infinity when called with 0), either natively
+> if supported, or by setting a specific flag  via mtherr().
+>
+> 3/ I am currently using the source from Netlib. Do you recommend using
+> the source from your website instead ?
+>
+> Thank you very much for your attention,
+> and, more importantly, for the time and effort your poured into Cephes.
+>
+> Best regards,
+>
+> Julien Cornebise, Ph.D.
+> London, UK
+> http://www.cornebise.com/julien
+> -----BEGIN PGP SIGNATURE-----
+> Version: GnuPG v1.4.14 (Darwin)
+> Comment: GPGTools - http://gpgtools.org
+> Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/
+>
+> iEYEARECAAYFAlJlOtEACgkQKYR3gC0rw/gIpQCfZKu6+iDh9ghhm6QfsLXnldKN
+> BuIAn2zZHu1c/IrRAevhjM7N7xGg0LHO
+> =WeP5
+> -----END PGP SIGNATURE-----
+
+
+==== LICENSE ====
+Copyright (c) 2013, Julien Cornebise
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    * Neither the name of the organization nor the
+      names of its contributors may be used to endorse or promote products
+      derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--- a/external/cephes/src/chbevl.c
+++ b/external/cephes/src/chbevl.c
@ -0,0 +1,82 @@
+/*							chbevl.c
+ *
+ *	Evaluate Chebyshev series
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int N;
+ * double x, y, coef[N], chebevl();
+ *
+ * y = chbevl( x, coef, N );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Evaluates the series
+ *
+ *        N-1
+ *         - '
+ *  y  =   >   coef[i] T (x/2)
+ *         -            i
+ *        i=0
+ *
+ * of Chebyshev polynomials Ti at argument x/2.
+ *
+ * Coefficients are stored in reverse order, i.e. the zero
+ * order term is last in the array.  Note N is the number of
+ * coefficients, not the order.
+ *
+ * If coefficients are for the interval a to b, x must
+ * have been transformed to x -> 2(2x - b - a)/(b-a) before
+ * entering the routine.  This maps x from (a, b) to (-1, 1),
+ * over which the Chebyshev polynomials are defined.
+ *
+ * If the coefficients are for the inverted interval, in
+ * which (a, b) is mapped to (1/b, 1/a), the transformation
+ * required is x -> 2(2ab/x - b - a)/(b-a).  If b is infinity,
+ * this becomes x -> 4a/x - 1.
+ *
+ *
+ *
+ * SPEED:
+ *
+ * Taking advantage of the recurrence properties of the
+ * Chebyshev polynomials, the routine requires one more
+ * addition per loop than evaluating a nested polynomial of
+ * the same degree.
+ *
+ */
+/*							chbevl.c	*/
+
+/*
+Cephes Math Library Release 2.0:  April, 1987
+Copyright 1985, 1987 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+double chbevl( x, array, n )
+double x;
+double array[];
+int n;
+{
+double b0, b1, b2, *p;
+int i;
+
+p = array;
+b0 = *p++;
+b1 = 0.0;
+i = n - 1;
+
+do
+	{
+	b2 = b1;
+	b1 = b0;
+	b0 = x * b1  -  b2  + *p++;
+	}
+while( --i );
+
+return( 0.5*(b0-b2) );
+}
--- a/external/cephes/src/i0.c
+++ b/external/cephes/src/i0.c
@ -0,0 +1,193 @@
+/*							i0.c
+ *
+ *	Modified Bessel function of order zero
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, i0();
+ *
+ * y = i0( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns modified Bessel function of order zero of the
+ * argument.
+ *
+ * The function is defined as i0(x) = j0( ix ).
+ *
+ * The range is partitioned into the two intervals [0,8] and
+ * (8, infinity).  Chebyshev polynomial expansions are employed
+ * in each interval.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       0,30         6000       8.2e-17     1.9e-17
+ *    IEEE      0,30        30000       5.8e-16     1.4e-16
+ *
+ */
+/*							i0e.c
+ *
+ *	Modified Bessel function of order zero,
+ *	exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, i0e();
+ *
+ * y = i0e( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of order zero of the argument.
+ *
+ * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,30        30000       5.4e-16     1.2e-16
+ * See i0().
+ *
+ */
+
+/*							i0.c		*/
+
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*/
+
+#include <math.h>
+
+/* Chebyshev coefficients for exp(-x) I0(x)
+ * in the interval [0,8].
+ *
+ * lim(x->0){ exp(-x) I0(x) } = 1.
+ */
+
+static double A[] =
+{
+-4.41534164647933937950E-18,
+ 3.33079451882223809783E-17,
+-2.43127984654795469359E-16,
+ 1.71539128555513303061E-15,
+-1.16853328779934516808E-14,
+ 7.67618549860493561688E-14,
+-4.85644678311192946090E-13,
+ 2.95505266312963983461E-12,
+-1.72682629144155570723E-11,
+ 9.67580903537323691224E-11,
+-5.18979560163526290666E-10,
+ 2.65982372468238665035E-9,
+-1.30002500998624804212E-8,
+ 6.04699502254191894932E-8,
+-2.67079385394061173391E-7,
+ 1.11738753912010371815E-6,
+-4.41673835845875056359E-6,
+ 1.64484480707288970893E-5,
+-5.75419501008210370398E-5,
+ 1.88502885095841655729E-4,
+-5.76375574538582365885E-4,
+ 1.63947561694133579842E-3,
+-4.32430999505057594430E-3,
+ 1.05464603945949983183E-2,
+-2.37374148058994688156E-2,
+ 4.93052842396707084878E-2,
+-9.49010970480476444210E-2,
+ 1.71620901522208775349E-1,
+-3.04682672343198398683E-1,
+ 6.76795274409476084995E-1
+};
+
+
+/* Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
+ * in the inverted interval [8,infinity].
+ *
+ * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
+ */
+
+static double B[] =
+{
+-7.23318048787475395456E-18,
+-4.83050448594418207126E-18,
+ 4.46562142029675999901E-17,
+ 3.46122286769746109310E-17,
+-2.82762398051658348494E-16,
+-3.42548561967721913462E-16,
+ 1.77256013305652638360E-15,
+ 3.81168066935262242075E-15,
+-9.55484669882830764870E-15,
+-4.15056934728722208663E-14,
+ 1.54008621752140982691E-14,
+ 3.85277838274214270114E-13,
+ 7.18012445138366623367E-13,
+-1.79417853150680611778E-12,
+-1.32158118404477131188E-11,
+-3.14991652796324136454E-11,
+ 1.18891471078464383424E-11,
+ 4.94060238822496958910E-10,
+ 3.39623202570838634515E-9,
+ 2.26666899049817806459E-8,
+ 2.04891858946906374183E-7,
+ 2.89137052083475648297E-6,
+ 6.88975834691682398426E-5,
+ 3.36911647825569408990E-3,
+ 8.04490411014108831608E-1
+};
+
+
+extern double chbevl ( double, void *, int );
+
+double i0(x)
+double x;
+{
+double y;
+
+if( x < 0 )
+	x = -x;
+if( x <= 8.0 )
+	{
+	y = (x/2.0) - 2.0;
+	return( exp(x) * chbevl( y, A, 30 ) );
+	}
+
+return(  exp(x) * chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) );
+
+}
+
+
+
+
+double i0e( x )
+double x;
+{
+double y;
+
+if( x < 0 )
+	x = -x;
+if( x <= 8.0 )
+	{
+	y = (x/2.0) - 2.0;
+	return( chbevl( y, A, 30 ) );
+	}
+
+return(  chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) );
+
+}
--- 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
@ -242,7 +246,7 @@ target_sources(sfizz_internal PRIVATE ${SFIZZ_HEADERS} ${SFIZZ_SOURCES} ${FAUST_
 target_include_directories(sfizz_internal PUBLIC "." "sfizz")
 target_link_libraries(sfizz_internal
    PUBLIC absl::strings absl::span sfizz::filesystem sfizz::atomic_queue
-    PRIVATE sfizz::parser sfizz::messaging absl::flat_hash_map Threads::Threads st_audiofile sfizz::pugixml sfizz::spline sfizz::tunings sfizz::hiir sfizz::kissfft sfizz::cpuid sfizz::threadpool sfizz::jsl sfizz::atomic)
+    PRIVATE sfizz::parser sfizz::messaging absl::flat_hash_map Threads::Threads st_audiofile sfizz::pugixml sfizz::spline sfizz::tunings sfizz::hiir sfizz::kissfft sfizz::cephes sfizz::cpuid sfizz::threadpool sfizz::jsl sfizz::atomic)
 if(SFIZZ_USE_SNDFILE)
    target_compile_definitions(sfizz_internal PUBLIC "SFIZZ_USE_SNDFILE=1")
    target_link_libraries(sfizz_internal PUBLIC st_audiofile)
--- a/src/sfizz/Config.h
+++ b/src/sfizz/Config.h
@ -67,7 +67,7 @@ namespace config {
    constexpr int chunkSize { 1024 };
    constexpr unsigned int defaultAlignment { 16 };
    constexpr int filtersInPool { maxVoices * 2 };
-    constexpr int excessFileFrames { 8 };
+    constexpr int excessFileFrames { 64 };
    constexpr int maxLFOSubs { 8 };
    constexpr int maxLFOSteps { 128 };
    /**
--- a/src/sfizz/Interpolators.cpp
+++ 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
--- 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
 *
--- 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,138 @@ 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 <size_t Points>
+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<Points, TableSize> globalInstance(
+            SincInterpolatorDetail::getBetaForNumPoints(Points));
+        windowedSinc = &globalInstance;
+    }
+
+    static const FixedWindowedSinc<Points, TableSize>* windowedSinc;
+};
+
+template <size_t Points>
+const FixedWindowedSinc<Points, SincInterpolatorTraits<Points>::TableSize>*
+SincInterpolatorTraits<Points>::windowedSinc = nullptr;
+
+///
+template <class R, size_t Points>
+class SincInterpolator;
+
+//------------------------------------------------------------------------------
+// Windowed sinc any order, SSE specialization
+#if SFIZZ_HAVE_SSE2
+template <size_t Points>
+class SincInterpolator<float, Points>
+{
+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<Points>::windowedSinc;
+
+        constexpr int j0 = 1 - int(Points) / 2;
+        float x0 = j0 - coeff;
+
+        __m128 y = _mm_set1_ps(0.0f);
+        __m128 x = _mm_add_ps(_mm_set1_ps(x0), _mm_setr_ps(0, 1, 2, 3));
+        size_t i = 0;
+        do {
+            __m128 h = ws.getUncheckedX4(x);
+            y = _mm_add_ps(y, _mm_mul_ps(h, _mm_loadu_ps(&values[j0 + i])));
+            x = _mm_add_ps(x, _mm_set1_ps(4.0f));
+            i += 4;
+        } while (i < Points);
+
+        // 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 R, size_t Points>
+class SincInterpolator
+{
+public:
+    static inline R process(const R* values, R coeff)
+    {
+        const auto &ws = *SincInterpolatorTraits<Points>::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 R>
+class Interpolator<kInterpolatorSinc8, R> : public SincInterpolator<R, 8> {};
+template <class R>
+class Interpolator<kInterpolatorSinc12, R> : public SincInterpolator<R, 12> {};
+template <class R>
+class Interpolator<kInterpolatorSinc16, R> : public SincInterpolator<R, 16> {};
+template <class R>
+class Interpolator<kInterpolatorSinc24, R> : public SincInterpolator<R, 24> {};
+template <class R>
+class Interpolator<kInterpolatorSinc36, R> : public SincInterpolator<R, 36> {};
+template <class R>
+class Interpolator<kInterpolatorSinc48, R> : public SincInterpolator<R, 48> {};
+template <class R>
+class Interpolator<kInterpolatorSinc60, R> : public SincInterpolator<R, 60> {};
+template <class R>
+class Interpolator<kInterpolatorSinc72, R> : public SincInterpolator<R, 72> {};
+
 } // namespace sfz
--- a/src/sfizz/MathHelpers.h
+++ b/src/sfizz/MathHelpers.h
@ -25,6 +25,13 @@
 #include <xmmintrin.h>
 #endif

+#if __cplusplus >= 201703L
+static double i0(double x) { return std::cyl_bessel_i(0.0, x); }
+#else
+// external Bessel function from cephes
+extern "C" double i0(double x);
+#endif
+
 template <class T>
 constexpr T max(T op1, T op2)
 {
@ -465,6 +472,54 @@ bool isReasonableAudio(absl::Span<Type> span)
    return true;
 }

+/**
+ * @brief Compute the Kaiser window
+ *
+ * @param b Kaiser parameter beta
+ * @param window Span of real which receives the window
+ */
+template <class T>
+void kaiserWindow(double b, absl::Span<T> window)
+{
+    double i0b = i0(b);
+    for (size_t i = 0, n = window.size(); i < n; ++i) {
+        double x = i / static_cast<double>(n - 1);
+        double t = x + x - 1.0;
+        window[i] = static_cast<T>(i0(b * std::sqrt(1.0 - t * t)) / i0b);
+    }
+}
+
+/**
+ * @brief Compute a single point of the Kaiser window
+ * This is less efficient than calculating the whole window at once.
+ *
+ * @param b Kaiser parameter beta
+ * @param x Point to evaluate, normalized in 0 to 1
+ */
+inline double kaiserWindowSinglePoint(double b, double x)
+{
+    double t = x + x - 1.0;
+    return i0(b * std::sqrt(1.0 - t * t)) / i0(b);
+}
+
+/**
+ * @brief Compute the cardinal sine
+ */
+template <class T>
+T sinc(T x)
+{
+    return (x == T(0)) ? T(1) : (std::sin(x) / x);
+}
+
+/**
+ * @brief Compute the normalized cardinal sine
+ */
+template <class T>
+T normalizedSinc(T x)
+{
+    return sinc(pi<T>() * x);
+}
+
 /**
 * @brief Finds the minimum size of 2 spans
 *
--- 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 <absl/algorithm/container.h>
 #include <absl/memory/memory.h>
 #include <absl/strings/str_replace.h>
@ -48,6 +49,7 @@ Synth::~Synth()
 Synth::Impl::Impl()
 {
    initializeSIMDDispatchers();
+    initializeInterpolators();

    const std::lock_guard<SpinMutex> disableCallback { callbackGuard_ };
    parser_.setListener(this);
--- a/src/sfizz/Voice.cpp
+++ b/src/sfizz/Voice.cpp
@ -1149,18 +1149,20 @@ void Voice::Impl::fillInterpolatedWithQuality(
    absl::Span<const int> indices, absl::Span<const float> coeffs,
    absl::Span<const float> 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<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
    case 1:
        {
            constexpr auto itp = kInterpolatorLinear;
            fillInterpolated<itp, Adding>(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<itp, Adding>(source, dest, indices, coeffs, addingGains);
        }
        break;
+    case 3:
+        {
+            constexpr auto itp = kInterpolatorSinc8;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 4:
+        {
+            constexpr auto itp = kInterpolatorSinc12;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 5:
+        {
+            constexpr auto itp = kInterpolatorSinc16;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 6:
+        {
+            constexpr auto itp = kInterpolatorSinc24;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 7:
+        {
+            constexpr auto itp = kInterpolatorSinc36;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 8:
+        {
+            constexpr auto itp = kInterpolatorSinc48;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 9:
+        {
+            constexpr auto itp = kInterpolatorSinc60;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 10:
+        {
+            constexpr auto itp = kInterpolatorSinc72;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
    }
 }

--- a/src/sfizz/WindowedSinc.cpp
+++ 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 <absl/memory/memory.h>
+
+namespace sfz {
+
+void WindowedSincDetail::calculateTable(absl::Span<float> table, size_t sincExtent, double beta, size_t extra)
+{
+    size_t tableSize = table.size();
+
+    auto window = absl::make_unique<float[]>(tableSize);
+    kaiserWindow(beta, absl::MakeSpan(window.get(), tableSize));
+
+    // table domain [-N/2:+N/2]
+    double scale = sincExtent / static_cast<double>(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
--- a/src/sfizz/WindowedSinc.h
+++ b/src/sfizz/WindowedSinc.h
@ -0,0 +1,120 @@
+// 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 "SIMDConfig.h"
+#include <absl/types/span.h>
+#include <memory>
+#if SFIZZ_HAVE_SSE2
+#include <xmmintrin.h>
+#endif
+
+namespace sfz {
+
+namespace WindowedSincDetail {
+    void calculateTable(absl::Span<float> table, size_t sincExtent, double beta, size_t extra);
+    double calculateExact(double x, size_t sincExtent, double beta);
+};
+
+template <class T>
+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;
+
+#if SFIZZ_HAVE_SSE2
+    // interpolate f(x), 4 values at once
+    __m128 getUncheckedX4(__m128 x) const noexcept;
+#endif
+
+    // 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), and SSE, provided for safety
+    enum { TableExtra = 4 };
+
+private:
+    double beta_ {};
+};
+
+/**
+ * @brief Windowed-sinc using fixed compile-time parameters
+ * This can help to save some instructions in a resampler loop.
+ */
+template <size_t Points, size_t TableSize>
+class FixedWindowedSinc final :
+        public AbstractWindowedSinc<FixedWindowedSinc<Points, TableSize>> {
+public:
+    using Self = FixedWindowedSinc<Points, TableSize>;
+    using Super = AbstractWindowedSinc<Self>;
+
+    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<const float> 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<WindowedSinc>
+{
+public:
+    using Self = WindowedSinc;
+    using Super = AbstractWindowedSinc<WindowedSinc>;
+
+    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<const float> getTableSpan() const noexcept { return absl::MakeConstSpan(table_.get(), tableSize_); }
+
+private:
+    size_t points_ {};
+    size_t tableSize_ {};
+    std::unique_ptr<float[]> table_;
+};
+
+} // namespace sfz
+
+#include "WindowedSinc.hpp"
--- a/src/sfizz/WindowedSinc.hpp
+++ b/src/sfizz/WindowedSinc.hpp
@ -0,0 +1,76 @@
+// 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"
+#include <cstdint>
+
+namespace sfz {
+
+template <class T>
+inline void AbstractWindowedSinc<T>::fillTable() noexcept
+{
+    float* table = const_cast<float*>(static_cast<T*>(this)->getTablePointer());
+    size_t points = static_cast<T*>(this)->getNumPoints();
+    size_t tableSize = static_cast<T*>(this)->getTableSize();
+
+    WindowedSincDetail::calculateTable(
+        absl::MakeSpan(table, tableSize), points, beta_, TableExtra);
+}
+
+template <class T>
+inline float AbstractWindowedSinc<T>::getUnchecked(float x) const noexcept
+{
+    const float* table = static_cast<const T*>(this)->getTablePointer();
+    size_t points = static_cast<const T*>(this)->getNumPoints();
+    size_t tableSize = static_cast<const T*>(this)->getTableSize();
+
+    float ix = (x + points / 2.0f) * ((tableSize - 1) / points);
+    intptr_t i0 = static_cast<intptr_t>(ix);
+    float mu = ix - i0;
+    float y0 = table[i0];
+    float dy = table[i0 + 1] - y0;
+    return y0 + mu * dy;
+}
+
+#if SFIZZ_HAVE_SSE2
+template <class T>
+inline __m128 AbstractWindowedSinc<T>::getUncheckedX4(__m128 x) const noexcept
+{
+    const float* table = static_cast<const T*>(this)->getTablePointer();
+    size_t points = static_cast<const T*>(this)->getNumPoints();
+    size_t tableSize = static_cast<const T*>(this)->getTableSize();
+
+    __m128 ix = _mm_mul_ps(
+        _mm_add_ps(x, _mm_set1_ps(points / 2.0f)),
+        _mm_set1_ps((tableSize - 1) / points));
+    alignas(__m128i) int j0[4];
+    __m128i i0 = _mm_cvttps_epi32(ix);
+    _mm_store_si128((__m128i*)j0, i0);
+    __m128 mu = _mm_sub_ps(ix, _mm_cvtepi32_ps(i0));
+
+    // reference: Interpolated table lookups using SSE2 [2/2]
+    // https://rawstudio.org/blog/?p=482
+    __m128 p0p1 = _mm_castsi128_ps(_mm_loadl_epi64((__m128i*)&table[j0[0]]));
+    __m128 p2p3 = _mm_castsi128_ps(_mm_loadl_epi64((__m128i*)&table[j0[2]]));
+    p0p1 = _mm_loadh_pi(p0p1, (__m64*)&table[j0[1]]);
+    p2p3 = _mm_loadh_pi(p2p3, (__m64*)&table[j0[3]]);
+    __m128 y0 = _mm_shuffle_ps(p0p1, p2p3, _MM_SHUFFLE(2, 0, 2, 0));
+    __m128 y1 = _mm_shuffle_ps(p0p1, p2p3, _MM_SHUFFLE(3, 1, 3, 1));
+
+    __m128 dy = _mm_sub_ps(y1, y0);
+    return _mm_add_ps(y0, _mm_mul_ps(mu, dy));
+}
+#endif
+
+template <class T>
+inline double AbstractWindowedSinc<T>::getExact(double x) const noexcept
+{
+    size_t points = static_cast<const T*>(this)->getNumPoints();
+    return WindowedSincDetail::calculateExact(x, points, beta_);
+}
+
+} // namespace sfz
--- a/tests/InterpolatorsT.cpp
+++ b/tests/InterpolatorsT.cpp
@ -70,3 +70,60 @@ TEST_CASE("[Interpolators] Squares")
            == Approx(expected).margin(1e-2));
    }
 }
+
+template <class WS>
+static std::pair<double, double> 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<double, double> 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));
+}