Skip to content

Commit

Permalink
Uses integer arithmetics in SDL_ResampleAudio
Browse files Browse the repository at this point in the history 
- Revert resampler workaround
- 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
  • Loading branch information
@Qrox @slouken
Qrox authored and slouken committed Mar 10, 2023
1 parent 926db1b commit 71c80c4
Show file tree
Hide file tree
Showing 2 changed files with 183 additions and 33 deletions.
62 changes: 30 additions & 32 deletions src/audio/SDL_audiocvt.c
View file
Expand Up @@ -261,13 +261,17 @@ static void SDLCALL SDL_ConvertMonoToStereo_SSE(SDL_AudioCVT *cvt, SDL_AudioForm

#include "SDL_audio_resampler_filter.h"

static int GetResamplerPadding(const int inrate, const int outrate)
static Sint32 GetResamplerPadding(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;
}
Expand All @@ -278,65 +282,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 = GetResamplerPadding(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);
Expand Down
154 changes: 153 additions & 1 deletion test/testautomation_audio.c
View file
Expand Up @@ -8,6 +8,7 @@
#define _CRT_SECURE_NO_WARNINGS
#endif

#include <math.h>
#include <stdio.h>

#include <SDL3/SDL.h>
Expand Down Expand Up @@ -989,6 +990,153 @@ static int audio_openCloseAudioDeviceConnected(void *arg)
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 * PI + phase)
* = sin(mod(i / rate * freq, 1) * 2 * PI + phase)
* = sin(mod(i * freq, rate) / rate * 2 * PI + phase) */
return SDL_sin(((double)(idx * freq % rate)) / ((double)rate) * (SDL_PI_D * 2) + phase);
}

/**
* \brief Check signal-to-noise ratio and maximum error of audio resampling.
*
* \sa https://wiki.libsdl.org/SDL_CreateAudioStream
* \sa https://wiki.libsdl.org/SDL_DestroyAudioStream
* \sa https://wiki.libsdl.org/SDL_PutAudioStreamData
* \sa https://wiki.libsdl.org/SDL_FlushAudioStream
* \sa https://wiki.libsdl.org/SDL_GetAudioStreamData
*/
static int audio_resampleLoss(void *arg)
{
/* 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, SDL_PI_D / 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;
int i = 0;
int ret = 0;
SDL_AudioStream *stream = NULL;
float *buf_in = NULL;
float *buf_out = NULL;
int len_out = 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);

stream = SDL_CreateAudioStream(AUDIO_F32, 1, spec->rate_in, AUDIO_F32, 1, spec->rate_out);
SDLTest_AssertPass("Call to SDL_CreateAudioStream(AUDIO_F32, 1, %i, AUDIO_F32, 1, %i)", spec->rate_in, spec->rate_out);
SDLTest_AssertCheck(stream != NULL, "Expected SDL_CreateAudioStream to succeed.");
if (stream == NULL) {
return TEST_ABORTED;
}

buf_in = (float *)SDL_malloc(len_in);
SDLTest_AssertCheck(buf_in != NULL, "Expected input buffer to be created.");
if (buf_in == NULL) {
SDL_DestroyAudioStream(stream);
return TEST_ABORTED;
}

for (i = 0; i < frames_in; ++i) {
*(buf_in + i) = (float)sine_wave_sample(i, spec->rate_in, spec->freq, spec->phase);
}

tick_beg = SDL_GetPerformanceCounter();

ret = SDL_PutAudioStreamData(stream, buf_in, len_in);
SDLTest_AssertPass("Call to SDL_PutAudioStreamData(stream, buf_in, %i)", len_in);
SDLTest_AssertCheck(ret == 0, "Expected SDL_PutAudioStreamData to succeed.");
SDL_free(buf_in);
if (ret != 0) {
SDL_DestroyAudioStream(stream);
return TEST_ABORTED;
}

ret = SDL_FlushAudioStream(stream);
SDLTest_AssertPass("Call to SDL_FlushAudioStream(stream)");
SDLTest_AssertCheck(ret == 0, "Expected SDL_FlushAudioStream to succeed");
if (ret != 0) {
SDL_DestroyAudioStream(stream);
return TEST_ABORTED;
}

buf_out = (float *)SDL_malloc(len_target);
SDLTest_AssertCheck(buf_out != NULL, "Expected output buffer to be created.");
if (buf_out == NULL) {
SDL_DestroyAudioStream(stream);
return TEST_ABORTED;
}

len_out = SDL_GetAudioStreamData(stream, buf_out, len_target);
SDLTest_AssertPass("Call to SDL_GetAudioStreamData(stream, buf_out, %i)", len_target);
/** !!! FIXME: SDL_AudioStream does not return output of the same length as
** !!! FIXME: the input even if SDL_FlushAudioStream is called. */
SDLTest_AssertCheck(len_out <= len_target, "Expected output length to be no larger than %i, got %i.",
len_target, len_out);
SDL_DestroyAudioStream(stream);
if (len_out > len_target) {
SDL_free(buf_out);
return TEST_ABORTED;
}

tick_end = SDL_GetPerformanceCounter();
SDLTest_Log("Resampling used %f seconds.", ((double)(tick_end - tick_beg)) / SDL_GetPerformanceFrequency());

for (i = 0; i < len_out / (int)sizeof(float); ++i) {
const float output = *(buf_out + 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(buf_out);
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 */
Expand Down Expand Up @@ -1058,11 +1206,15 @@ static const SDLTest_TestCaseReference audioTest15 = {
audio_pauseUnpauseAudio, "audio_pauseUnpauseAudio", "Pause and Unpause audio for various audio specs while testing callback.", TEST_ENABLED
};

static const SDLTest_TestCaseReference audioTest16 = {
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) */
Expand Down

0 comments on commit 71c80c4

Please sign in to comment.