Replace the Lanczos window with Kaiser for the sinc resampler

4 files changed, 93 insertions, 23 deletions

diff --git a/Alc/mixer.c b/Alc/mixer.c
index 4917477a..e7a924bf 100644
--- a/Alc/mixer.c
+++ b/Alc/mixer.c
@@ -133,17 +133,87 @@ static inline ResamplerFunc SelectResampler(enum Resampler resampler)
     return Resample_point32_C;
 }
 
+
+/* The sinc resampler makes use of a Kaiser window to limit the needed sample
+ * points to 4 and 8, respectively.
+ */
+
 #ifndef M_PI
 #define M_PI                         (3.14159265358979323846)
 #endif
-static float lanc(double r, double x)
+static inline double Sinc(double x)
 {
-    if(x == 0.0) return 1.0f;
-    if(fabs(x) >= r) return 0.0f;
-    return (float)(r*sin(x*M_PI)*sin(x*M_PI/r) /
-                   (M_PI*M_PI * x*x));
+    if(x == 0.0) return 1.0;
+    return sin(x*M_PI) / (x*M_PI);
 }
 
+/* The zero-order modified Bessel function of the first kind, used for the
+ * Kaiser window.
+ *
+ *   I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
+ *          = sum_{k=0}^inf ((x / 2)^k / k!)^2
+ */
+static double BesselI_0(double x)
+{
+    double term, sum, x2, y, last_sum;
+    int k;
+
+    /* Start at k=1 since k=0 is trivial. */
+    term = 1.0;
+    sum = 1.0;
+    x2 = x / 2.0;
+    k = 1;
+
+    /* Let the integration converge until the term of the sum is no longer
+     * significant.
+     */
+    do {
+        y = x2 / k;
+        k ++;
+        last_sum = sum;
+        term *= y * y;
+        sum += term;
+    } while(sum != last_sum);
+    return sum;
+}
+
+/* Calculate a Kaiser window from the given beta value and a normalized k
+ * [-1, 1].
+ *
+ *   w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B),  -1 <= k <= 1
+ *          { 0,                              elsewhere.
+ *
+ * Where k can be calculated as:
+ *
+ *   k = i / l,         where -l <= i <= l.
+ *
+ * or:
+ *
+ *   k = 2 i / M - 1,   where 0 <= i <= M.
+ */
+static inline double Kaiser(double b, double k)
+{
+    if(k <= -1.0 || k >= 1.0) return 0.0;
+    return BesselI_0(b * sqrt(1.0 - (k*k))) / BesselI_0(b);
+}
+
+static inline double CalcKaiserBeta(double rejection)
+{
+    if(rejection > 50.0)
+        return 0.1102 * (rejection - 8.7);
+    if(rejection >= 21.0)
+        return (0.5842 * pow(rejection - 21.0, 0.4)) +
+               (0.07886 * (rejection - 21.0));
+    return 0.0;
+}
+
+static float SincKaiser(double r, double x)
+{
+    /* Limit rippling to -90dB. */
+    return Kaiser(CalcKaiserBeta(90.0), x / r) * Sinc(x);
+}
+
+
 void aluInitMixer(void)
 {
     enum Resampler resampler = ResamplerDefault;
@@ -180,23 +250,23 @@ void aluInitMixer(void)
         for(i = 0;i < FRACTIONONE;i++)
         {
             ALdouble mu = (ALdouble)i / FRACTIONONE;
-            ResampleCoeffs.FIR8[i][0] = lanc(4.0, mu - -3.0);
-            ResampleCoeffs.FIR8[i][1] = lanc(4.0, mu - -2.0);
-            ResampleCoeffs.FIR8[i][2] = lanc(4.0, mu - -1.0);
-            ResampleCoeffs.FIR8[i][3] = lanc(4.0, mu -  0.0);
-            ResampleCoeffs.FIR8[i][4] = lanc(4.0, mu -  1.0);
-            ResampleCoeffs.FIR8[i][5] = lanc(4.0, mu -  2.0);
-            ResampleCoeffs.FIR8[i][6] = lanc(4.0, mu -  3.0);
-            ResampleCoeffs.FIR8[i][7] = lanc(4.0, mu -  4.0);
+            ResampleCoeffs.FIR8[i][0] = SincKaiser(4.0, mu - -3.0);
+            ResampleCoeffs.FIR8[i][1] = SincKaiser(4.0, mu - -2.0);
+            ResampleCoeffs.FIR8[i][2] = SincKaiser(4.0, mu - -1.0);
+            ResampleCoeffs.FIR8[i][3] = SincKaiser(4.0, mu -  0.0);
+            ResampleCoeffs.FIR8[i][4] = SincKaiser(4.0, mu -  1.0);
+            ResampleCoeffs.FIR8[i][5] = SincKaiser(4.0, mu -  2.0);
+            ResampleCoeffs.FIR8[i][6] = SincKaiser(4.0, mu -  3.0);
+            ResampleCoeffs.FIR8[i][7] = SincKaiser(4.0, mu -  4.0);
         }
     else if(resampler == FIR4Resampler)
         for(i = 0;i < FRACTIONONE;i++)
         {
             ALdouble mu = (ALdouble)i / FRACTIONONE;
-            ResampleCoeffs.FIR4[i][0] = lanc(2.0, mu - -1.0);
-            ResampleCoeffs.FIR4[i][1] = lanc(2.0, mu -  0.0);
-            ResampleCoeffs.FIR4[i][2] = lanc(2.0, mu -  1.0);
-            ResampleCoeffs.FIR4[i][3] = lanc(2.0, mu -  2.0);
+            ResampleCoeffs.FIR4[i][0] = SincKaiser(2.0, mu - -1.0);
+            ResampleCoeffs.FIR4[i][1] = SincKaiser(2.0, mu -  0.0);
+            ResampleCoeffs.FIR4[i][2] = SincKaiser(2.0, mu -  1.0);
+            ResampleCoeffs.FIR4[i][3] = SincKaiser(2.0, mu -  2.0);
         }
 
     MixHrtfSamples = SelectHrtfMixer();
diff --git a/ChangeLog b/ChangeLog
index 40ff34ed..e0fdd4cd 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -6,7 +6,7 @@ openal-soft-1.17.0:
 
     Implemented the ALC_SOFT_HRTF extension.
 
-    Implemented C, SSE3, and SSE4.1 based 4- and 8-point Sinc-Lanczos resamplers.
+    Implemented C, SSE3, and SSE4.1 based 4- and 8-point Sinc resamplers.
 
     Implemented B-Format output support for the wave file writer. This creates
     FuMa-style first-order Ambisonics wave files (AMB format).
diff --git a/alsoftrc.sample b/alsoftrc.sample
index 84e57f4c..70bfbcf2 100644
--- a/alsoftrc.sample
+++ b/alsoftrc.sample
@@ -129,8 +129,8 @@
 #  Selects the resampler used when mixing sources. Valid values are:
 #  point - nearest sample, no interpolation
 #  linear - extrapolates samples using a linear slope between samples
-#  sinc4 - extrapolates samples using a 4-point Sinc-Lanczos filter
-#  sinc8 - extrapolates samples using an 8-point Sinc-Lanczos filter
+#  sinc4 - extrapolates samples using a 4-point Sinc filter
+#  sinc8 - extrapolates samples using an 8-point Sinc filter
 #  Specifying other values will result in using the default (linear).
 #resampler = linear
 
diff --git a/utils/alsoft-config/mainwindow.cpp b/utils/alsoft-config/mainwindow.cpp
index dff8b930..75703f03 100644
--- a/utils/alsoft-config/mainwindow.cpp
+++ b/utils/alsoft-config/mainwindow.cpp
@@ -89,8 +89,8 @@ static const struct {
     { "Default", "" },
     { "Point (low quality, fast)", "point" },
     { "Linear (basic quality, fast)", "linear" },
-    { "Sinc-Lanczos 4pt (good quality)", "sinc4" },
-    { "Sinc-Lanczos 8pt (high quality, slow)", "sinc8" },
+    { "4-Point Sinc (good quality)", "sinc4" },
+    { "8-Point Sinc (high quality, slow)", "sinc8" },
 
     { "", "" }
 }, stereoModeList[] = {
@@ -355,7 +355,7 @@ void MainWindow::loadConfig(const QString &fname)
     if(resampler.isEmpty() == false)
     {
         /* The "cubic" resampler is no longer supported. It's been replaced by
-         * "sinc4" (4-point Sinc-Lanczos). */
+         * "sinc4". */
         if(resampler == "cubic")
             resampler = "sinc4";