Refactor DSP interpolation code

CMakeLists.txt

@@ -464,28 +464,18 @@ endif ( enable-floats )
 unset ( WITH_PROFILING CACHE )
 if ( enable-profiling )
     set ( WITH_PROFILING 1 )
-    if ( CMAKE_C_COMPILER_ID STREQUAL "Clang" )
+    if ( CMAKE_C_COMPILER_ID MATCHES "Clang" )
       set ( OPT_FLAGS "-Rpass=loop-vectorize -Rpass-analysis=loop-vectorize" )
-      find_program( CLANG_TIDY
-                    NAMES "clang-tidy"
-                    DOC "Path to clang-tidy executable" )
-
-      if ( CLANG_TIDY )
-        message ( STATUS "Found clang-tidy at ${CLANG_TIDY}" )
-        execute_process ( COMMAND ${CLANG_TIDY} "--version" )
-        set ( CMAKE_C_CLANG_TIDY ${CLANG_TIDY} )
-      endif ( CLANG_TIDY )
-    elseif ( CMAKE_C_COMPILER_ID STREQUAL "Intel" )
+    elseif ( CMAKE_C_COMPILER_ID MATCHES "Intel" )
       set ( OPT_FLAGS "-qopt-report=3" )
-    elseif ( CMAKE_C_COMPILER_ID STREQUAL "GNU" )
+    elseif ( CMAKE_C_COMPILER_ID MATCHES "GNU" )
       set ( OPT_FLAGS "-fopt-info -fopt-info-vec-missed" )
-	elseif ( CMAKE_C_COMPILER_ID STREQUAL "MSVC" )
+    elseif ( CMAKE_C_COMPILER_ID MATCHES "MSVC" )
 	  set ( OPT_FLAGS "/Qvec-report:2" )
     endif ( )
 
     set ( CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OPT_FLAGS}" )
     set ( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OPT_FLAGS}" )
-
 endif ( enable-profiling )
 
 unset ( ENABLE_TRAPONFPE CACHE )

src/rvoice/fluid_iir_filter_impl.cpp

@@ -156,7 +156,7 @@ static inline void fluid_iir_filter_calculate_coefficients(R fres,
  */
 template<bool GAIN_NORM, bool AMPLIFY, enum fluid_iir_filter_type TYPE>
 static void
-fluid_iir_filter_apply_local(fluid_iir_filter_t *iir_filter, fluid_real_t *dsp_buf, unsigned int count)
+fluid_iir_filter_apply_local(fluid_iir_filter_t *iir_filter, fluid_real_t *FLUID_RESTRICT dsp_buf, unsigned int count)
 {
     // FLUID_IIR_Q_LINEAR may switch the filter off by setting Q==0
     // Due to the linear smoothing, last_q may not exactly become zero.

src/rvoice/fluid_rvoice_dsp_sinc.hpp

@@ -0,0 +1,146 @@
+/* FluidSynth - A Software Synthesizer
+ *
+ * Copyright (C) 2003  Peter Hanappe and others.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1 of
+ * the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, see
+ * <https://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#define __STDCPP_WANT_MATH_SPEC_FUNCS__
+#include "fluid_sys.h"
+#include <array>
+#include <cmath>
+
+#define USE_KAISER_WINDOW 0
+
+#if USE_KAISER_WINDOW
+// Code borrowed from sfizz: https://github.com/sfztools/sfizz/blob/f5c6e29f23b8057867c08e88f5f6ac6738baa30b/src/sfizz/Interpolators.hpp#L146-L157
+constexpr fluid_real_t get_beta(int sinc_order)
+{
+    constexpr size_t PointsMin = 7;
+    constexpr size_t PointsMax = 72;
+
+    /* Kaiser shape parameter. beta ~5.658 targets ~60 dB stopband attenuation,
+     * a good compromise for short kernels: deeper stopbands can't be realized
+     * with few taps and only widen the transition band (dulling the passband). */
+    constexpr double BetaMin = 5.658;
+    constexpr double BetaMax = 10.0;
+
+    sinc_order = sinc_order < PointsMin ? PointsMin : sinc_order;
+
+    return BetaMin + (BetaMax - BetaMin) * (double(sinc_order - PointsMin) / double(PointsMax - PointsMin));
+}
+
+template<int SINC_ORDER>
+static inline fluid_real_t kaiser(fluid_real_t i_shifted)
+{
+    constexpr fluid_real_t beta = get_beta(SINC_ORDER);
+
+    /* Precompute 1 / I0(beta) once per call (normalizes the Kaiser window peak to 1). */
+    static const fluid_real_t inv_i0_beta = (fluid_real_t)1.0 / std::cyl_bessel_i(0.f, beta);
+
+    /* Half-width of the window in "tap" units. The kernel spans SINC_ORDER taps,
+     * so the window runs from -N/2 .. +N/2 around the reconstruction point. */
+    constexpr fluid_real_t halfD = (fluid_real_t)(SINC_ORDER / 2.0);
+
+    /* Kaiser window, centered on the reconstruction point.
+     * w(t) = I0(beta * sqrt(1 - (t/halfN)^2)) / I0(beta),  for |t| <= halfN.
+     * ratio can slightly exceed 1 for the outermost tap at some fractional
+     * phases; clamp the radicand to zero so we never take sqrt of a negative.
+     */
+    fluid_real_t ratio = i_shifted / halfD;
+    fluid_real_t rad = (fluid_real_t)1.0 - ratio * ratio;
+    fluid_real_t w = (rad > 0.0)
+                     ? std::cyl_bessel_i(0.f, beta * std::sqrt(rad)) * inv_i0_beta
+                     : (fluid_real_t)0.0;
+
+    return w;
+}
+#endif // USE_KAISER_WINDOW
+
+/**
+ * Normalized windowed-sinc interpolation kernel of order SINC_ORDER.
+ *
+ * The caller fills an array s[] of SINC_ORDER samples where:
+ *   s[j] == waveform[ dsp_phase_index + (j - half) ]
+ *   half == SINC_ORDER / 2
+ *
+ * @param s  Input sample array of length SINC_ORDER.
+ * @param x  Fractional phase in [0, 1].  This is the raw fractional part of
+ *           dsp_phase *after* the +0.5-sample advance applied by
+ *           fluid_rvoice_dsp_interpolate_sinc_local().  The true playback
+ *           position within s[] is therefore:
+ *
+ *               center = half + (x - 0.5)  =  half - 0.5 + x
+ *
+ *           When x == 0.5 the center falls exactly on s[half] (for odd orders only), so
+ *           fluid_interp_sinc_kernel returns s[half] unmodified (sinc(0) = 1,
+ *           sinc(n*pi) = 0 for all non-zero integers n).
+ */
+template<int SINC_ORDER>
+static inline fluid_real_t fluid_interp_sinc_kernel(const std::array<fluid_real_t, SINC_ORDER> &s, fluid_real_t x)
+{
+    static_assert(SINC_ORDER >= 1 && SINC_ORDER != 2, "SINC_ORDER must be at least 1 and not equal to 2");
+
+    constexpr int half = SINC_ORDER / 2;
+
+    const auto center = (SINC_ORDER % 2 == 0)
+                        ? (fluid_real_t)(half - 1) + x // == half + (x - 1)
+                        : (fluid_real_t)half - 0.5f + x; // == half + (x - 0.5f)
+
+
+    std::array<fluid_real_t, SINC_ORDER> coeffs;
+    fluid_real_t sum = 0.0f;
+
+    for(int i = 0; i < SINC_ORDER; i++)
+    {
+        fluid_real_t v;
+        const fluid_real_t i_shifted = (fluid_real_t)i - center;
+        const fluid_real_t arg = FLUID_M_PI * i_shifted;
+
+        if(std::fabs(arg) > 1e-6f)
+        {
+            v = std::sin(arg) / arg;
+#if USE_KAISER_WINDOW
+            const fluid_real_t k = kaiser<SINC_ORDER>(i_shifted);
+            v *= k;
+#else
+            /* Hanning window: */
+            // 0.5f * (1.0f + std::cos(arg * (fluid_real_t)(2.0 / SINC_ORDER))) == cos²(arg / SINC_ORDER)
+            const fluid_real_t hann = std::cos(arg / SINC_ORDER);
+            v *= hann * hann;
+#endif
+        }
+        else
+        {
+            v = 1.0f;
+        }
+
+        coeffs[i] = v;
+    }
+
+    fluid_real_t result = 0.0f;
+
+    for(int i = 0; i < SINC_ORDER; i++)
+    {
+        result += coeffs[i] * s[i];
+        sum += coeffs[i];
+    }
+
+    /* Normalize so coefficients always sum to 1.0, preventing amplitude
+     * modulation artifacts (harmonic distortion) as fractional phase varies. */
+    return result / sum;
+}

test/CMakeLists.txt

@@ -369,10 +369,10 @@ if(NOT CMAKE_GENERATOR MATCHES "Visual Studio")
     endforeach()
 
     add_custom_target(renderdspInterp
-        COMMAND ${MANUAL_TEST_FLUIDSYNTH} -R 0 -C 0 -g 0.6 -F "${DSPINTERP_RENDER_DIR}/sample interpolation test.${FEXT}" "sample interpolation test.mid" "sample interpolation test.sf2"
+        COMMAND ${MANUAL_TEST_FLUIDSYNTH} -R 0 -C 0 -g 0.6 -o audio.file.format=s24 -o audio.file.type=flac -r 48000 -F "${DSPINTERP_RENDER_DIR}/sample interpolation test.flac" "sample interpolation test.mid" "sample interpolation test.sf2"
         COMMENT "Rendering Interpolation examples"
         DEPENDS ${MANUAL_TEST_DEPENDS} create_out_dir
-        WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/manual/dsp_interp/
+        WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/manual/SoundFont-Spec-Test/sample interpolation test/"
         VERBATIM
     )
 

test/Wave_Sine_093.h

@@ -0,0 +1,27 @@
+// Auto-generated from : Wave_Sine_093.wav
+// Sample rate         : 44000 Hz
+// Bit depth           : 16 bit
+// Samples             : 71
+//
+// Loop points
+//   [0] type=forward  start=38  end=62  play_count=0
+
+#include <stdint.h>
+
+enum
+{
+    AUDIO_SAMPLES_LOOP0_START       = 38U,
+    AUDIO_SAMPLES_LOOP0_END         = 62U, // inclusive!
+    AUDIO_SAMPLES_LOOP0_TYPE        = 0U,  /* forward */
+    AUDIO_SAMPLES_LOOP0_PLAY_COUNT  = 0U   /* 0 = infinite */
+};
+
+const int16_t wave_sine_093[71] =
+{
+    12427,  19085,  24533,  28457,  30572,  30783,  29056,  25488,  20352,  13893,   6605,  -1131,
+    -8781, -15874, -21993, -26714, -29782, -30933, -30183, -27509, -23099, -17264, -10324,  -2764,
+     4998,  12428,  19077,  24546,  28443,  30581,  30783,  29046,  25504,  20336,  13902,   6606,
+    -1142,  -8764, -15893, -21976, -26728, -29772, -30939, -30182, -27505, -23109, -17249, -10343,
+    -2744,   4981,  12437,  19078,  24536,  28459,  30564,  30797,  29038,  25507,  20336,  13904,
+     6600,  -1131,  -8777, -15882, -21981, -26730, -29763, -30953, -30166, -27520, -23097,
+};