diff --git a/src/audio/SDL_audiocvt.c b/src/audio/SDL_audiocvt.c index 7fde2b93c..3dc60183c 100644 --- a/src/audio/SDL_audiocvt.c +++ b/src/audio/SDL_audiocvt.c @@ -377,10 +377,24 @@ SDL_Convert51To71(SDL_AudioCVT * cvt, SDL_AudioFormat format) /* SDL's resampler uses a "bandlimited interpolation" algorithm: https://ccrma.stanford.edu/~jos/resample/ */ -#define RESAMPLER_ZERO_CROSSINGS 5 +#define RESAMPLER_ZERO_CROSSINGS 10 #define RESAMPLER_BITS_PER_SAMPLE 16 #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) +/* +This is ResamplerFilter's "native" cutoff frequency, expressed as a fraction of the sample rate we apply the kernel to (so, between 0 and 0.5). +This would be the cutoff frequency if we used every value in ResamplerFilter without any lerping / skipping array elements. +*/ +#define RESAMPLER_NATIVE_CUTOFF_FREQ (1.0 / (2.0 * (double)RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) +/* +Desired bandwidth (between 0 and 1), used to set the cutoff frequency of the interpolated filter kernel, +where 1 means the filter cutoff is at the nyquist rate. + +This is used to make the filter start rolling off before the nyquist rate of the dest sample rate, +to suppress aliasing. +*/ +#define RESAMPLER_BANDWIDTH 0.80 + /* This is a "modified" bessel function, so you can't use POSIX j0() */ static double @@ -426,11 +440,91 @@ kaiser_and_sinc(float *table, float *diffs, const int tablelen, const double bet diffs[lenm1] = 0.0f; } - static SDL_SpinLock ResampleFilterSpinlock = 0; static float *ResamplerFilter = NULL; static float *ResamplerFilterDifference = NULL; +/* + Returns a value from a filter kernel that performs lowpass at f_c. + You must make repeated calls with i=0, 1, 2, ... to get the whole kernel. + + f_c is a fraction of the sample rate of the data with which you convolve the returned + filter kernel. (0 < f_c <= 0.5, since it can't be greater than the nyquist frequency) + + as a special feature, you can request `i` at a fractional index. + This phase-shifts the kernel by between 0 and 1 samples. + + returns an approximation of: (ignoring the window function in ResamplerFilter): + const float x = f_c * 2.0 * i * M_PI + return sin(x)/x + + */ +static inline float +kernel_sample(double f_c, float i) +{ + /* + ResamplerFilter[i] has a cutoff of RESAMPLER_NATIVE_CUTOFF_FREQ + ResamplerFilter[(f_c/RESAMPLER_NATIVE_CUTOFF_FREQ) * i] has a cutoff of f_c + + next, we want to evaluate ResamplerFilter at (scale*i), which is a float. + + in practice we can only evaluate it at integers, so we'll do linear interpolation + between adjacent kernel entries. + */ + const double scale = (f_c / RESAMPLER_NATIVE_CUTOFF_FREQ); + const double filter_index_float = scale * i; + const int filter_index_int = (int)filter_index_float; + const double filter_index_fraction = (filter_index_float - filter_index_int); /* between 0 and 1 */ + float lerped; + + // FIXME: this is currently necessary, which means kernel_size() is probably wrong + if (filter_index_int >= RESAMPLER_FILTER_SIZE) + return 0; + + lerped = ResamplerFilter[filter_index_int] + (filter_index_fraction * ResamplerFilterDifference[filter_index_int]); + return lerped; +} + +static inline int +kernel_size(double desired_f_c) +{ + const double scale = (desired_f_c / RESAMPLER_NATIVE_CUTOFF_FREQ); + + /* for: + const double filter_index_float = scale * i; + + the max element we need in ResamplerFilter is: + int max_index = floor(filter_index_float) + 1; + + solve for max_index == (RESAMPLER_FILTER_SIZE - 1): + + RESAMPLER_FILTER_SIZE - 2 = floor(filter_index_float) + RESAMPLER_FILTER_SIZE - 2 = floor(scale * i) + RESAMPLER_FILTER_SIZE - 2 + R = scale * i, set R == 1 + (RESAMPLER_FILTER_SIZE - 1) / scale = i + */ + + double maxarg = (double)(RESAMPLER_FILTER_SIZE - 1) / scale; + maxarg = ceil(maxarg); + return (int)maxarg; +} + +float kernel_sum(double f_c) +{ + const int kernelsize = kernel_size(f_c); + float sum = 0; + int i; + + /* center point */ + sum = kernel_sample(f_c, 0); + + /* do left and right wings in one loop */ + for (i = 1; i <= kernelsize; i++) { + sum += 2 * kernel_sample(f_c, i); + } + return sum; +} + int SDL_PrepareResampleFilter(void) { @@ -474,12 +568,14 @@ ResamplerPadding(const int inrate, const int outrate) { if (inrate == outrate) { return 0; - } else if (inrate > outrate) { - return (int) SDL_ceil(((float) (RESAMPLER_SAMPLES_PER_ZERO_CROSSING * inrate) / ((float) outrate))); } - return RESAMPLER_SAMPLES_PER_ZERO_CROSSING; + // FIXME: must match below + /* filter cutoff is a fraction of the input samplerate */ + const double f_c = (RESAMPLER_BANDWIDTH * 0.5 * SDL_min(inrate, outrate)) / (double)inrate; + + /* +1 because "const int srcframe = srcindex + 1 + j;" */ + return kernel_size(f_c) + 1; } - /* lpadding and rpadding are expected to be buffers of (ResamplePadding(inrate, outrate) * chans * sizeof (float)) bytes. */ static int SDL_ResampleAudio(const int chans, const int inrate, const int outrate, @@ -487,50 +583,50 @@ SDL_ResampleAudio(const int chans, const int inrate, const int outrate, const float *inbuf, const int inbuflen, float *outbuf, const int outbuflen) { - const double finrate = (double) inrate; - const double outtimeincr = 1.0 / ((float) outrate); - const double ratio = ((float) outrate) / ((float) inrate); + const double inv_ratio = ((double) inrate) / ((double) outrate); + const double ratio = ((double) outrate) / ((double) 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. */ const int maxoutframes = outbuflen / framelen; const int outframes = SDL_min(wantedoutframes, maxoutframes); + const double f_c = (RESAMPLER_BANDWIDTH * 0.5 * SDL_min(inrate, outrate)) / (double)inrate; /* desired filter cutoff as a fraction of the input samplerate */ + const int filtsize = kernel_size(f_c); + const float kernelsum = kernel_sum(f_c); float *dst = outbuf; - double outtime = 0.0; int i, j, chan; - + for (i = 0; i < outframes; i++) { - const int srcindex = (int) (outtime * inrate); - const double intime = ((double) srcindex) / finrate; - const double innexttime = ((double) (srcindex + 1)) / finrate; - const double interpolation1 = 1.0 - ((innexttime - outtime) / (innexttime - intime)); - const int filterindex1 = (int) (interpolation1 * RESAMPLER_SAMPLES_PER_ZERO_CROSSING); - const double interpolation2 = 1.0 - interpolation1; - const int filterindex2 = (int) (interpolation2 * RESAMPLER_SAMPLES_PER_ZERO_CROSSING); + const double outtime = i * inv_ratio; /* this is in fractional input samples */ + const int srcindex = (int)floor(outtime); /* integer input same index on or before the output sample */ + const double interpolation = outtime - floor(outtime); /* in [0..1) */ 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. */ /* !!! FIXME: do both wings in one loop */ - for (j = 0; (filterindex1 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) < RESAMPLER_FILTER_SIZE; j++) { + for (j = 0; j < filtsize; j++) { 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)]))); + const float kernel_coeff = kernel_sample(f_c, (float)j + interpolation); + outsample += insample * kernel_coeff; } - for (j = 0; (filterindex2 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) < RESAMPLER_FILTER_SIZE; j++) { + for (j = 0; j < filtsize; j++) { 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 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)] + (interpolation2 * ResamplerFilterDifference[filterindex2 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)]))); + const float kernel_coeff = kernel_sample(f_c, (float)j + (1.0 - interpolation)); + outsample += insample * kernel_coeff; } + + /* scale the kernel so it sums to 1 */ + outsample /= kernelsum; *(dst++) = outsample; } - - outtime += outtimeincr; } return outframes * chans * sizeof (float);