Replace the cubic resampler with a 4-point sinc/lanczos filter

6 files changed, 44 insertions, 31 deletions

diff --git a/Alc/ALc.c b/Alc/ALc.c
index 20787101..fcb6504b 100644
--- a/Alc/ALc.c
+++ b/Alc/ALc.c
@@ -1009,14 +1009,18 @@ static void alc_initconfig(void)
             DefaultResampler = PointResampler;
         else if(strcasecmp(str, "linear") == 0)
             DefaultResampler = LinearResampler;
+        else if(strcasecmp(str, "sinc4") == 0)
+            DefaultResampler = FIR4Resampler;
         else if(strcasecmp(str, "cubic") == 0)
-            DefaultResampler = CubicResampler;
+        {
+            ERR("Resampler option \"cubic\" is deprecated, using sinc4\n");
+            DefaultResampler = FIR4Resampler;
+        }
         else
         {
             char *end;
-
             n = strtol(str, &end, 0);
-            if(*end == '\0' && (n == PointResampler || n == LinearResampler || n == CubicResampler))
+            if(*end == '\0' && (n == PointResampler || n == LinearResampler || n == FIR4Resampler))
                 DefaultResampler = n;
             else
                 WARN("Invalid resampler: %s\n", str);
diff --git a/Alc/mixer.c b/Alc/mixer.c
index 93119e6a..b2bd1dce 100644
--- a/Alc/mixer.c
+++ b/Alc/mixer.c
@@ -43,7 +43,7 @@ static_assert((INT_MAX>>FRACTIONBITS)/MAX_PITCH > BUFFERSIZE,
 
 extern inline void InitiatePositionArrays(ALuint frac, ALuint increment, ALuint *frac_arr, ALuint *pos_arr, ALuint size);
 
-alignas(16) ALfloat CubicLUT[FRACTIONONE][4];
+alignas(16) ALfloat ResampleCoeffs[FRACTIONONE][4];
 
 
 static HrtfMixerFunc MixHrtfSamples = MixHrtf_C;
@@ -94,16 +94,16 @@ static inline ResamplerFunc SelectResampler(enum Resampler resampler)
                 return Resample_lerp32_SSE2;
 #endif
             return Resample_lerp32_C;
-        case CubicResampler:
+        case FIR4Resampler:
 #ifdef HAVE_SSE4_1
             if((CPUCapFlags&CPU_CAP_SSE4_1))
-                return Resample_cubic32_SSE41;
+                return Resample_fir4_32_SSE41;
 #endif
 #ifdef HAVE_SSE2
             if((CPUCapFlags&CPU_CAP_SSE2))
-                return Resample_cubic32_SSE2;
+                return Resample_fir4_32_SSE2;
 #endif
-            return Resample_cubic32_C;
+            return Resample_fir4_32_C;
         case ResamplerMax:
             /* Shouldn't happen */
             break;
@@ -113,17 +113,26 @@ static inline ResamplerFunc SelectResampler(enum Resampler resampler)
 }
 
 
+static float lanc2(float x)
+{
+    if(x == 0.0f)
+        return 1.0f;
+    if(fabsf(x) >= 2.0f)
+        return 0.0f;
+    return 2.0f*sinf(x*F_PI)*sinf(x*F_PI/2.0f) /
+           (F_PI*F_PI * x*x);
+}
+
 void aluInitMixer(void)
 {
     ALuint i;
     for(i = 0;i < FRACTIONONE;i++)
     {
         ALfloat mu = (ALfloat)i / FRACTIONONE;
-        ALfloat mu2 = mu*mu, mu3 = mu*mu*mu;
-        CubicLUT[i][0] = -0.5f*mu3 +       mu2 + -0.5f*mu;
-        CubicLUT[i][1] =  1.5f*mu3 + -2.5f*mu2            + 1.0f;
-        CubicLUT[i][2] = -1.5f*mu3 +  2.0f*mu2 +  0.5f*mu;
-        CubicLUT[i][3] =  0.5f*mu3 + -0.5f*mu2;
+        ResampleCoeffs[i][0] = lanc2(mu - -1.0f);
+        ResampleCoeffs[i][1] = lanc2(mu -  0.0f);
+        ResampleCoeffs[i][2] = lanc2(mu -  1.0f);
+        ResampleCoeffs[i][3] = lanc2(mu -  2.0f);
     }
 
     MixHrtfSamples = SelectHrtfMixer();
diff --git a/Alc/mixer_c.c b/Alc/mixer_c.c
index d06e52b4..6164eb00 100644
--- a/Alc/mixer_c.c
+++ b/Alc/mixer_c.c
@@ -12,8 +12,8 @@ static inline ALfloat point32(const ALfloat *vals, ALuint UNUSED(frac))
 { return vals[0]; }
 static inline ALfloat lerp32(const ALfloat *vals, ALuint frac)
 { return lerp(vals[0], vals[1], frac * (1.0f/FRACTIONONE)); }
-static inline ALfloat cubic32(const ALfloat *vals, ALuint frac)
-{ return cubic(vals[-1], vals[0], vals[1], vals[2], frac); }
+static inline ALfloat fir4_32(const ALfloat *vals, ALuint frac)
+{ return resample_fir4(vals[-1], vals[0], vals[1], vals[2], frac); }
 
 const ALfloat *Resample_copy32_C(const ALfloat *src, ALuint UNUSED(frac),
   ALuint UNUSED(increment), ALfloat *restrict dst, ALuint numsamples)
@@ -45,7 +45,7 @@ const ALfloat *Resample_##Sampler##_C(const ALfloat *src, ALuint frac,        \
 
 DECL_TEMPLATE(point32)
 DECL_TEMPLATE(lerp32)
-DECL_TEMPLATE(cubic32)
+DECL_TEMPLATE(fir4_32)
 
 #undef DECL_TEMPLATE
 
diff --git a/Alc/mixer_defs.h b/Alc/mixer_defs.h
index 00176ad4..e577fafc 100644
--- a/Alc/mixer_defs.h
+++ b/Alc/mixer_defs.h
@@ -15,7 +15,7 @@ struct HrtfState;
 const ALfloat *Resample_copy32_C(const ALfloat *src, ALuint frac, ALuint increment, ALfloat *restrict dst, ALuint dstlen);
 const ALfloat *Resample_point32_C(const ALfloat *src, ALuint frac, ALuint increment, ALfloat *restrict dst, ALuint dstlen);
 const ALfloat *Resample_lerp32_C(const ALfloat *src, ALuint frac, ALuint increment, ALfloat *restrict dst, ALuint dstlen);
-const ALfloat *Resample_cubic32_C(const ALfloat *src, ALuint frac, ALuint increment, ALfloat *restrict dst, ALuint dstlen);
+const ALfloat *Resample_fir4_32_C(const ALfloat *src, ALuint frac, ALuint increment, ALfloat *restrict dst, ALuint dstlen);
 
 
 /* C mixers */
@@ -54,9 +54,9 @@ const ALfloat *Resample_lerp32_SSE2(const ALfloat *src, ALuint frac, ALuint incr
 const ALfloat *Resample_lerp32_SSE41(const ALfloat *src, ALuint frac, ALuint increment,
                                      ALfloat *restrict dst, ALuint numsamples);
 
-const ALfloat *Resample_cubic32_SSE2(const ALfloat *src, ALuint frac, ALuint increment,
+const ALfloat *Resample_fir4_32_SSE2(const ALfloat *src, ALuint frac, ALuint increment,
                                      ALfloat *restrict dst, ALuint numsamples);
-const ALfloat *Resample_cubic32_SSE41(const ALfloat *src, ALuint frac, ALuint increment,
+const ALfloat *Resample_fir4_32_SSE41(const ALfloat *src, ALuint frac, ALuint increment,
                                       ALfloat *restrict dst, ALuint numsamples);
 
 /* Neon mixers */
diff --git a/Alc/mixer_sse2.c b/Alc/mixer_sse2.c
index 39109245..09a797c3 100644
--- a/Alc/mixer_sse2.c
+++ b/Alc/mixer_sse2.c
@@ -77,7 +77,7 @@ const ALfloat *Resample_lerp32_SSE2(const ALfloat *src, ALuint frac, ALuint incr
     return dst;
 }
 
-const ALfloat *Resample_cubic32_SSE2(const ALfloat *src, ALuint frac, ALuint increment,
+const ALfloat *Resample_fir4_32_SSE2(const ALfloat *src, ALuint frac, ALuint increment,
                                      ALfloat *restrict dst, ALuint numsamples)
 {
     const __m128i increment4 = _mm_set1_epi32(increment*4);
@@ -100,10 +100,10 @@ const ALfloat *Resample_cubic32_SSE2(const ALfloat *src, ALuint frac, ALuint inc
         const __m128 val1 = _mm_loadu_ps(&src[pos_.i[1]]);
         const __m128 val2 = _mm_loadu_ps(&src[pos_.i[2]]);
         const __m128 val3 = _mm_loadu_ps(&src[pos_.i[3]]);
-        __m128 k0 = _mm_load_ps(CubicLUT[frac_.i[0]]);
-        __m128 k1 = _mm_load_ps(CubicLUT[frac_.i[1]]);
-        __m128 k2 = _mm_load_ps(CubicLUT[frac_.i[2]]);
-        __m128 k3 = _mm_load_ps(CubicLUT[frac_.i[3]]);
+        __m128 k0 = _mm_load_ps(ResampleCoeffs[frac_.i[0]]);
+        __m128 k1 = _mm_load_ps(ResampleCoeffs[frac_.i[1]]);
+        __m128 k2 = _mm_load_ps(ResampleCoeffs[frac_.i[2]]);
+        __m128 k3 = _mm_load_ps(ResampleCoeffs[frac_.i[3]]);
         __m128 out;
 
         k0 = _mm_mul_ps(k0, val0);
@@ -130,7 +130,7 @@ const ALfloat *Resample_cubic32_SSE2(const ALfloat *src, ALuint frac, ALuint inc
 
     for(;i < numsamples;i++)
     {
-        dst[i] = cubic(src[pos], src[pos+1], src[pos+2], src[pos+3], frac);
+        dst[i] = resample_fir4(src[pos], src[pos+1], src[pos+2], src[pos+3], frac);
 
         frac += increment;
         pos  += frac>>FRACTIONBITS;
diff --git a/Alc/mixer_sse41.c b/Alc/mixer_sse41.c
index 1c859bc1..c85a733d 100644
--- a/Alc/mixer_sse41.c
+++ b/Alc/mixer_sse41.c
@@ -81,7 +81,7 @@ const ALfloat *Resample_lerp32_SSE41(const ALfloat *src, ALuint frac, ALuint inc
     return dst;
 }
 
-const ALfloat *Resample_cubic32_SSE41(const ALfloat *src, ALuint frac, ALuint increment,
+const ALfloat *Resample_fir4_32_SSE41(const ALfloat *src, ALuint frac, ALuint increment,
                                       ALfloat *restrict dst, ALuint numsamples)
 {
     const __m128i increment4 = _mm_set1_epi32(increment*4);
@@ -104,10 +104,10 @@ const ALfloat *Resample_cubic32_SSE41(const ALfloat *src, ALuint frac, ALuint in
         const __m128 val1 = _mm_loadu_ps(&src[pos_.i[1]]);
         const __m128 val2 = _mm_loadu_ps(&src[pos_.i[2]]);
         const __m128 val3 = _mm_loadu_ps(&src[pos_.i[3]]);
-        __m128 k0 = _mm_load_ps(CubicLUT[frac_.i[0]]);
-        __m128 k1 = _mm_load_ps(CubicLUT[frac_.i[1]]);
-        __m128 k2 = _mm_load_ps(CubicLUT[frac_.i[2]]);
-        __m128 k3 = _mm_load_ps(CubicLUT[frac_.i[3]]);
+        __m128 k0 = _mm_load_ps(ResampleCoeffs[frac_.i[0]]);
+        __m128 k1 = _mm_load_ps(ResampleCoeffs[frac_.i[1]]);
+        __m128 k2 = _mm_load_ps(ResampleCoeffs[frac_.i[2]]);
+        __m128 k3 = _mm_load_ps(ResampleCoeffs[frac_.i[3]]);
         __m128 out;
 
         k0 = _mm_mul_ps(k0, val0);
@@ -140,7 +140,7 @@ const ALfloat *Resample_cubic32_SSE41(const ALfloat *src, ALuint frac, ALuint in
 
     for(;i < numsamples;i++)
     {
-        dst[i] = cubic(src[pos], src[pos+1], src[pos+2], src[pos+3], frac);
+        dst[i] = resample_fir4(src[pos], src[pos+1], src[pos+2], src[pos+3], frac);
 
         frac += increment;
         pos  += frac>>FRACTIONBITS;