From 20e17559e545c5d3cfe86c1c4772365e70090779 Mon Sep 17 00:00:00 2001 From: Qrox Date: Thu, 9 Mar 2023 17:34:51 +0800 Subject: [PATCH] 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 | 56 +++++++++-------- test/testautomation_audio.c | 118 +++++++++++++++++++++++++++++++++++- 2 files changed, 148 insertions(+), 26 deletions(-) diff --git a/src/audio/SDL_audiocvt.c b/src/audio/SDL_audiocvt.c index de97beacc..a30ae5c9e 100644 --- a/src/audio/SDL_audiocvt.c +++ b/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 -ResamplerPadding(const int inrate, const int outrate) +static Sint32 +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 - almost twice the speed on Intel processors, and orders of magnitude more - on CPUs that need a software fallback for double calculations. */ - typedef float ResampleFloatType; - - const ResampleFloatType finrate = (ResampleFloatType) inrate; - const ResampleFloatType ratio = ((float) outrate) / ((float) inrate); + /* This function uses integer arithmetics to avoid precision loss caused + * by large floating point numbers. For some operations, Sint32 or Sint64 + * are needed for the large number multiplications. The input integers are + * assumed to be non-negative so that division rounds by truncation and + * modulo is always non-negative. Note that the operator order is important + * for these integer divisions. */ 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 ResampleFloatType intime = ((ResampleFloatType) srcindex) / finrate; - const ResampleFloatType innexttime = ((ResampleFloatType) (srcindex + 1)) / finrate; - const ResampleFloatType indeltatime = innexttime - intime; - const ResampleFloatType interpolation1 = (indeltatime == 0.0f) ? 1.0f : (1.0f - ((innexttime - outtime) / indeltatime)); - const int filterindex1 = (int) (interpolation1 * RESAMPLER_SAMPLES_PER_ZERO_CROSSING); - const ResampleFloatType interpolation2 = 1.0f - interpolation1; - const int filterindex2 = (int) (interpolation2 * RESAMPLER_SAMPLES_PER_ZERO_CROSSING); + const int srcindex = ((Sint64) i) * inrate / outrate; + /* Calculating the following way avoids subtraction or modulo of large + * floats which have low result precision. + * interpolation1 + * = (i / outrate * inrate) - floor(i / outrate * inrate) + * = mod(i / outrate * inrate, 1) + * = mod(i * inrate, outrate) / outrate */ + 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); diff --git a/test/testautomation_audio.c b/test/testautomation_audio.c index e6127b573..22ee2e6b3 100644 --- a/test/testautomation_audio.c +++ b/test/testautomation_audio.c @@ -8,6 +8,7 @@ # define _CRT_SECURE_NO_WARNINGS #endif +#include #include #include @@ -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) */