audio: Replaced the resampler. Again.
authorRyan C. Gordon <icculus@icculus.org>
Thu, 21 Sep 2017 02:51:14 -0400
changeset 11508a8382e3d0b54
parent 11507 071deb801690
child 11509 d069a4aefa4e
audio: Replaced the resampler. Again.

This time it's using real math from a real whitepaper instead of my previous
amateur, fast-but-low-quality attempt. The new resampler does "bandlimited
interpolation," as described here: https://ccrma.stanford.edu/~jos/resample/

The output appears to sound cleaner, especially at high frequencies, and of
course works with non-power-of-two rate conversions.

There are some obvious optimizations to be done to this still, and there is
other fallout: this doesn't resample a buffer in-place, the 2-channels-Sint16
fast path is gone because this resampler does a _lot_ of floating point math.
There is a nasty hack to make it work with SDL_AudioCVT.

It's possible these issues are solvable, but they aren't solved as of yet.
Still, I hope this effort is slouching in the right direction.
src/audio/SDL_audio.c
src/audio/SDL_audio_c.h
src/audio/SDL_audiocvt.c
src/audio/kaiser_window.pl
     1.1 --- a/src/audio/SDL_audio.c	Thu Sep 21 14:01:12 2017 +0900
     1.2 +++ b/src/audio/SDL_audio.c	Thu Sep 21 02:51:14 2017 -0400
     1.3 @@ -1543,6 +1543,8 @@
     1.4  #ifdef HAVE_LIBSAMPLERATE_H
     1.5      UnloadLibSampleRate();
     1.6  #endif
     1.7 +
     1.8 +    SDL_FreeResampleFilter();
     1.9  }
    1.10  
    1.11  #define NUM_FORMATS 10
     2.1 --- a/src/audio/SDL_audio_c.h	Thu Sep 21 14:01:12 2017 +0900
     2.2 +++ b/src/audio/SDL_audio_c.h	Thu Sep 21 02:51:14 2017 -0400
     2.3 @@ -69,6 +69,11 @@
     2.4  extern SDL_AudioFilter SDL_Convert_F32_to_U16;
     2.5  extern SDL_AudioFilter SDL_Convert_F32_to_S32;
     2.6  
     2.7 +/* You need to call SDL_PrepareResampleFilter() before using the internal resampler.
     2.8 +   SDL_AudioQuit() calls SDL_FreeResamplerFilter(), you should never call it yourself. */
     2.9 +int SDL_PrepareResampleFilter(void);
    2.10 +void SDL_FreeResampleFilter(void);
    2.11 +
    2.12  
    2.13  /* SDL_AudioStream is a new audio conversion interface. It
    2.14      might eventually become a public API.
     3.1 --- a/src/audio/SDL_audiocvt.c	Thu Sep 21 14:01:12 2017 +0900
     3.2 +++ b/src/audio/SDL_audiocvt.c	Thu Sep 21 02:51:14 2017 -0400
     3.3 @@ -369,228 +369,157 @@
     3.4      }
     3.5  }
     3.6  
     3.7 +/* SDL's resampler uses a "bandlimited interpolation" algorithm:
     3.8 +     https://ccrma.stanford.edu/~jos/resample/ */
     3.9 +
    3.10 +#define RESAMPLER_ZERO_CROSSINGS 5
    3.11 +#define RESAMPLER_BITS_PER_SAMPLE 16
    3.12 +#define RESAMPLER_SAMPLES_PER_ZERO_CROSSING  (1 << ((RESAMPLER_BITS_PER_SAMPLE / 2) + 1))
    3.13 +#define RESAMPLER_FILTER_SIZE ((RESAMPLER_SAMPLES_PER_ZERO_CROSSING * RESAMPLER_ZERO_CROSSINGS) + 1)
    3.14 +
    3.15 +/* This is a "modified" bessel function, so you can't use POSIX j0() */
    3.16 +static double
    3.17 +bessel(const double x)
    3.18 +{
    3.19 +    const double xdiv2 = x / 2.0;
    3.20 +    double i0 = 1.0f;
    3.21 +    double f = 1.0f;
    3.22 +    int i = 1;
    3.23 +
    3.24 +    while (SDL_TRUE) {
    3.25 +        const double diff = SDL_pow(xdiv2, i * 2) / SDL_pow(f, 2);
    3.26 +        if (diff < 1.0e-21f) {
    3.27 +            break;
    3.28 +        }
    3.29 +        i0 += diff;
    3.30 +        i++;
    3.31 +        f *= (double) i;
    3.32 +    }
    3.33 +
    3.34 +    return i0;
    3.35 +}
    3.36 +
    3.37 +/* build kaiser table with cardinal sine applied to it, and array of differences between elements. */
    3.38 +static void
    3.39 +kaiser_and_sinc(float *table, float *diffs, const int tablelen, const double beta)
    3.40 +{
    3.41 +    const int lenm1 = tablelen - 1;
    3.42 +    const int lenm1div2 = lenm1 / 2;
    3.43 +    int i;
    3.44 +
    3.45 +    table[0] = 1.0f;
    3.46 +    for (i = 1; i < tablelen; i++) {
    3.47 +        const double kaiser = bessel(beta * SDL_sqrt(1.0 - SDL_pow(((i - lenm1) / 2.0) / lenm1div2, 2.0))) / bessel(beta);
    3.48 +        table[tablelen - i] = (float) kaiser;
    3.49 +    }
    3.50 +
    3.51 +    for (i = 1; i < tablelen; i++) {
    3.52 +        const float x = (((float) i) / ((float) RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) * ((float) M_PI);
    3.53 +        table[i] *= SDL_sinf(x) / x;
    3.54 +        diffs[i - 1] = table[i] - table[i - 1];
    3.55 +    }
    3.56 +    diffs[lenm1] = 0.0f;
    3.57 +}
    3.58 +
    3.59 +
    3.60 +static SDL_SpinLock ResampleFilterSpinlock = 0;
    3.61 +static float *ResamplerFilter = NULL;
    3.62 +static float *ResamplerFilterDifference = NULL;
    3.63 +
    3.64 +int
    3.65 +SDL_PrepareResampleFilter(void)
    3.66 +{
    3.67 +    SDL_AtomicLock(&ResampleFilterSpinlock);
    3.68 +    if (!ResamplerFilter) {
    3.69 +        /* if dB > 50, beta=(0.1102 * (dB - 8.7)), according to Matlab. */
    3.70 +        const double dB = 80.0;
    3.71 +        const double beta = 0.1102 * (dB - 8.7);
    3.72 +        const size_t alloclen = RESAMPLER_FILTER_SIZE * sizeof (float);
    3.73 +
    3.74 +        ResamplerFilter = (float *) SDL_malloc(alloclen);
    3.75 +        if (!ResamplerFilter) {
    3.76 +            SDL_AtomicUnlock(&ResampleFilterSpinlock);
    3.77 +            return SDL_OutOfMemory();
    3.78 +        }
    3.79 +
    3.80 +        ResamplerFilterDifference = (float *) SDL_malloc(alloclen);
    3.81 +        if (!ResamplerFilterDifference) {
    3.82 +            SDL_free(ResamplerFilter);
    3.83 +            ResamplerFilter = NULL;
    3.84 +            SDL_AtomicUnlock(&ResampleFilterSpinlock);
    3.85 +            return SDL_OutOfMemory();
    3.86 +        }
    3.87 +        kaiser_and_sinc(ResamplerFilter, ResamplerFilterDifference, RESAMPLER_FILTER_SIZE, beta);
    3.88 +    }
    3.89 +    SDL_AtomicUnlock(&ResampleFilterSpinlock);
    3.90 +    return 0;
    3.91 +}
    3.92 +
    3.93 +void
    3.94 +SDL_FreeResampleFilter(void)
    3.95 +{
    3.96 +    SDL_free(ResamplerFilter);
    3.97 +    SDL_free(ResamplerFilterDifference);
    3.98 +    ResamplerFilter = NULL;
    3.99 +    ResamplerFilterDifference = NULL;
   3.100 +}
   3.101 +
   3.102  
   3.103  static int
   3.104 -SDL_ResampleAudioSimple(const int chans, const double rate_incr,
   3.105 +SDL_ResampleAudio(const int chans, const int inrate, const int outrate,
   3.106                          float *last_sample, const float *inbuf,
   3.107                          const int inbuflen, float *outbuf, const int outbuflen)
   3.108  {
   3.109 +    const float outtimeincr = 1.0f / ((float) outrate);
   3.110 +    const float ratio = ((float) outrate) / ((float) inrate);
   3.111 +    /*const int padding_len = (ratio < 1.0f) ? (int) SDL_ceilf(((float) (RESAMPLER_SAMPLES_PER_ZERO_CROSSING * inrate) / ((float) outrate))) : RESAMPLER_SAMPLES_PER_ZERO_CROSSING;*/
   3.112      const int framelen = chans * (int)sizeof (float);
   3.113 -    const int total = (inbuflen / framelen);
   3.114 -    const int finalpos = (total * chans) - chans;
   3.115 -    const int dest_samples = (int)(((double)total) * rate_incr);
   3.116 -    const double src_incr = 1.0 / rate_incr;
   3.117 -    float *dst;
   3.118 -    double idx;
   3.119 -    int i;
   3.120 +    const int inframes = inbuflen / framelen;
   3.121 +    const int wantedoutframes = (int) ((inbuflen / framelen) * ratio);  /* outbuflen isn't total to write, it's total available. */
   3.122 +    const int maxoutframes = outbuflen / framelen;
   3.123 +    const int outframes = (wantedoutframes < maxoutframes) ? wantedoutframes : maxoutframes;
   3.124 +    float *dst = outbuf;
   3.125 +    float outtime = 0.0f;
   3.126 +    int i, j, chan;
   3.127  
   3.128 -    SDL_assert((dest_samples * framelen) <= outbuflen);
   3.129 -    SDL_assert((inbuflen % framelen) == 0);
   3.130 +    for (i = 0; i < outframes; i++) {
   3.131 +        const int srcindex = (int) (outtime * inrate);
   3.132 +        const float finrate = (float) inrate;
   3.133 +        const float intime = ((float) srcindex) / finrate;
   3.134 +        const float innexttime = ((float) (srcindex + 1)) / finrate;
   3.135  
   3.136 -    if (rate_incr > 1.0) {  /* upsample */
   3.137 -        float *target = (outbuf + chans);
   3.138 -        dst = outbuf + (dest_samples * chans);
   3.139 -        idx = (double) total;
   3.140 +        const float interpolation1 = 1.0f - (innexttime - outtime) / (innexttime - intime);
   3.141 +        const int filterindex1 = (int) (interpolation1 * RESAMPLER_SAMPLES_PER_ZERO_CROSSING);
   3.142 +        const float interpolation2 = 1.0f - interpolation1;
   3.143 +        const int filterindex2 = interpolation2 * RESAMPLER_SAMPLES_PER_ZERO_CROSSING;
   3.144  
   3.145 -        if (chans == 1) {
   3.146 -            const float final_sample = inbuf[finalpos];
   3.147 -            float earlier_sample = inbuf[finalpos];
   3.148 -            while (dst > target) {
   3.149 -                const int pos = ((int) idx) * chans;
   3.150 -                const float *src = &inbuf[pos];
   3.151 -                const float val = *(--src);
   3.152 -                SDL_assert(pos >= 0.0);
   3.153 -                *(--dst) = (val + earlier_sample) * 0.5f;
   3.154 -                earlier_sample = val;
   3.155 -                idx -= src_incr;
   3.156 +        for (chan = 0; chan < chans; chan++) {
   3.157 +            float outsample = 0.0f;
   3.158 +
   3.159 +            /* do this twice to calculate the sample, once for the "left wing" and then same for the right. */
   3.160 +            /* !!! FIXME: do both wings in one loop */
   3.161 +            for (j = 0; (filterindex1 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) < RESAMPLER_FILTER_SIZE; j++) {
   3.162 +                /* !!! FIXME: insample uses zero for padding samples, but it should use prior state from last_sample. */
   3.163 +                const int srcframe = srcindex - j;
   3.164 +                const float insample = (srcframe < 0) ? 0.0f : inbuf[(srcframe * chans) + chan];  /* !!! FIXME: we can bubble this conditional out of here by doing a pre loop. */
   3.165 +                outsample += (insample * (ResamplerFilter[filterindex1 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)] + (interpolation1 * ResamplerFilterDifference[filterindex1 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)])));
   3.166              }
   3.167 -            /* do last sample, interpolated against previous run's state. */
   3.168 -            *(--dst) = (inbuf[0] + last_sample[0]) * 0.5f;
   3.169 -            *last_sample = final_sample;
   3.170 -        } else if (chans == 2) {
   3.171 -            const float final_sample2 = inbuf[finalpos+1];
   3.172 -            const float final_sample1 = inbuf[finalpos];
   3.173 -            float earlier_sample2 = inbuf[finalpos];
   3.174 -            float earlier_sample1 = inbuf[finalpos-1];
   3.175 -            while (dst > target) {
   3.176 -                const int pos = ((int) idx) * chans;
   3.177 -                const float *src = &inbuf[pos];
   3.178 -                const float val2 = *(--src);
   3.179 -                const float val1 = *(--src);
   3.180 -                SDL_assert(pos >= 0.0);
   3.181 -                *(--dst) = (val2 + earlier_sample2) * 0.5f;
   3.182 -                *(--dst) = (val1 + earlier_sample1) * 0.5f;
   3.183 -                earlier_sample2 = val2;
   3.184 -                earlier_sample1 = val1;
   3.185 -                idx -= src_incr;
   3.186 +
   3.187 +            for (j = 0; (filterindex2 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)) < RESAMPLER_FILTER_SIZE; j++) {
   3.188 +                const int srcframe = srcindex + 1 + j;
   3.189 +                /* !!! FIXME: insample uses zero for padding samples, but it should use prior state from last_sample. */
   3.190 +                const float insample = (srcframe >= inframes) ? 0.0f : inbuf[(srcframe * chans) + chan];  /* !!! FIXME: we can bubble this conditional out of here by doing a post loop. */
   3.191 +                outsample += (insample * (ResamplerFilter[filterindex2 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)] + (interpolation2 * ResamplerFilterDifference[filterindex2 + (j * RESAMPLER_SAMPLES_PER_ZERO_CROSSING)])));
   3.192              }
   3.193 -            /* do last sample, interpolated against previous run's state. */
   3.194 -            *(--dst) = (inbuf[1] + last_sample[1]) * 0.5f;
   3.195 -            *(--dst) = (inbuf[0] + last_sample[0]) * 0.5f;
   3.196 -            last_sample[1] = final_sample2;
   3.197 -            last_sample[0] = final_sample1;
   3.198 -        } else {
   3.199 -            const float *earlier_sample = &inbuf[finalpos];
   3.200 -            float final_sample[8];
   3.201 -            SDL_memcpy(final_sample, &inbuf[finalpos], framelen);
   3.202 -            while (dst > target) {
   3.203 -                const int pos = ((int) idx) * chans;
   3.204 -                const float *src = &inbuf[pos];
   3.205 -                SDL_assert(pos >= 0.0);
   3.206 -                for (i = chans - 1; i >= 0; i--) {
   3.207 -                    const float val = *(--src);
   3.208 -                    *(--dst) = (val + earlier_sample[i]) * 0.5f;
   3.209 -                }
   3.210 -                earlier_sample = src;
   3.211 -                idx -= src_incr;
   3.212 -            }
   3.213 -            /* do last sample, interpolated against previous run's state. */
   3.214 -            for (i = chans - 1; i >= 0; i--) {
   3.215 -                const float val = inbuf[i];
   3.216 -                *(--dst) = (val + last_sample[i]) * 0.5f;
   3.217 -            }
   3.218 -            SDL_memcpy(last_sample, final_sample, framelen);
   3.219 +            *(dst++) = outsample;
   3.220          }
   3.221  
   3.222 -        dst = (outbuf + (dest_samples * chans));
   3.223 -    } else {  /* downsample */
   3.224 -        float *target = (outbuf + (dest_samples * chans));
   3.225 -        dst = outbuf;
   3.226 -        idx = 0.0;
   3.227 -        if (chans == 1) {
   3.228 -            float last = *last_sample;
   3.229 -            while (dst < target) {
   3.230 -                const int pos = ((int) idx) * chans;
   3.231 -                const float val = inbuf[pos];
   3.232 -                SDL_assert(pos <= finalpos);
   3.233 -                *(dst++) = (val + last) * 0.5f;
   3.234 -                last = val;
   3.235 -                idx += src_incr;
   3.236 -            }
   3.237 -            *last_sample = last;
   3.238 -        } else if (chans == 2) {
   3.239 -            float last1 = last_sample[0];
   3.240 -            float last2 = last_sample[1];
   3.241 -            while (dst < target) {
   3.242 -                const int pos = ((int) idx) * chans;
   3.243 -                const float val1 = inbuf[pos];
   3.244 -                const float val2 = inbuf[pos+1];
   3.245 -                SDL_assert(pos <= finalpos);
   3.246 -                *(dst++) = (val1 + last1) * 0.5f;
   3.247 -                *(dst++) = (val2 + last2) * 0.5f;
   3.248 -                last1 = val1;
   3.249 -                last2 = val2;
   3.250 -                idx += src_incr;
   3.251 -            }
   3.252 -            last_sample[0] = last1;
   3.253 -            last_sample[1] = last2;
   3.254 -        } else {
   3.255 -            while (dst < target) {
   3.256 -                const int pos = ((int) idx) * chans;
   3.257 -                const float *src = &inbuf[pos];
   3.258 -                SDL_assert(pos <= finalpos);
   3.259 -                for (i = 0; i < chans; i++) {
   3.260 -                    const float val = *(src++);
   3.261 -                    *(dst++) = (val + last_sample[i]) * 0.5f;
   3.262 -                    last_sample[i] = val;
   3.263 -                }
   3.264 -                idx += src_incr;
   3.265 -            }
   3.266 -        }
   3.267 +        outtime += outtimeincr;
   3.268      }
   3.269  
   3.270 -    return (int) ((dst - outbuf) * ((int) sizeof (float)));
   3.271 +    return outframes * chans * sizeof (float);
   3.272  }
   3.273  
   3.274 -/* We keep one special-case fast path around for an extremely common audio format. */
   3.275 -static int
   3.276 -SDL_ResampleAudioSimple_si16_c2(const double rate_incr,
   3.277 -                        Sint16 *last_sample, const Sint16 *inbuf,
   3.278 -                        const int inbuflen, Sint16 *outbuf, const int outbuflen)
   3.279 -{
   3.280 -    const int chans = 2;
   3.281 -    const int framelen = 4;  /* stereo 16 bit */
   3.282 -    const int total = (inbuflen / framelen);
   3.283 -    const int finalpos = (total * chans) - chans;
   3.284 -    const int dest_samples = (int)(((double)total) * rate_incr);
   3.285 -    const double src_incr = 1.0 / rate_incr;
   3.286 -    Sint16 *dst;
   3.287 -    double idx;
   3.288 -
   3.289 -    SDL_assert((dest_samples * framelen) <= outbuflen);
   3.290 -    SDL_assert((inbuflen % framelen) == 0);
   3.291 -
   3.292 -    if (rate_incr > 1.0) {
   3.293 -        Sint16 *target = (outbuf + chans);
   3.294 -        const Sint16 final_right = inbuf[finalpos+1];
   3.295 -        const Sint16 final_left = inbuf[finalpos];
   3.296 -        Sint16 earlier_right = inbuf[finalpos-1];
   3.297 -        Sint16 earlier_left = inbuf[finalpos-2];
   3.298 -        dst = outbuf + (dest_samples * chans);
   3.299 -        idx = (double) total;
   3.300 -
   3.301 -        while (dst > target) {
   3.302 -            const int pos = ((int) idx) * chans;
   3.303 -            const Sint16 *src = &inbuf[pos];
   3.304 -            const Sint16 right = *(--src);
   3.305 -            const Sint16 left = *(--src);
   3.306 -            SDL_assert(pos >= 0.0);
   3.307 -            *(--dst) = (((Sint32) right) + ((Sint32) earlier_right)) >> 1;
   3.308 -            *(--dst) = (((Sint32) left) + ((Sint32) earlier_left)) >> 1;
   3.309 -            earlier_right = right;
   3.310 -            earlier_left = left;
   3.311 -            idx -= src_incr;
   3.312 -        }
   3.313 -
   3.314 -        /* do last sample, interpolated against previous run's state. */
   3.315 -        *(--dst) = (((Sint32) inbuf[1]) + ((Sint32) last_sample[1])) >> 1;
   3.316 -        *(--dst) = (((Sint32) inbuf[0]) + ((Sint32) last_sample[0])) >> 1;
   3.317 -        last_sample[1] = final_right;
   3.318 -        last_sample[0] = final_left;
   3.319 -
   3.320 -        dst = (outbuf + (dest_samples * chans));
   3.321 -    } else {
   3.322 -        Sint16 *target = (outbuf + (dest_samples * chans));
   3.323 -        dst = outbuf;
   3.324 -        idx = 0.0;
   3.325 -        while (dst < target) {
   3.326 -            const int pos = ((int) idx) * chans;
   3.327 -            const Sint16 *src = &inbuf[pos];
   3.328 -            const Sint16 left = *(src++);
   3.329 -            const Sint16 right = *(src++);
   3.330 -            SDL_assert(pos <= finalpos);
   3.331 -            *(dst++) = (((Sint32) left) + ((Sint32) last_sample[0])) >> 1;
   3.332 -            *(dst++) = (((Sint32) right) + ((Sint32) last_sample[1])) >> 1;
   3.333 -            last_sample[0] = left;
   3.334 -            last_sample[1] = right;
   3.335 -            idx += src_incr;
   3.336 -        }
   3.337 -    }
   3.338 -
   3.339 -    return (int) ((dst - outbuf) * ((int) sizeof (Sint16)));
   3.340 -}
   3.341 -
   3.342 -static void SDLCALL
   3.343 -SDL_ResampleCVT_si16_c2(SDL_AudioCVT *cvt, SDL_AudioFormat format)
   3.344 -{
   3.345 -    const Sint16 *src = (const Sint16 *) cvt->buf;
   3.346 -    const int srclen = cvt->len_cvt;
   3.347 -    Sint16 *dst = (Sint16 *) cvt->buf;
   3.348 -    const int dstlen = (cvt->len * cvt->len_mult);
   3.349 -    Sint16 state[2];
   3.350 -
   3.351 -    state[0] = src[0];
   3.352 -    state[1] = src[1];
   3.353 -
   3.354 -    SDL_assert(format == AUDIO_S16SYS);
   3.355 -
   3.356 -    cvt->len_cvt = SDL_ResampleAudioSimple_si16_c2(cvt->rate_incr, state, src, srclen, dst, dstlen);
   3.357 -    if (cvt->filters[++cvt->filter_index]) {
   3.358 -        cvt->filters[cvt->filter_index](cvt, format);
   3.359 -    }
   3.360 -}
   3.361 -
   3.362 -
   3.363  int
   3.364  SDL_ConvertAudio(SDL_AudioCVT * cvt)
   3.365  {
   3.366 @@ -761,17 +690,28 @@
   3.367  static void
   3.368  SDL_ResampleCVT(SDL_AudioCVT *cvt, const int chans, const SDL_AudioFormat format)
   3.369  {
   3.370 +    /* !!! FIXME in 2.1: there are ten slots in the filter list, and the theoretical maximum we use is six (seven with NULL terminator).
   3.371 +       !!! FIXME in 2.1:   We need to store data for this resampler, because the cvt structure doesn't store the original sample rates,
   3.372 +       !!! FIXME in 2.1:   so we steal the ninth and tenth slot.  :( */
   3.373 +    const int srcrate = (int) (size_t) cvt->filters[SDL_AUDIOCVT_MAX_FILTERS-1];
   3.374 +    const int dstrate = (int) (size_t) cvt->filters[SDL_AUDIOCVT_MAX_FILTERS];
   3.375      const float *src = (const float *) cvt->buf;
   3.376      const int srclen = cvt->len_cvt;
   3.377 -    float *dst = (float *) cvt->buf;
   3.378 -    const int dstlen = (cvt->len * cvt->len_mult);
   3.379 +    /*float *dst = (float *) cvt->buf;
   3.380 +    const int dstlen = (cvt->len * cvt->len_mult);*/
   3.381 +    /* !!! FIXME: remove this if we can get the resampler to work in-place again. */
   3.382 +    float *dst = (float *) (cvt->buf + srclen);
   3.383 +    const int dstlen = (cvt->len * cvt->len_mult) - srclen;
   3.384      float state[8];
   3.385  
   3.386      SDL_assert(format == AUDIO_F32SYS);
   3.387  
   3.388 -    SDL_memcpy(state, src, chans*sizeof(*src));
   3.389 +    SDL_zero(state);
   3.390  
   3.391 -    cvt->len_cvt = SDL_ResampleAudioSimple(chans, cvt->rate_incr, state, src, srclen, dst, dstlen);
   3.392 +    cvt->len_cvt = SDL_ResampleAudio(chans, srcrate, dstrate, state, src, srclen, dst, dstlen);
   3.393 +
   3.394 +    SDL_memcpy(cvt->buf, dst, cvt->len_cvt);  /* !!! FIXME: remove this if we can get the resampler to work in-place again. */
   3.395 +
   3.396      if (cvt->filters[++cvt->filter_index]) {
   3.397          cvt->filters[cvt->filter_index](cvt, format);
   3.398      }
   3.399 @@ -823,10 +763,24 @@
   3.400          return SDL_SetError("No conversion available for these rates");
   3.401      }
   3.402  
   3.403 +    if (SDL_PrepareResampleFilter() < 0) {
   3.404 +        return -1;
   3.405 +    }
   3.406 +
   3.407      /* Update (cvt) with filter details... */
   3.408      if (SDL_AddAudioCVTFilter(cvt, filter) < 0) {
   3.409          return -1;
   3.410      }
   3.411 +
   3.412 +    /* !!! FIXME in 2.1: there are ten slots in the filter list, and the theoretical maximum we use is six (seven with NULL terminator).
   3.413 +       !!! FIXME in 2.1:   We need to store data for this resampler, because the cvt structure doesn't store the original sample rates,
   3.414 +       !!! FIXME in 2.1:   so we steal the ninth and tenth slot.  :( */
   3.415 +    if (cvt->filter_index >= (SDL_AUDIOCVT_MAX_FILTERS-2)) {
   3.416 +        return SDL_SetError("Too many filters needed for conversion, exceeded maximum of %d", SDL_AUDIOCVT_MAX_FILTERS-2);
   3.417 +    }
   3.418 +    cvt->filters[SDL_AUDIOCVT_MAX_FILTERS-1] = (SDL_AudioFilter) (size_t) src_rate;
   3.419 +    cvt->filters[SDL_AUDIOCVT_MAX_FILTERS] = (SDL_AudioFilter) (size_t) dst_rate;
   3.420 +
   3.421      if (src_rate < dst_rate) {
   3.422          const double mult = ((double) dst_rate) / ((double) src_rate);
   3.423          cvt->len_mult *= (int) SDL_ceil(mult);
   3.424 @@ -835,6 +789,11 @@
   3.425          cvt->len_ratio /= ((double) src_rate) / ((double) dst_rate);
   3.426      }
   3.427  
   3.428 +    /* !!! FIXME: remove this if we can get the resampler to work in-place again. */
   3.429 +    /* the buffer is big enough to hold the destination now, but
   3.430 +       we need it large enough to hold a separate scratch buffer. */
   3.431 +    cvt->len_mult *= 2;
   3.432 +
   3.433      return 1;               /* added a converter. */
   3.434  }
   3.435  
   3.436 @@ -922,7 +881,7 @@
   3.437      cvt->dst_format = dst_fmt;
   3.438      cvt->needed = 0;
   3.439      cvt->filter_index = 0;
   3.440 -    cvt->filters[0] = NULL;
   3.441 +    SDL_zero(cvt->filters);
   3.442      cvt->len_mult = 1;
   3.443      cvt->len_ratio = 1.0;
   3.444      cvt->rate_incr = ((double) dst_rate) / ((double) src_rate);
   3.445 @@ -930,32 +889,6 @@
   3.446      /* Make sure we've chosen audio conversion functions (MMX, scalar, etc.) */
   3.447      SDL_ChooseAudioConverters();
   3.448  
   3.449 -    /* SDL now favors float32 as its preferred internal format, and considers
   3.450 -       everything else to be a degenerate case that we might have to make
   3.451 -       multiple passes over the data to convert to and from float32 as
   3.452 -       necessary. That being said, we keep one special case around for
   3.453 -       efficiency: stereo data in Sint16 format, in the native byte order,
   3.454 -       that only needs resampling. This is likely to be the most popular
   3.455 -       legacy format, that apps, hardware and the OS are likely to be able
   3.456 -       to process directly, so we handle this one case directly without
   3.457 -       unnecessary conversions. This means that apps on embedded devices
   3.458 -       without floating point hardware should consider aiming for this
   3.459 -       format as well. */
   3.460 -    if ((src_channels == 2) && (dst_channels == 2) && (src_fmt == AUDIO_S16SYS) && (dst_fmt == AUDIO_S16SYS) && (src_rate != dst_rate)) {
   3.461 -        cvt->needed = 1;
   3.462 -        if (SDL_AddAudioCVTFilter(cvt, SDL_ResampleCVT_si16_c2) < 0) {
   3.463 -            return -1;
   3.464 -        }
   3.465 -        if (src_rate < dst_rate) {
   3.466 -            const double mult = ((double) dst_rate) / ((double) src_rate);
   3.467 -            cvt->len_mult *= (int) SDL_ceil(mult);
   3.468 -            cvt->len_ratio *= mult;
   3.469 -        } else {
   3.470 -            cvt->len_ratio /= ((double) src_rate) / ((double) dst_rate);
   3.471 -        }
   3.472 -        return 1;
   3.473 -    }
   3.474 -
   3.475      /* Type conversion goes like this now:
   3.476          - byteswap to CPU native format first if necessary.
   3.477          - convert to native Float32 if necessary.
   3.478 @@ -1282,30 +1215,23 @@
   3.479  
   3.480      SDL_assert(chans <= SDL_arraysize(state->resampler_state.f));
   3.481  
   3.482 +    if (inbuf == ((const float *) outbuf)) {  /* !!! FIXME can't work in-place (for now!). */
   3.483 +        Uint8 *ptr = EnsureStreamBufferSize(stream, inbuflen + outbuflen);
   3.484 +        if (ptr == NULL) {
   3.485 +            SDL_OutOfMemory();
   3.486 +            return 0;
   3.487 +        }
   3.488 +        SDL_memcpy(ptr + outbuflen, ptr, inbuflen);
   3.489 +        inbuf = (const float *) (ptr + outbuflen);
   3.490 +        outbuf = (float *) ptr;
   3.491 +    }
   3.492 +
   3.493      if (!state->resampler_seeded) {
   3.494 -        SDL_memcpy(state->resampler_state.f, inbuf, chans * sizeof (float));
   3.495 +        SDL_zero(state->resampler_state.f);
   3.496          state->resampler_seeded = SDL_TRUE;
   3.497      }
   3.498  
   3.499 -    return SDL_ResampleAudioSimple(chans, stream->rate_incr, state->resampler_state.f, inbuf, inbuflen, outbuf, outbuflen);
   3.500 -}
   3.501 -
   3.502 -static int
   3.503 -SDL_ResampleAudioStream_si16_c2(SDL_AudioStream *stream, const void *_inbuf, const int inbuflen, void *_outbuf, const int outbuflen)
   3.504 -{
   3.505 -    const Sint16 *inbuf = (const Sint16 *) _inbuf;
   3.506 -    Sint16 *outbuf = (Sint16 *) _outbuf;
   3.507 -    SDL_AudioStreamResamplerState *state = (SDL_AudioStreamResamplerState*)stream->resampler_state;
   3.508 -
   3.509 -    SDL_assert(((int)stream->pre_resample_channels) <= SDL_arraysize(state->resampler_state.si16));
   3.510 -
   3.511 -    if (!state->resampler_seeded) {
   3.512 -        state->resampler_state.si16[0] = inbuf[0];
   3.513 -        state->resampler_state.si16[1] = inbuf[1];
   3.514 -        state->resampler_seeded = SDL_TRUE;
   3.515 -    }
   3.516 -
   3.517 -    return SDL_ResampleAudioSimple_si16_c2(stream->rate_incr, state->resampler_state.si16, inbuf, inbuflen, outbuf, outbuflen);
   3.518 +    return SDL_ResampleAudio(chans, stream->src_rate, stream->dst_rate, state->resampler_state.f, inbuf, inbuflen, outbuf, outbuflen);
   3.519  }
   3.520  
   3.521  static void
   3.522 @@ -1332,9 +1258,6 @@
   3.523      const int packetlen = 4096;  /* !!! FIXME: good enough for now. */
   3.524      Uint8 pre_resample_channels;
   3.525      SDL_AudioStream *retval;
   3.526 -#ifndef HAVE_LIBSAMPLERATE_H
   3.527 -    const SDL_bool SRC_available = SDL_FALSE;
   3.528 -#endif
   3.529  
   3.530      retval = (SDL_AudioStream *) SDL_calloc(1, sizeof (SDL_AudioStream));
   3.531      if (!retval) {
   3.532 @@ -1366,18 +1289,6 @@
   3.533              SDL_FreeAudioStream(retval);
   3.534              return NULL;  /* SDL_BuildAudioCVT should have called SDL_SetError. */
   3.535          }
   3.536 -    /* fast path special case for stereo Sint16 data that just needs resampling. */
   3.537 -    } else if ((!SRC_available) && (src_channels == 2) && (dst_channels == 2) && (src_format == AUDIO_S16SYS) && (dst_format == AUDIO_S16SYS)) {
   3.538 -        SDL_assert(src_rate != dst_rate);
   3.539 -        retval->resampler_state = SDL_calloc(1, sizeof(SDL_AudioStreamResamplerState));
   3.540 -        if (!retval->resampler_state) {
   3.541 -            SDL_FreeAudioStream(retval);
   3.542 -            SDL_OutOfMemory();
   3.543 -            return NULL;
   3.544 -        }
   3.545 -        retval->resampler_func = SDL_ResampleAudioStream_si16_c2;
   3.546 -        retval->reset_resampler_func = SDL_ResetAudioStreamResampler;
   3.547 -        retval->cleanup_resampler_func = SDL_CleanupAudioStreamResampler;
   3.548      } else {
   3.549          /* Don't resample at first. Just get us to Float32 format. */
   3.550          /* !!! FIXME: convert to int32 on devices without hardware float. */
   3.551 @@ -1397,6 +1308,14 @@
   3.552                  SDL_OutOfMemory();
   3.553                  return NULL;
   3.554              }
   3.555 +
   3.556 +            if (SDL_PrepareResampleFilter() < 0) {
   3.557 +                SDL_free(retval->resampler_state);
   3.558 +                retval->resampler_state = NULL;
   3.559 +                SDL_FreeAudioStream(retval);
   3.560 +                return NULL;
   3.561 +            }
   3.562 +
   3.563              retval->resampler_func = SDL_ResampleAudioStream;
   3.564              retval->reset_resampler_func = SDL_ResetAudioStreamResampler;
   3.565              retval->cleanup_resampler_func = SDL_CleanupAudioStreamResampler;
     4.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.2 +++ b/src/audio/kaiser_window.pl	Thu Sep 21 02:51:14 2017 -0400
     4.3 @@ -0,0 +1,210 @@
     4.4 +#!/usr/bin/perl -w
     4.5 +
     4.6 +use warnings;
     4.7 +use strict;
     4.8 +
     4.9 +# The resampling algorithm: https://ccrma.stanford.edu/~jos/resample/
    4.10 +# https://www.mathworks.com/help/signal/ref/kaiser.html
    4.11 +# "Thus kaiser(L,beta) is equivalent to
    4.12 +#    besseli(0,beta*sqrt(1-(((0:L-1)-(L-1)/2)/((L-1)/2)).^2))/besseli(0,beta)."
    4.13 +# Matlab kaiser calls besseli():
    4.14 +# https://www.mathworks.com/help/matlab/ref/besseli.htm
    4.15 +# https://en.wikipedia.org/wiki/Bessel_function
    4.16 +
    4.17 +sub print_table {
    4.18 +    my $tableref = shift;
    4.19 +    my $name = shift;
    4.20 +    my @table = @{$tableref};
    4.21 +    my $comma = '';
    4.22 +    my $count = 0;
    4.23 +    print("static const float $name = {\n    ");
    4.24 +    foreach (@table) {
    4.25 +        print("$comma$_");
    4.26 +        #print(sprintf("%.6f\n", $_));
    4.27 +        if (++$count > 4) {
    4.28 +            $count = 0;
    4.29 +            print(",\n    ");
    4.30 +            $comma = '';
    4.31 +        } else {
    4.32 +            $comma = ', ';
    4.33 +        }
    4.34 +    }
    4.35 +    print("\n};\n\n");
    4.36 +}
    4.37 +
    4.38 +
    4.39 +use POSIX ();
    4.40 +
    4.41 +# This is a "modified" bessel function, so you can't use POSIX j0()
    4.42 +sub bessel {
    4.43 +    my $x = shift;
    4.44 +
    4.45 +    my $i0 = 1;
    4.46 +    my $f = 1;
    4.47 +    my $i = 1;
    4.48 +
    4.49 +    while (1) {
    4.50 +        my $diff = POSIX::pow($x / 2.0, $i * 2) / POSIX::pow($f, 2);
    4.51 +        last if ($diff < 1.0e-21);
    4.52 +        $i0 += $diff;
    4.53 +        $i++;
    4.54 +        $f *= $i;
    4.55 +    }
    4.56 +
    4.57 +    return $i0;
    4.58 +}
    4.59 +
    4.60 +sub kaiser {
    4.61 +    my $L = shift;
    4.62 +    my $beta = shift;
    4.63 +    my @retval;
    4.64 +
    4.65 +    #print("L=$L, beta=$beta\n"); exit(0);
    4.66 +
    4.67 +    for (my $i = 0; $i < $L; $i++) {
    4.68 +        my $val = bessel($beta * sqrt(1.0 - 
    4.69 +        POSIX::pow(
    4.70 +            (
    4.71 +                (
    4.72 +                    ($i-($L-1.0))
    4.73 +                ) / 2.0
    4.74 +            ) / (($L-1)/2.0), 2.0 ))
    4.75 +        ) / bessel($beta);
    4.76 +
    4.77 +        unshift @retval, $val;
    4.78 +    }
    4.79 +    return @retval;
    4.80 +}
    4.81 +
    4.82 +
    4.83 +my $zero_crossings = 5;
    4.84 +my $bits_per_sample = 16;
    4.85 +my $samples_per_zero_crossing = 1 << (($bits_per_sample / 2) + 1);
    4.86 +my $kaiser_window_table_size = ($samples_per_zero_crossing * $zero_crossings) + 1;
    4.87 +
    4.88 +# if dB > 50: 0.1102 * ($db - 8.7)
    4.89 +my $db = 80.0;
    4.90 +my $beta = 0.1102 * ($db - 8.7);
    4.91 +
    4.92 +my @table = kaiser($kaiser_window_table_size, $beta);
    4.93 +
    4.94 +print_table(\@table, 'kaiser_window');
    4.95 +
    4.96 +# Kaiser window has "sinc function" ("cardinal sine") applied to it:
    4.97 +#   sin(pi * x) / (pi * x)
    4.98 +# "For example, to use the ideal lowpass filter, the table would contain
    4.99 +#  h(l) = sinc(l/L)."
   4.100 +
   4.101 +use Math::Trig ':pi';
   4.102 +for (my $i = 1; $i < $kaiser_window_table_size; $i++) {
   4.103 +    my $x = $i / $samples_per_zero_crossing;
   4.104 +    $table[$i] *= sin($x * pi) / ($x * pi);
   4.105 +}
   4.106 +
   4.107 +print_table(\@table, 'with_sinc');
   4.108 +
   4.109 +# "Our implementation also stores a table of differences ¯h(l) = h(l + 1) − h(l) between successive
   4.110 +# FIR sample values in order to speed up the linear interpolation. The length of each table is
   4.111 +# Nh = LNz + 1, including the endpoint definition ¯h(Nh) = 0."
   4.112 +
   4.113 +my @differences = ();
   4.114 +for (my $i = 1; $i < $kaiser_window_table_size; $i++) {
   4.115 +    push @differences, $table[$i] - $table[$i - 1];
   4.116 +}
   4.117 +push @differences, 0;
   4.118 +
   4.119 +print_table(\@differences, 'differences');
   4.120 +
   4.121 +
   4.122 +# Might as well use this code as a test harness...
   4.123 +
   4.124 +use autodie;
   4.125 +my $fnamein = shift @ARGV;
   4.126 +my $fnameout = shift @ARGV;
   4.127 +my $inrate = shift @ARGV;
   4.128 +my $outrate = shift @ARGV;
   4.129 +
   4.130 +print("Resampling $fnamein (freq=$inrate) to $fnameout (freq=$outrate).\n");
   4.131 +
   4.132 +open(IN, '<:raw', $fnamein);
   4.133 +my @src = ();
   4.134 +
   4.135 +# this assumes mono Sint16 raw data since we aren't parsing .wav files.
   4.136 +# !!! FIXME: deal with multichannel audio.
   4.137 +my $channels = 1;
   4.138 +
   4.139 +# this is inefficient, but this is just throwaway code...
   4.140 +while (read(IN, my $bytes, 2) == 2) {
   4.141 +    my ($samp) = unpack('s', $bytes);
   4.142 +    push @src, $samp;
   4.143 +}
   4.144 +
   4.145 +close(IN);
   4.146 +
   4.147 +my $ratio = $outrate / $inrate;
   4.148 +my $sample_frames_in = scalar(@src) / $channels;
   4.149 +my $sample_frames_out = $sample_frames_in * $ratio;
   4.150 +
   4.151 +my $outsamples = $sample_frames_out * $channels;
   4.152 +#my @dst = (0) x ($outsamples);
   4.153 +my @dst = ();
   4.154 +print("Resampling $sample_frames_in input frames to $sample_frames_out output (ratio=$ratio).\n");
   4.155 +
   4.156 +
   4.157 +my $inv_spzc = int(POSIX::ceil(($samples_per_zero_crossing * $inrate) / $outrate));
   4.158 +my $padding_len;
   4.159 +if ($ratio < 1.0) {
   4.160 +    $padding_len = int(POSIX::ceil(($samples_per_zero_crossing * $inrate) / $outrate));
   4.161 +} else {
   4.162 +    $padding_len = $samples_per_zero_crossing;
   4.163 +}
   4.164 +
   4.165 +# You need to pad the input or we'll get buffer overflows.
   4.166 +# !!! FIXME: deal with multichannel audio.
   4.167 +for (my $i = 0; $i < $padding_len; $i++) {
   4.168 +    push @src, 0;
   4.169 +    unshift @src, 0;
   4.170 +}
   4.171 +
   4.172 +# !!! FIXME: deal with multichannel audio.
   4.173 +my $time = 0.0;
   4.174 +for (my $i = 0; $i < $outsamples; $i++) {
   4.175 +    my $srcindex = int($time * $inrate);  # !!! FIXME: truncate or round?
   4.176 +
   4.177 +    my $ftime = $srcindex / $inrate;  # this would be $time if we didn't convert $srcindex to int.
   4.178 +    my $fnexttime = ($srcindex + 1) / $inrate;
   4.179 +
   4.180 +    # do this twice to calculate the sample, once for the "left wing" and then same for the right.
   4.181 +    my $sample = 0;
   4.182 +    my $interpolation = 1.0 - ($fnexttime - $time) / ($fnexttime - $ftime);
   4.183 +    my $filterindex = int($interpolation * $samples_per_zero_crossing);
   4.184 +
   4.185 +    $srcindex += $padding_len;
   4.186 +
   4.187 +    for (my $j = 0; ($filterindex + ($j * $samples_per_zero_crossing)) < $kaiser_window_table_size; $j++) {
   4.188 +        $sample += int($src[$srcindex - $j] * ($table[$filterindex + $j * $samples_per_zero_crossing] + $interpolation * $differences[$filterindex + $j * $samples_per_zero_crossing]));
   4.189 +    }
   4.190 +
   4.191 +    $interpolation = 1 - $interpolation;
   4.192 +    $filterindex = $interpolation * $samples_per_zero_crossing;
   4.193 +    for (my $j = 0; ($filterindex + ($j * $samples_per_zero_crossing)) < $kaiser_window_table_size; $j++) {
   4.194 +        $sample += int($src[$srcindex + 1 + $j] * ($table[$filterindex + $j * $samples_per_zero_crossing] + $interpolation * $differences[$filterindex + $j * $samples_per_zero_crossing]));
   4.195 +    }
   4.196 +
   4.197 +    push @dst, $sample;
   4.198 +
   4.199 +    # "After each output sample is computed, the time register is incremented by 2nl+nη /ρ (i.e., time is incremented by 1/ρ in fixed-point format)."
   4.200 +    $time += 1.0 / $outrate;
   4.201 +}
   4.202 +
   4.203 +open(OUT, '>:raw', $fnameout);
   4.204 +
   4.205 +# this is inefficient, but this is just throwaway code...
   4.206 +foreach (@dst) {
   4.207 +    print OUT pack('s', $_);
   4.208 +}
   4.209 +
   4.210 +close(OUT);
   4.211 +
   4.212 +print("Done.\n");
   4.213 +