Uses integer arithmetics in SDL_ResampleAudio

- Avoids precision loss caused by large floating point numbers.
- Adds unit test to test the signal-to-noise ratio and maximum error of resampler.
- Code cleanup

src/audio/SDL_audiocvt.c

@@ -177,14 +177,18 @@ SDL_ConvertMonoToStereo_SSE(SDL_AudioCVT * cvt, SDL_AudioFormat format)
 
 #include "SDL_audio_resampler_filter.h"
 
-static int
+static Sint32
-ResamplerPadding(const int inrate, const int outrate)
+ResamplerPadding(const Sint32 inrate, const Sint32 outrate)
 {
+    /* This function uses integer arithmetics to avoid precision loss caused
+     * by large floating point numbers. Sint32 is needed for the large number
+     * multiplication. The integers are assumed to be non-negative so that
+     * division rounds by truncation. */
     if (inrate == outrate) {
         return 0;
     }
     if (inrate > outrate) {
-        return (int) SDL_ceilf(((float) (RESAMPLER_SAMPLES_PER_ZERO_CROSSING * inrate) / ((float) outrate)));
+        return (RESAMPLER_SAMPLES_PER_ZERO_CROSSING * inrate + outrate - 1) / outrate;
     }
     return RESAMPLER_SAMPLES_PER_ZERO_CROSSING;
 }
@@ -196,57 +200,59 @@ SDL_ResampleAudio(const int chans, const int inrate, const int outrate,
                         const float *inbuf, const int inbuflen,
                         float *outbuf, const int outbuflen)
 {
-    /* Note that this used to be double, but it looks like we can get by with float in most cases at
+    /* This function uses integer arithmetics to avoid precision loss caused
-       almost twice the speed on Intel processors, and orders of magnitude more
+     * by large floating point numbers. For some operations, Sint32 or Sint64
-       on CPUs that need a software fallback for double calculations. */
+     * are needed for the large number multiplications. The input integers are
-    typedef float ResampleFloatType;
+     * assumed to be non-negative so that division rounds by truncation and
-
+     * modulo is always non-negative. Note that the operator order is important
-    const ResampleFloatType finrate = (ResampleFloatType) inrate;
+     * for these integer divisions. */
-    const ResampleFloatType ratio = ((float) outrate) / ((float) inrate);
     const int paddinglen = ResamplerPadding(inrate, outrate);
     const int framelen = chans * (int)sizeof (float);
     const int inframes = inbuflen / framelen;
-    const int wantedoutframes = (int) ((inbuflen / framelen) * ratio);  /* outbuflen isn't total to write, it's total available. */
+    /* outbuflen isn't total to write, it's total available. */
+    const int wantedoutframes = ((Sint64) inframes) * outrate / inrate;
     const int maxoutframes = outbuflen / framelen;
     const int outframes = SDL_min(wantedoutframes, maxoutframes);
-    ResampleFloatType outtime = 0.0f;
     float *dst = outbuf;
     int i, j, chan;
 
     for (i = 0; i < outframes; i++) {
-        const int srcindex = (int) (outtime * inrate);
+        const int srcindex = ((Sint64) i) * inrate / outrate;
-        const ResampleFloatType intime = ((ResampleFloatType) srcindex) / finrate;
+        /* Calculating the following way avoids subtraction or modulo of large
-        const ResampleFloatType innexttime = ((ResampleFloatType) (srcindex + 1)) / finrate;
+         * floats which have low result precision.
-        const ResampleFloatType indeltatime = innexttime - intime;
+         *   interpolation1
-        const ResampleFloatType interpolation1 = (indeltatime == 0.0f) ? 1.0f : (1.0f - ((innexttime - outtime) / indeltatime));
+         * = (i / outrate * inrate) - floor(i / outrate * inrate)
-        const int filterindex1 = (int) (interpolation1 * RESAMPLER_SAMPLES_PER_ZERO_CROSSING);
+         * = mod(i / outrate * inrate, 1)
-        const ResampleFloatType interpolation2 = 1.0f - interpolation1;
+         * = mod(i * inrate, outrate) / outrate */
-        const int filterindex2 = (int) (interpolation2 * RESAMPLER_SAMPLES_PER_ZERO_CROSSING);
+        const int srcfraction = ((Sint64) i) * inrate % outrate;
+        const float interpolation1 = ((float) srcfraction) / ((float) outrate);
+        const int filterindex1 = ((Sint32) srcfraction) * RESAMPLER_SAMPLES_PER_ZERO_CROSSING / outrate;
+        const float interpolation2 = 1.0f - interpolation1;
+        const int filterindex2 = ((Sint32) (outrate - srcfraction)) * RESAMPLER_SAMPLES_PER_ZERO_CROSSING / outrate;
 
         for (chan = 0; chan < chans; chan++) {
             float outsample = 0.0f;
 
             /* do this twice to calculate the sample, once for the "left wing" and then same for the right. */
             for (j = 0; (filterindex1 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) < RESAMPLER_FILTER_SIZE; j++) {
+                const int filt_ind = filterindex1 + j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING;
                 const int srcframe = srcindex - j;
                 /* !!! FIXME: we can bubble this conditional out of here by doing a pre loop. */
                 const float insample = (srcframe < 0) ? lpadding[((paddinglen + srcframe) * chans) + chan] : inbuf[(srcframe * chans) + chan];
-                outsample += (float)(insample * (ResamplerFilter[filterindex1 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)] + (interpolation1 * ResamplerFilterDifference[filterindex1 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)])));
+                outsample += (float)(insample * (ResamplerFilter[filt_ind] + (interpolation1 * ResamplerFilterDifference[filt_ind])));
             }
 
             /* Do the right wing! */
             for (j = 0; (filterindex2 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) < RESAMPLER_FILTER_SIZE; j++) {
-                const int jsamples = j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING;
+                const int filt_ind = filterindex2 + j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING;
                 const int srcframe = srcindex + 1 + j;
                 /* !!! FIXME: we can bubble this conditional out of here by doing a post loop. */
                 const float insample = (srcframe >= inframes) ? rpadding[((srcframe - inframes) * chans) + chan] : inbuf[(srcframe * chans) + chan];
-                outsample += (float)(insample * (ResamplerFilter[filterindex2 + jsamples] + (interpolation2 * ResamplerFilterDifference[filterindex2 + jsamples])));
+                outsample += (float)(insample * (ResamplerFilter[filt_ind] + (interpolation2 * ResamplerFilterDifference[filt_ind])));
             }
 
             *(dst++) = outsample;
         }
-
-        outtime = ((ResampleFloatType) i) / ((ResampleFloatType) outrate);
     }
 
     return outframes * chans * sizeof (float);

test/testautomation_audio.c

@@ -8,6 +8,7 @@
 # define _CRT_SECURE_NO_WARNINGS
 #endif
 
+#include <math.h>
 #include <stdio.h>
 #include <string.h>
 
@@ -969,6 +970,118 @@ int audio_openCloseAudioDeviceConnected()
    return TEST_COMPLETED;
 }
 
+static double sine_wave_sample(const Sint64 idx, const Sint64 rate, const Sint64 freq, const double phase)
+{
+  /* Using integer modulo to avoid precision loss caused by large floating
+   * point numbers. Sint64 is needed for the large integer multiplication.
+   * The integers are assumed to be non-negative so that modulo is always
+   * non-negative.
+   *   sin(i / rate * freq * 2 * M_PI + phase)
+   * = sin(mod(i / rate * freq, 1) * 2 * M_PI + phase)
+   * = sin(mod(i * freq, rate) / rate * 2 * M_PI + phase) */
+  return SDL_sin(((double) (idx * freq % rate)) / ((double) rate) * (M_PI * 2) + phase);
+}
+
+/**
+ * \brief Check signal-to-noise ratio and maximum error of audio resampling.
+ *
+ * \sa https://wiki.libsdl.org/SDL_BuildAudioCVT
+ * \sa https://wiki.libsdl.org/SDL_ConvertAudio
+ */
+int audio_resampleLoss()
+{
+  /* Note: always test long input time (>= 5s from experience) in some test
+   * cases because an improper implementation may suffer from low resampling
+   * precision with long input due to e.g. doing subtraction with large floats. */
+  struct test_spec_t {
+    int time;
+    int freq;
+    double phase;
+    int rate_in;
+    int rate_out;
+    double signal_to_noise;
+    double max_error;
+  } test_specs[] = {
+    { 50, 440, 0, 44100, 48000, 60, 0.0025 },
+    { 50, 5000, M_PI / 2, 20000, 10000, 65, 0.0010 },
+    { 0 }
+  };
+
+  int spec_idx = 0;
+
+  for (spec_idx = 0; test_specs[spec_idx].time > 0; ++spec_idx) {
+    const struct test_spec_t *spec = &test_specs[spec_idx];
+    const int frames_in = spec->time * spec->rate_in;
+    const int frames_target = spec->time * spec->rate_out;
+    const int len_in = frames_in * (int) sizeof (float);
+    const int len_target = frames_target * (int) sizeof (float);
+
+    Uint64 tick_beg = 0;
+    Uint64 tick_end = 0;
+    SDL_AudioCVT cvt;
+    int i = 0;
+    int ret = 0;
+    double max_error = 0;
+    double sum_squared_error = 0;
+    double sum_squared_value = 0;
+    double signal_to_noise = 0;
+
+    SDLTest_AssertPass("Test resampling of %i s %i Hz %f phase sine wave from sampling rate of %i Hz to %i Hz",
+                       spec->time, spec->freq, spec->phase, spec->rate_in, spec->rate_out);
+
+    ret = SDL_BuildAudioCVT(&cvt, AUDIO_F32, 1, spec->rate_in, AUDIO_F32, 1, spec->rate_out);
+    SDLTest_AssertPass("Call to SDL_BuildAudioCVT(&cvt, AUDIO_F32, 1, %i, AUDIO_F32, 1, %i)", spec->rate_in, spec->rate_out);
+    SDLTest_AssertCheck(ret == 1, "Expected SDL_BuildAudioCVT to succeed and conversion to be needed.");
+    if (ret != 1) {
+      return TEST_ABORTED;
+    }
+
+    cvt.buf = (Uint8 *) SDL_malloc(len_in * cvt.len_mult);
+    SDLTest_AssertCheck(cvt.buf != NULL, "Expected input buffer to be created.");
+    if (cvt.buf == NULL) {
+      return TEST_ABORTED;
+    }
+
+    cvt.len = len_in;
+    for (i = 0; i < frames_in; ++i) {
+      *(((float *) cvt.buf) + i) = (float) sine_wave_sample(i, spec->rate_in, spec->freq, spec->phase);
+    }
+
+    tick_beg = SDL_GetPerformanceCounter();
+    ret = SDL_ConvertAudio(&cvt);
+    tick_end = SDL_GetPerformanceCounter();
+    SDLTest_AssertPass("Call to SDL_ConvertAudio(&cvt)");
+    SDLTest_AssertCheck(ret == 0, "Expected SDL_ConvertAudio to succeed.");
+    SDLTest_AssertCheck(cvt.len_cvt == len_target, "Expected output length %i, got %i.", len_target, cvt.len_cvt);
+    if (ret != 0 || cvt.len_cvt != len_target) {
+      SDL_free(cvt.buf);
+      return TEST_ABORTED;
+    }
+    SDLTest_Log("Resampling used %f seconds.", ((double) (tick_end - tick_beg)) / SDL_GetPerformanceFrequency());
+
+    for (i = 0; i < frames_target; ++i) {
+        const float output = *(((float *) cvt.buf) + i);
+        const double target = sine_wave_sample(i, spec->rate_out, spec->freq, spec->phase);
+        const double error = SDL_fabs(target - output);
+        max_error = SDL_max(max_error, error);
+        sum_squared_error += error * error;
+        sum_squared_value += target * target;
+    }
+    SDL_free(cvt.buf);
+    signal_to_noise = 10 * SDL_log10(sum_squared_value / sum_squared_error); /* decibel */
+    SDLTest_AssertCheck(isfinite(sum_squared_value), "Sum of squared target should be finite.");
+    SDLTest_AssertCheck(isfinite(sum_squared_error), "Sum of squared error should be finite.");
+    /* Infinity is theoretically possible when there is very little to no noise */
+    SDLTest_AssertCheck(!isnan(signal_to_noise), "Signal-to-noise ratio should not be NaN.");
+    SDLTest_AssertCheck(isfinite(max_error), "Maximum conversion error should be finite.");
+    SDLTest_AssertCheck(signal_to_noise >= spec->signal_to_noise, "Conversion signal-to-noise ratio %f dB should be no less than %f dB.",
+                        signal_to_noise, spec->signal_to_noise);
+    SDLTest_AssertCheck(max_error <= spec->max_error, "Maximum conversion error %f should be no more than %f.",
+                        max_error, spec->max_error);
+  }
+
+  return TEST_COMPLETED;
+}
 
 
 /* ================= Test Case References ================== */
@@ -1024,11 +1137,14 @@ static const SDLTest_TestCaseReference audioTest14 =
 static const SDLTest_TestCaseReference audioTest15 =
         { (SDLTest_TestCaseFp)audio_pauseUnpauseAudio, "audio_pauseUnpauseAudio", "Pause and Unpause audio for various audio specs while testing callback.", TEST_ENABLED };
 
+static const SDLTest_TestCaseReference audioTest16 =
+        { (SDLTest_TestCaseFp)audio_resampleLoss, "audio_resampleLoss", "Check signal-to-noise ratio and maximum error of audio resampling.", TEST_ENABLED };
+
 /* Sequence of Audio test cases */
 static const SDLTest_TestCaseReference *audioTests[] =  {
     &audioTest1, &audioTest2, &audioTest3, &audioTest4, &audioTest5, &audioTest6,
     &audioTest7, &audioTest8, &audioTest9, &audioTest10, &audioTest11,
-    &audioTest12, &audioTest13, &audioTest14, &audioTest15, NULL
+    &audioTest12, &audioTest13, &audioTest14, &audioTest15, &audioTest16, NULL
 };
 
 /* Audio test suite (global) */