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

(cherry-picked from commit 20e17559e545c5d3cfe86c1c4772365e70090779)
This commit is contained in:
Qrox 2023-03-09 17:34:51 +08:00 committed by Ozkan Sezer
parent 5e35309913
commit ef347776c2
2 changed files with 149 additions and 33 deletions

View file

@ -181,13 +181,17 @@ static void SDLCALL SDL_ConvertMonoToStereo_SSE(SDL_AudioCVT *cvt, SDL_AudioForm
#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;
}
@ -198,65 +202,59 @@ static int SDL_ResampleAudio(const int chans, const int inrate, const int outrat
const float *inbuf, const int inbuflen,
float *outbuf, const int outbuflen)
{
/* !!! FIXME: this produces artifacts if we don't work at double precision, but this turns out to
be a big performance hit. Until we can resolve this better, we force this to double
for amd64 CPUs, which should be able to take the hit for now, vs small embedded
things that might end up in a software fallback here. */
/* 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. */
#if defined(_M_X64) || defined(__x86_64__)
typedef double ResampleFloatType;
#else
typedef float ResampleFloatType;
#endif
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)(inframes * 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);

View file

@ -8,6 +8,7 @@
#define _CRT_SECURE_NO_WARNINGS
#endif
#include <math.h>
#include <stdio.h>
#include <string.h>
@ -965,6 +966,119 @@ 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 ================== */
/* Audio test cases */
@ -1033,11 +1147,15 @@ 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) */