From 75c169217d54ba649a1d520bb9e6f941246afc23 Mon Sep 17 00:00:00 2001 From: Jean Pierre Cimalando Date: Sat, 30 Jan 2021 18:26:34 +0100 Subject: [PATCH 1/6] Add the kaiser window helper --- cmake/SfizzDeps.cmake | 6 ++ external/cephes/LICENSE.txt | 119 +++++++++++++++++++++ external/cephes/src/chbevl.c | 82 +++++++++++++++ external/cephes/src/i0.c | 193 +++++++++++++++++++++++++++++++++++ src/CMakeLists.txt | 2 +- src/sfizz/MathHelpers.h | 55 ++++++++++ 6 files changed, 456 insertions(+), 1 deletion(-) create mode 100644 external/cephes/LICENSE.txt create mode 100644 external/cephes/src/chbevl.c create mode 100644 external/cephes/src/i0.c diff --git a/cmake/SfizzDeps.cmake b/cmake/SfizzDeps.cmake index 1348b558..3dcb1916 100644 --- 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" diff --git a/external/cephes/LICENSE.txt b/external/cephes/LICENSE.txt new file mode 100644 index 00000000..666ddb13 --- /dev/null +++ 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: +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 ; 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 ; 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 +X-X-Sender: steve@d510 +To: Julien Cornebise +Subject: Re: Cephes: permission to wrap+distribute for Lua +In-Reply-To: <52653AD3.1010004@cornebise.com> +Message-ID: +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 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. diff --git a/external/cephes/src/chbevl.c b/external/cephes/src/chbevl.c new file mode 100644 index 00000000..53938816 --- /dev/null +++ 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) ); +} diff --git a/external/cephes/src/i0.c b/external/cephes/src/i0.c new file mode 100644 index 00000000..307398e9 --- /dev/null +++ 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 + +/* 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) ); + +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index df68a26b..47c24210 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -242,7 +242,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) diff --git a/src/sfizz/MathHelpers.h b/src/sfizz/MathHelpers.h index 50b012bd..d04793fc 100644 --- a/src/sfizz/MathHelpers.h +++ b/src/sfizz/MathHelpers.h @@ -25,6 +25,13 @@ #include #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 constexpr T max(T op1, T op2) { @@ -465,6 +472,54 @@ bool isReasonableAudio(absl::Span span) return true; } +/** + * @brief Compute the Kaiser window + * + * @param b Kaiser parameter beta + * @param window Span of real which receives the window + */ +template +void kaiserWindow(double b, absl::Span window) +{ + double i0b = i0(b); + for (size_t i = 0, n = window.size(); i < n; ++i) { + double x = i / static_cast(n - 1); + double t = x + x - 1.0; + window[i] = static_cast(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 +T sinc(T x) +{ + return (x == T(0)) ? T(1) : (std::sin(x) / x); +} + +/** + * @brief Compute the normalized cardinal sine + */ +template +T normalizedSinc(T x) +{ + return sinc(pi() * x); +} + /** * @brief Finds the minimum size of 2 spans * From bf3f41482fa7996c402623391407ea907e3977cc Mon Sep 17 00:00:00 2001 From: Jean Pierre Cimalando Date: Sun, 31 Jan 2021 14:04:30 +0100 Subject: [PATCH 2/6] Increase file padding, for longer windowed sinc kernels --- src/sfizz/Config.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sfizz/Config.h b/src/sfizz/Config.h index 349a9c12..2a597c67 100644 --- 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 }; /** From 4dbe8bd6020f0e0cb0a051c17c138b3437826123 Mon Sep 17 00:00:00 2001 From: Jean Pierre Cimalando Date: Sun, 31 Jan 2021 15:00:56 +0100 Subject: [PATCH 3/6] Add windowed sinc interpolation --- common.mk | 4 +- src/CMakeLists.txt | 4 ++ src/sfizz/Interpolators.cpp | 23 +++++++ src/sfizz/Interpolators.h | 28 ++++++++ src/sfizz/Interpolators.hpp | 132 ++++++++++++++++++++++++++++++++++++ src/sfizz/Synth.cpp | 2 + src/sfizz/Voice.cpp | 62 +++++++++++++++-- src/sfizz/WindowedSinc.cpp | 39 +++++++++++ src/sfizz/WindowedSinc.h | 111 ++++++++++++++++++++++++++++++ src/sfizz/WindowedSinc.hpp | 45 ++++++++++++ tests/InterpolatorsT.cpp | 57 ++++++++++++++++ 11 files changed, 500 insertions(+), 7 deletions(-) create mode 100644 src/sfizz/Interpolators.cpp create mode 100644 src/sfizz/WindowedSinc.cpp create mode 100644 src/sfizz/WindowedSinc.h create mode 100644 src/sfizz/WindowedSinc.hpp 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)); +} From d5bb51488d34531f75de273d44bc2393c077d6a1 Mon Sep 17 00:00:00 2001 From: Jean Pierre Cimalando Date: Sun, 31 Jan 2021 22:56:48 +0100 Subject: [PATCH 4/6] Improve the x64 code output slightly --- src/sfizz/WindowedSinc.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/sfizz/WindowedSinc.hpp b/src/sfizz/WindowedSinc.hpp index 47d2c89d..7908b5aa 100644 --- a/src/sfizz/WindowedSinc.hpp +++ b/src/sfizz/WindowedSinc.hpp @@ -6,6 +6,7 @@ #pragma once #include "WindowedSinc.h" +#include namespace sfz { @@ -28,7 +29,7 @@ inline float AbstractWindowedSinc::getUnchecked(float x) const noexcept size_t tableSize = static_cast(this)->getTableSize(); float ix = (x + points / 2.0f) * ((tableSize - 1) / points); - int i0 = static_cast(ix); + intptr_t i0 = static_cast(ix); float mu = ix - i0; float y0 = table[i0]; float dy = table[i0 + 1] - y0; From 0629a6b30d560389610c42070e6f3ffb5ac518bf Mon Sep 17 00:00:00 2001 From: Jean Pierre Cimalando Date: Sun, 31 Jan 2021 23:52:39 +0100 Subject: [PATCH 5/6] Add the windowed-sinc benchmark --- benchmarks/BM_interpolators.cpp | 57 ++++++++++++++++----------------- 1 file changed, 28 insertions(+), 29 deletions(-) diff --git a/benchmarks/BM_interpolators.cpp b/benchmarks/BM_interpolators.cpp index 77337576..dd44900e 100644 --- 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 input, absl::Span 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 span = absl::MakeSpan(output); \ + doInterpolation(input, span); \ + } \ + } \ + BENCHMARK_REGISTER_F(Interpolators, Type) \ + ->RangeMultiplier(4)->Range(1 << 4, 1 << 12); - for (auto _ : state) { - doInterpolation(input, absl::MakeSpan(output)); - } -} - -BENCHMARK_DEFINE_F(Interpolators, Hermite3)(benchmark::State& state) -{ - ScopedFTZ ftz; - - for (auto _ : state) { - doInterpolation(input, absl::MakeSpan(output)); - } -} - -BENCHMARK_DEFINE_F(Interpolators, Bspline3)(benchmark::State& state) -{ - ScopedFTZ ftz; - - for (auto _ : state) { - doInterpolation(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) From c7a75ee2fc3f06492753e05623b1308d11b3f504 Mon Sep 17 00:00:00 2001 From: Jean Pierre Cimalando Date: Mon, 1 Feb 2021 00:53:21 +0100 Subject: [PATCH 6/6] Accelerate windowed sinc SSE2 --- src/sfizz/Interpolators.hpp | 21 ++++++++++++--------- src/sfizz/WindowedSinc.h | 13 +++++++++++-- src/sfizz/WindowedSinc.hpp | 30 ++++++++++++++++++++++++++++++ 3 files changed, 53 insertions(+), 11 deletions(-) diff --git a/src/sfizz/Interpolators.hpp b/src/sfizz/Interpolators.hpp index 9cef3d7d..a75c9608 100644 --- a/src/sfizz/Interpolators.hpp +++ b/src/sfizz/Interpolators.hpp @@ -190,7 +190,7 @@ class SincInterpolator; //------------------------------------------------------------------------------ // Windowed sinc any order, SSE specialization -#if SFIZZ_HAVE_SSE +#if SFIZZ_HAVE_SSE2 template class SincInterpolator { @@ -201,15 +201,18 @@ public: { const auto &ws = *SincInterpolatorTraits::windowedSinc; - int j0 = 1 - int(Points) / 2; + constexpr int j0 = 1 - int(Points) / 2; + float x0 = j0 - coeff; - __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]))); + __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; diff --git a/src/sfizz/WindowedSinc.h b/src/sfizz/WindowedSinc.h index 94a8a9de..b8beb582 100644 --- a/src/sfizz/WindowedSinc.h +++ b/src/sfizz/WindowedSinc.h @@ -5,8 +5,12 @@ // If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz #pragma once +#include "SIMDConfig.h" #include #include +#if SFIZZ_HAVE_SSE2 +#include +#endif namespace sfz { @@ -26,6 +30,11 @@ public: // 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; @@ -35,8 +44,8 @@ public: protected: void fillTable() noexcept; - // allows interpolating f(Points/2), provided for safety - enum { TableExtra = 1 }; + // allows interpolating f(Points/2), and SSE, provided for safety + enum { TableExtra = 4 }; private: double beta_ {}; diff --git a/src/sfizz/WindowedSinc.hpp b/src/sfizz/WindowedSinc.hpp index 7908b5aa..1bd4f19c 100644 --- a/src/sfizz/WindowedSinc.hpp +++ b/src/sfizz/WindowedSinc.hpp @@ -36,6 +36,36 @@ inline float AbstractWindowedSinc::getUnchecked(float x) const noexcept return y0 + mu * dy; } +#if SFIZZ_HAVE_SSE2 +template +inline __m128 AbstractWindowedSinc::getUncheckedX4(__m128 x) const noexcept +{ + const float* table = static_cast(this)->getTablePointer(); + size_t points = static_cast(this)->getNumPoints(); + size_t tableSize = static_cast(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 inline double AbstractWindowedSinc::getExact(double x) const noexcept {