Improved ResampleAudio

* filterindex2 was off-by-one
* Generate ResamplerFilter using doubles
* Transpose ResamplerFilter to improve access patterns
main
Brick 2023-08-19 20:12:40 +01:00 committed by Ryan C. Gordon
parent cdaa19869d
commit b9541b9eab
4 changed files with 1065 additions and 1066 deletions

View File

@ -41,23 +41,27 @@ gcc -o genfilter build-scripts/gen_audio_resampler_filter.c -lm && ./genfilter >
#include <stdio.h> #include <stdio.h>
#include <math.h> #include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define RESAMPLER_ZERO_CROSSINGS 5 #define RESAMPLER_ZERO_CROSSINGS 5
#define RESAMPLER_BITS_PER_SAMPLE 16 #define RESAMPLER_BITS_PER_SAMPLE 16
#define RESAMPLER_SAMPLES_PER_ZERO_CROSSING (1 << ((RESAMPLER_BITS_PER_SAMPLE / 2) + 1)) #define RESAMPLER_SAMPLES_PER_ZERO_CROSSING (1 << ((RESAMPLER_BITS_PER_SAMPLE / 2) + 1))
#define RESAMPLER_FILTER_SIZE ((RESAMPLER_SAMPLES_PER_ZERO_CROSSING * RESAMPLER_ZERO_CROSSINGS) + 1) #define RESAMPLER_FILTER_SIZE (RESAMPLER_SAMPLES_PER_ZERO_CROSSING * RESAMPLER_ZERO_CROSSINGS)
/* This is a "modified" bessel function, so you can't use POSIX j0() */ /* This is a "modified" bessel function, so you can't use POSIX j0() */
static double static double
bessel(const double x) bessel(const double x)
{ {
const double xdiv2 = x / 2.0; const double xdiv2 = x / 2.0;
double i0 = 1.0f; double i0 = 1.0;
double f = 1.0f; double f = 1.0;
int i = 1; int i = 1;
while (1) { while (1) {
const double diff = pow(xdiv2, i * 2) / pow(f, 2); const double diff = pow(xdiv2, i * 2) / pow(f, 2);
if (diff < 1.0e-21f) { if (diff < 1.0e-21) {
break; break;
} }
i0 += diff; i0 += diff;
@ -70,30 +74,25 @@ bessel(const double x)
/* build kaiser table with cardinal sine applied to it, and array of differences between elements. */ /* build kaiser table with cardinal sine applied to it, and array of differences between elements. */
static void static void
kaiser_and_sinc(float *table, float *diffs, const int tablelen, const double beta) kaiser_and_sinc(double *table, double *diffs, const int tablelen, const double beta)
{ {
const int lenm1 = tablelen - 1;
const int lenm1div2 = lenm1 / 2;
const double bessel_beta = bessel(beta); const double bessel_beta = bessel(beta);
int i; int i;
table[0] = 1.0f; table[0] = 1.0;
for (i = 1; i < tablelen; i++) { diffs[tablelen - 1] = 0.0;
const double kaiser = bessel(beta * sqrt(1.0 - pow(((i - lenm1) / 2.0) / lenm1div2, 2.0))) / bessel_beta;
table[tablelen - i] = (float) kaiser;
}
for (i = 1; i < tablelen; i++) { for (i = 1; i < tablelen; i++) {
const float x = (((float) i) / ((float) RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) * ((float) M_PI); const double kaiser = bessel(beta * sqrt(1.0 - pow((double)i / (double)(tablelen - 1), 2.0))) / bessel_beta;
table[i] *= sinf(x) / x; const double x = (((double) i) / ((double) RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) * M_PI;
table[i] = kaiser * (sin(x) / x);
diffs[i - 1] = table[i] - table[i - 1]; diffs[i - 1] = table[i] - table[i - 1];
} }
diffs[lenm1] = 0.0f;
} }
static float ResamplerFilter[RESAMPLER_FILTER_SIZE]; static double ResamplerFilter[RESAMPLER_FILTER_SIZE + 1];
static float ResamplerFilterDifference[RESAMPLER_FILTER_SIZE]; static double ResamplerFilterDifference[RESAMPLER_FILTER_SIZE + 1];
static void static void
PrepareResampleFilter(void) PrepareResampleFilter(void)
@ -101,12 +100,12 @@ PrepareResampleFilter(void)
/* if dB > 50, beta=(0.1102 * (dB - 8.7)), according to Matlab. */ /* if dB > 50, beta=(0.1102 * (dB - 8.7)), according to Matlab. */
const double dB = 80.0; const double dB = 80.0;
const double beta = 0.1102 * (dB - 8.7); const double beta = 0.1102 * (dB - 8.7);
kaiser_and_sinc(ResamplerFilter, ResamplerFilterDifference, RESAMPLER_FILTER_SIZE, beta); kaiser_and_sinc(ResamplerFilter, ResamplerFilterDifference, RESAMPLER_FILTER_SIZE + 1, beta);
} }
int main(void) int main(void)
{ {
int i; int i, j;
PrepareResampleFilter(); PrepareResampleFilter();
@ -137,21 +136,21 @@ int main(void)
"#define RESAMPLER_ZERO_CROSSINGS %d\n" "#define RESAMPLER_ZERO_CROSSINGS %d\n"
"#define RESAMPLER_BITS_PER_SAMPLE %d\n" "#define RESAMPLER_BITS_PER_SAMPLE %d\n"
"#define RESAMPLER_SAMPLES_PER_ZERO_CROSSING (1 << ((RESAMPLER_BITS_PER_SAMPLE / 2) + 1))\n" "#define RESAMPLER_SAMPLES_PER_ZERO_CROSSING (1 << ((RESAMPLER_BITS_PER_SAMPLE / 2) + 1))\n"
"#define RESAMPLER_FILTER_SIZE ((RESAMPLER_SAMPLES_PER_ZERO_CROSSING * RESAMPLER_ZERO_CROSSINGS) + 1)\n" "#define RESAMPLER_FILTER_SIZE (RESAMPLER_SAMPLES_PER_ZERO_CROSSING * RESAMPLER_ZERO_CROSSINGS)\n"
"\n", RESAMPLER_ZERO_CROSSINGS, RESAMPLER_BITS_PER_SAMPLE "\n", RESAMPLER_ZERO_CROSSINGS, RESAMPLER_BITS_PER_SAMPLE
); );
printf("static const float ResamplerFilter[RESAMPLER_FILTER_SIZE] = {\n"); printf("static const float ResamplerFilter[RESAMPLER_FILTER_SIZE] = {");
printf(" %.9ff", ResamplerFilter[0]); for (i = 0; i < RESAMPLER_FILTER_SIZE; i++) {
for (i = 0; i < RESAMPLER_FILTER_SIZE-1; i++) { j = (i % RESAMPLER_ZERO_CROSSINGS) * RESAMPLER_SAMPLES_PER_ZERO_CROSSING + (i / RESAMPLER_ZERO_CROSSINGS);
printf("%s%.9ff", ((i % 5) == 4) ? ",\n " : ", ", ResamplerFilter[i+1]); printf("%s%12.9ff,", (i % RESAMPLER_ZERO_CROSSINGS) ? "" : "\n ", ResamplerFilter[j]);
} }
printf("\n};\n\n"); printf("\n};\n\n");
printf("static const float ResamplerFilterDifference[RESAMPLER_FILTER_SIZE] = {\n"); printf("static const float ResamplerFilterDifference[RESAMPLER_FILTER_SIZE] = {");
printf(" %.9ff", ResamplerFilterDifference[0]); for (i = 0; i < RESAMPLER_FILTER_SIZE; i++) {
for (i = 0; i < RESAMPLER_FILTER_SIZE-1; i++) { j = (i % RESAMPLER_ZERO_CROSSINGS) * RESAMPLER_SAMPLES_PER_ZERO_CROSSING + (i / RESAMPLER_ZERO_CROSSINGS);
printf("%s%.9ff", ((i % 5) == 4) ? ",\n " : ", ", ResamplerFilterDifference[i+1]); printf("%s%12.9ff,", (i % RESAMPLER_ZERO_CROSSINGS) ? "" : "\n ", ResamplerFilterDifference[j]);
} }
printf("\n};\n\n"); printf("\n};\n\n");

File diff suppressed because it is too large Load Diff

View File

@ -95,16 +95,16 @@ static void ResampleAudio(const int chans, const int inrate, const int outrate,
* = mod(i * inrate, outrate) / outrate */ * = mod(i * inrate, outrate) / outrate */
const int srcfraction = (int)(srcpos % outrate); const int srcfraction = (int)(srcpos % outrate);
const float interpolation1 = ((float)srcfraction) / ((float)outrate); const float interpolation1 = ((float)srcfraction) / ((float)outrate);
const int filterindex1 = ((Sint32)srcfraction) * RESAMPLER_SAMPLES_PER_ZERO_CROSSING / outrate; const int filterindex1 = (((Sint32)srcfraction) * RESAMPLER_SAMPLES_PER_ZERO_CROSSING / outrate) * RESAMPLER_ZERO_CROSSINGS;
const float interpolation2 = 1.0f - interpolation1; const float interpolation2 = 1.0f - interpolation1;
const int filterindex2 = ((Sint32)(outrate - srcfraction)) * RESAMPLER_SAMPLES_PER_ZERO_CROSSING / outrate; const int filterindex2 = RESAMPLER_FILTER_SIZE - filterindex1 - RESAMPLER_ZERO_CROSSINGS;
for (chan = 0; chan < chans; chan++) { for (chan = 0; chan < chans; chan++) {
float outsample = 0.0f; float outsample = 0.0f;
// do this twice to calculate the sample, once for the "left wing" and then same for the right. // 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++) { for (j = 0; j < RESAMPLER_ZERO_CROSSINGS; j++) {
const int filt_ind = filterindex1 + j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING; const int filt_ind = j + filterindex1;
const int srcframe = srcindex - j; const int srcframe = srcindex - j;
SDL_assert(paddinglen + srcframe >= 0); SDL_assert(paddinglen + srcframe >= 0);
/* !!! FIXME: we can bubble this conditional out of here by doing a pre loop. */ /* !!! FIXME: we can bubble this conditional out of here by doing a pre loop. */
@ -113,8 +113,8 @@ static void ResampleAudio(const int chans, const int inrate, const int outrate,
} }
// Do the right wing! // Do the right wing!
for (j = 0; (filterindex2 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) < RESAMPLER_FILTER_SIZE; j++) { for (j = 0; j < RESAMPLER_ZERO_CROSSINGS; j++) {
const int filt_ind = filterindex2 + j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING; const int filt_ind = j + filterindex2;
const int srcframe = srcindex + 1 + j; const int srcframe = srcindex + 1 + j;
SDL_assert(srcframe - inframes < paddinglen); SDL_assert(srcframe - inframes < paddinglen);
// !!! FIXME: we can bubble this conditional out of here by doing a post loop. // !!! FIXME: we can bubble this conditional out of here by doing a post loop.

View File

@ -794,9 +794,10 @@ static int audio_resampleLoss(void *arg)
double signal_to_noise; double signal_to_noise;
double max_error; double max_error;
} test_specs[] = { } test_specs[] = {
{ 50, 440, 0, 44100, 48000, 60, 0.0025 }, { 50, 440, 0, 44100, 48000, 79, 0.0008 },
{ 50, 5000, SDL_PI_D / 2, 20000, 10000, 65, 0.0010 }, { 50, 5000, SDL_PI_D / 2, 20000, 10000, 999, 0.0001 },
{ 50, 440, 0, 22050, 96000, 60, 0.0120 }, /* I have no idea how to tune these values */ { 50, 440, 0, 22050, 96000, 76, 0.0120 }, /* I have no idea how to tune these values */
{ 50, 440, 0, 96000, 22050, 80, 0.0014 },
{ 0 } { 0 }
}; };