git.sesse.net Git - ffmpeg/blob - libswresample/resample.c

   1 /*
   2  * audio resampling
   3  * Copyright (c) 2004-2012 Michael Niedermayer <michaelni@gmx.at>
   4  * bessel function: Copyright (c) 2006 Xiaogang Zhang
   5  *
   6  * This file is part of FFmpeg.
   7  *
   8  * FFmpeg is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU Lesser General Public
  10  * License as published by the Free Software Foundation; either
  11  * version 2.1 of the License, or (at your option) any later version.
  12  *
  13  * FFmpeg is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16  * Lesser General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU Lesser General Public
  19  * License along with FFmpeg; if not, write to the Free Software
  20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  21  */
  22
  23 /**
  24  * @file
  25  * audio resampling
  26  * @author Michael Niedermayer <michaelni@gmx.at>
  27  */
  28
  29 #include "libavutil/avassert.h"
  30 #include "resample.h"
  31
  32 static inline double eval_poly(const double *coeff, int size, double x) {
  33     double sum = coeff[size-1];
  34     int i;
  35     for (i = size-2; i >= 0; --i) {
  36         sum *= x;
  37         sum += coeff[i];
  38     }
  39     return sum;
  40 }
  41
  42 /**
  43  * 0th order modified bessel function of the first kind.
  44  * Algorithm taken from the Boost project, source:
  45  * https://searchcode.com/codesearch/view/14918379/
  46  * Use, modification and distribution are subject to the
  47  * Boost Software License, Version 1.0 (see notice below).
  48  * Boost Software License - Version 1.0 - August 17th, 2003
  49 Permission is hereby granted, free of charge, to any person or organization
  50 obtaining a copy of the software and accompanying documentation covered by
  51 this license (the "Software") to use, reproduce, display, distribute,
  52 execute, and transmit the Software, and to prepare derivative works of the
  53 Software, and to permit third-parties to whom the Software is furnished to
  54 do so, all subject to the following:
  55
  56 The copyright notices in the Software and this entire statement, including
  57 the above license grant, this restriction and the following disclaimer,
  58 must be included in all copies of the Software, in whole or in part, and
  59 all derivative works of the Software, unless such copies or derivative
  60 works are solely in the form of machine-executable object code generated by
  61 a source language processor.
  62
  63 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  64 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  65 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
  66 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
  67 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
  68 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  69 DEALINGS IN THE SOFTWARE.
  70  */
  71
  72 static double bessel(double x) {
  73 // Modified Bessel function of the first kind of order zero
  74 // minimax rational approximations on intervals, see
  75 // Blair and Edwards, Chalk River Report AECL-4928, 1974
  76     static const double p1[] = {
  77         -2.2335582639474375249e+15,
  78         -5.5050369673018427753e+14,
  79         -3.2940087627407749166e+13,
  80         -8.4925101247114157499e+11,
  81         -1.1912746104985237192e+10,
  82         -1.0313066708737980747e+08,
  83         -5.9545626019847898221e+05,
  84         -2.4125195876041896775e+03,
  85         -7.0935347449210549190e+00,
  86         -1.5453977791786851041e-02,
  87         -2.5172644670688975051e-05,
  88         -3.0517226450451067446e-08,
  89         -2.6843448573468483278e-11,
  90         -1.5982226675653184646e-14,
  91         -5.2487866627945699800e-18,
  92     };
  93     static const double q1[] = {
  94         -2.2335582639474375245e+15,
  95          7.8858692566751002988e+12,
  96         -1.2207067397808979846e+10,
  97          1.0377081058062166144e+07,
  98         -4.8527560179962773045e+03,
  99          1.0,
 100     };
 101     static const double p2[] = {
 102         -2.2210262233306573296e-04,
 103          1.3067392038106924055e-02,
 104         -4.4700805721174453923e-01,
 105          5.5674518371240761397e+00,
 106         -2.3517945679239481621e+01,
 107          3.1611322818701131207e+01,
 108         -9.6090021968656180000e+00,
 109     };
 110     static const double q2[] = {
 111         -5.5194330231005480228e-04,
 112          3.2547697594819615062e-02,
 113         -1.1151759188741312645e+00,
 114          1.3982595353892851542e+01,
 115         -6.0228002066743340583e+01,
 116          8.5539563258012929600e+01,
 117         -3.1446690275135491500e+01,
 118         1.0,
 119     };
 120     double y, r, factor;
 121     if (x == 0)
 122         return 1.0;
 123     x = fabs(x);
 124     if (x <= 15) {
 125         y = x * x;
 126         return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
 127     }
 128     else {
 129         y = 1 / x - 1.0 / 15;
 130         r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
 131         factor = exp(x) / sqrt(x);
 132         return factor * r;
 133     }
 134 }
 135
 136 /**
 137  * builds a polyphase filterbank.
 138  * @param factor resampling factor
 139  * @param scale wanted sum of coefficients for each filter
 140  * @param filter_type  filter type
 141  * @param kaiser_beta  kaiser window beta
 142  * @return 0 on success, negative on error
 143  */
 144 static int build_filter(ResampleContext *c, void *filter, double factor, int tap_count, int alloc, int phase_count, int scale,
 145                         int filter_type, double kaiser_beta){
 146     int ph, i;
 147     int ph_nb = phase_count % 2 ? phase_count : phase_count / 2 + 1;
 148     double x, y, w, t, s;
 149     double *tab = av_malloc_array(tap_count+1,  sizeof(*tab));
 150     double *sin_lut = av_malloc_array(ph_nb, sizeof(*sin_lut));
 151     const int center= (tap_count-1)/2;
 152     int ret = AVERROR(ENOMEM);
 153
 154     if (!tab || !sin_lut)
 155         goto fail;
 156
 157     /* if upsampling, only need to interpolate, no filter */
 158     if (factor > 1.0)
 159         factor = 1.0;
 160
 161     if (factor == 1.0) {
 162         for (ph = 0; ph < ph_nb; ph++)
 163             sin_lut[ph] = sin(M_PI * ph / phase_count);
 164     }
 165     for(ph = 0; ph < ph_nb; ph++) {
 166         double norm = 0;
 167         s = sin_lut[ph];
 168         for(i=0;i<=tap_count;i++) {
 169             x = M_PI * ((double)(i - center) - (double)ph / phase_count) * factor;
 170             if (x == 0) y = 1.0;
 171             else if (factor == 1.0)
 172                 y = s / x;
 173             else
 174                 y = sin(x) / x;
 175             switch(filter_type){
 176             case SWR_FILTER_TYPE_CUBIC:{
 177                 const float d= -0.5; //first order derivative = -0.5
 178                 x = fabs(((double)(i - center) - (double)ph / phase_count) * factor);
 179                 if(x<1.0) y= 1 - 3*x*x + 2*x*x*x + d*(            -x*x + x*x*x);
 180                 else      y=                       d*(-4 + 8*x - 5*x*x + x*x*x);
 181                 break;}
 182             case SWR_FILTER_TYPE_BLACKMAN_NUTTALL:
 183                 w = 2.0*x / (factor*tap_count);
 184                 t = -cos(w);
 185                 y *= 0.3635819 - 0.4891775 * t + 0.1365995 * (2*t*t-1) - 0.0106411 * (4*t*t*t - 3*t);
 186                 break;
 187             case SWR_FILTER_TYPE_KAISER:
 188                 w = 2.0*x / (factor*tap_count*M_PI);
 189                 y *= bessel(kaiser_beta*sqrt(FFMAX(1-w*w, 0)));
 190                 break;
 191             default:
 192                 av_assert0(0);
 193             }
 194
 195             tab[i] = y;
 196             s = -s;
 197             if (i < tap_count)
 198                 norm += y;
 199         }
 200
 201         /* normalize so that an uniform color remains the same */
 202         switch(c->format){
 203         case AV_SAMPLE_FMT_S16P:
 204             for(i=0;i<tap_count;i++)
 205                 ((int16_t*)filter)[ph * alloc + i] = av_clip_int16(lrintf(tab[i] * scale / norm));
 206             if (phase_count % 2) break;
 207             if (tap_count % 2 == 0 || tap_count == 1) {
 208                 for (i = 0; i < tap_count; i++)
 209                     ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int16_t*)filter)[ph * alloc + i];
 210             }
 211             else {
 212                 for (i = 1; i <= tap_count; i++)
 213                     ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
 214                         av_clip_int16(lrintf(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
 215             }
 216             break;
 217         case AV_SAMPLE_FMT_S32P:
 218             for(i=0;i<tap_count;i++)
 219                 ((int32_t*)filter)[ph * alloc + i] = av_clipl_int32(llrint(tab[i] * scale / norm));
 220             if (phase_count % 2) break;
 221             if (tap_count % 2 == 0 || tap_count == 1) {
 222                 for (i = 0; i < tap_count; i++)
 223                     ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int32_t*)filter)[ph * alloc + i];
 224             }
 225             else {
 226                 for (i = 1; i <= tap_count; i++)
 227                     ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
 228                         av_clipl_int32(llrint(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
 229             }
 230             break;
 231         case AV_SAMPLE_FMT_FLTP:
 232             for(i=0;i<tap_count;i++)
 233                 ((float*)filter)[ph * alloc + i] = tab[i] * scale / norm;
 234             if (phase_count % 2) break;
 235             if (tap_count % 2 == 0 || tap_count == 1) {
 236                 for (i = 0; i < tap_count; i++)
 237                     ((float*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((float*)filter)[ph * alloc + i];
 238             }
 239             else {
 240                 for (i = 1; i <= tap_count; i++)
 241                     ((float*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
 242             }
 243             break;
 244         case AV_SAMPLE_FMT_DBLP:
 245             for(i=0;i<tap_count;i++)
 246                 ((double*)filter)[ph * alloc + i] = tab[i] * scale / norm;
 247             if (phase_count % 2) break;
 248             if (tap_count % 2 == 0 || tap_count == 1) {
 249                 for (i = 0; i < tap_count; i++)
 250                     ((double*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((double*)filter)[ph * alloc + i];
 251             }
 252             else {
 253                 for (i = 1; i <= tap_count; i++)
 254                     ((double*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
 255             }
 256             break;
 257         }
 258     }
 259 #if 0
 260     {
 261 #define LEN 1024
 262         int j,k;
 263         double sine[LEN + tap_count];
 264         double filtered[LEN];
 265         double maxff=-2, minff=2, maxsf=-2, minsf=2;
 266         for(i=0; i<LEN; i++){
 267             double ss=0, sf=0, ff=0;
 268             for(j=0; j<LEN+tap_count; j++)
 269                 sine[j]= cos(i*j*M_PI/LEN);
 270             for(j=0; j<LEN; j++){
 271                 double sum=0;
 272                 ph=0;
 273                 for(k=0; k<tap_count; k++)
 274                     sum += filter[ph * tap_count + k] * sine[k+j];
 275                 filtered[j]= sum / (1<<FILTER_SHIFT);
 276                 ss+= sine[j + center] * sine[j + center];
 277                 ff+= filtered[j] * filtered[j];
 278                 sf+= sine[j + center] * filtered[j];
 279             }
 280             ss= sqrt(2*ss/LEN);
 281             ff= sqrt(2*ff/LEN);
 282             sf= 2*sf/LEN;
 283             maxff= FFMAX(maxff, ff);
 284             minff= FFMIN(minff, ff);
 285             maxsf= FFMAX(maxsf, sf);
 286             minsf= FFMIN(minsf, sf);
 287             if(i%11==0){
 288                 av_log(NULL, AV_LOG_ERROR, "i:%4d ss:%f ff:%13.6e-%13.6e sf:%13.6e-%13.6e\n", i, ss, maxff, minff, maxsf, minsf);
 289                 minff=minsf= 2;
 290                 maxff=maxsf= -2;
 291             }
 292         }
 293     }
 294 #endif
 295
 296     ret = 0;
 297 fail:
 298     av_free(tab);
 299     av_free(sin_lut);
 300     return ret;
 301 }
 302
 303 static ResampleContext *resample_init(ResampleContext *c, int out_rate, int in_rate, int filter_size, int phase_shift, int linear,
 304                                     double cutoff0, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta,
 305                                     double precision, int cheby, int exact_rational)
 306 {
 307     double cutoff = cutoff0? cutoff0 : 0.97;
 308     double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
 309     int phase_count= 1<<phase_shift;
 310     int phase_count_compensation = phase_count;
 311
 312     if (exact_rational) {
 313         int phase_count_exact, phase_count_exact_den;
 314
 315         av_reduce(&phase_count_exact, &phase_count_exact_den, out_rate, in_rate, INT_MAX);
 316         if (phase_count_exact <= phase_count) {
 317             phase_count_compensation = phase_count_exact * (phase_count / phase_count_exact);
 318             phase_count = phase_count_exact;
 319         }
 320     }
 321
 322     if (!c || c->phase_count != phase_count || c->linear!=linear || c->factor != factor
 323            || c->filter_length != FFMAX((int)ceil(filter_size/factor), 1) || c->format != format
 324            || c->filter_type != filter_type || c->kaiser_beta != kaiser_beta) {
 325         c = av_mallocz(sizeof(*c));
 326         if (!c)
 327             return NULL;
 328
 329         c->format= format;
 330
 331         c->felem_size= av_get_bytes_per_sample(c->format);
 332
 333         switch(c->format){
 334         case AV_SAMPLE_FMT_S16P:
 335             c->filter_shift = 15;
 336             break;
 337         case AV_SAMPLE_FMT_S32P:
 338             c->filter_shift = 30;
 339             break;
 340         case AV_SAMPLE_FMT_FLTP:
 341         case AV_SAMPLE_FMT_DBLP:
 342             c->filter_shift = 0;
 343             break;
 344         default:
 345             av_log(NULL, AV_LOG_ERROR, "Unsupported sample format\n");
 346             av_assert0(0);
 347         }
 348
 349         if (filter_size/factor > INT32_MAX/256) {
 350             av_log(NULL, AV_LOG_ERROR, "Filter length too large\n");
 351             goto error;
 352         }
 353
 354         c->phase_count   = phase_count;
 355         c->linear        = linear;
 356         c->factor        = factor;
 357         c->filter_length = FFMAX((int)ceil(filter_size/factor), 1);
 358         c->filter_alloc  = FFALIGN(c->filter_length, 8);
 359         c->filter_bank   = av_calloc(c->filter_alloc, (phase_count+1)*c->felem_size);
 360         c->filter_type   = filter_type;
 361         c->kaiser_beta   = kaiser_beta;
 362         c->phase_count_compensation = phase_count_compensation;
 363         if (!c->filter_bank)
 364             goto error;
 365         if (build_filter(c, (void*)c->filter_bank, factor, c->filter_length, c->filter_alloc, phase_count, 1<<c->filter_shift, filter_type, kaiser_beta))
 366             goto error;
 367         memcpy(c->filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, c->filter_bank, (c->filter_alloc-1)*c->felem_size);
 368         memcpy(c->filter_bank + (c->filter_alloc*phase_count  )*c->felem_size, c->filter_bank + (c->filter_alloc - 1)*c->felem_size, c->felem_size);
 369     }
 370
 371     c->compensation_distance= 0;
 372     if(!av_reduce(&c->src_incr, &c->dst_incr, out_rate, in_rate * (int64_t)phase_count, INT32_MAX/2))
 373         goto error;
 374     while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
 375         c->dst_incr *= 2;
 376         c->src_incr *= 2;
 377     }
 378     c->ideal_dst_incr = c->dst_incr;
 379     c->dst_incr_div   = c->dst_incr / c->src_incr;
 380     c->dst_incr_mod   = c->dst_incr % c->src_incr;
 381
 382     c->index= -phase_count*((c->filter_length-1)/2);
 383     c->frac= 0;
 384
 385     swri_resample_dsp_init(c);
 386
 387     return c;
 388 error:
 389     av_freep(&c->filter_bank);
 390     av_free(c);
 391     return NULL;
 392 }
 393
 394 static void resample_free(ResampleContext **c){
 395     if(!*c)
 396         return;
 397     av_freep(&(*c)->filter_bank);
 398     av_freep(c);
 399 }
 400
 401 static int rebuild_filter_bank_with_compensation(ResampleContext *c)
 402 {
 403     uint8_t *new_filter_bank;
 404     int new_src_incr, new_dst_incr;
 405     int phase_count = c->phase_count_compensation;
 406     int ret;
 407
 408     if (phase_count == c->phase_count)
 409         return 0;
 410
 411     av_assert0(!c->frac && !c->dst_incr_mod && !c->compensation_distance);
 412
 413     new_filter_bank = av_calloc(c->filter_alloc, (phase_count + 1) * c->felem_size);
 414     if (!new_filter_bank)
 415         return AVERROR(ENOMEM);
 416
 417     ret = build_filter(c, new_filter_bank, c->factor, c->filter_length, c->filter_alloc,
 418                        phase_count, 1 << c->filter_shift, c->filter_type, c->kaiser_beta);
 419     if (ret < 0) {
 420         av_freep(&new_filter_bank);
 421         return ret;
 422     }
 423     memcpy(new_filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, new_filter_bank, (c->filter_alloc-1)*c->felem_size);
 424     memcpy(new_filter_bank + (c->filter_alloc*phase_count  )*c->felem_size, new_filter_bank + (c->filter_alloc - 1)*c->felem_size, c->felem_size);
 425
 426     if (!av_reduce(&new_src_incr, &new_dst_incr, c->src_incr,
 427                    c->dst_incr * (int64_t)(phase_count/c->phase_count), INT32_MAX/2))
 428     {
 429         av_freep(&new_filter_bank);
 430         return AVERROR(EINVAL);
 431     }
 432
 433     c->src_incr = new_src_incr;
 434     c->dst_incr = new_dst_incr;
 435     while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
 436         c->dst_incr *= 2;
 437         c->src_incr *= 2;
 438     }
 439     c->ideal_dst_incr = c->dst_incr;
 440     c->dst_incr_div   = c->dst_incr / c->src_incr;
 441     c->dst_incr_mod   = c->dst_incr % c->src_incr;
 442     c->index         *= phase_count / c->phase_count;
 443     c->phase_count    = phase_count;
 444     av_freep(&c->filter_bank);
 445     c->filter_bank = new_filter_bank;
 446     return 0;
 447 }
 448
 449 static int set_compensation(ResampleContext *c, int sample_delta, int compensation_distance){
 450     int ret;
 451
 452     if (compensation_distance && sample_delta) {
 453         ret = rebuild_filter_bank_with_compensation(c);
 454         if (ret < 0)
 455             return ret;
 456     }
 457
 458     c->compensation_distance= compensation_distance;
 459     if (compensation_distance)
 460         c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
 461     else
 462         c->dst_incr = c->ideal_dst_incr;
 463
 464     c->dst_incr_div   = c->dst_incr / c->src_incr;
 465     c->dst_incr_mod   = c->dst_incr % c->src_incr;
 466
 467     return 0;
 468 }
 469
 470 static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, AudioData *src, int src_size, int *consumed){
 471     int i;
 472     int av_unused mm_flags = av_get_cpu_flags();
 473     int need_emms = c->format == AV_SAMPLE_FMT_S16P && ARCH_X86_32 &&
 474                     (mm_flags & (AV_CPU_FLAG_MMX2 | AV_CPU_FLAG_SSE2)) == AV_CPU_FLAG_MMX2;
 475     int64_t max_src_size = (INT64_MAX/2 / c->phase_count) / c->src_incr;
 476
 477     if (c->compensation_distance)
 478         dst_size = FFMIN(dst_size, c->compensation_distance);
 479     src_size = FFMIN(src_size, max_src_size);
 480
 481     *consumed = 0;
 482
 483     if (c->filter_length == 1 && c->phase_count == 1) {
 484         int64_t index2= (1LL<<32)*c->frac/c->src_incr + (1LL<<32)*c->index;
 485         int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
 486         int new_size = (src_size * (int64_t)c->src_incr - c->frac + c->dst_incr - 1) / c->dst_incr;
 487
 488         dst_size = FFMAX(FFMIN(dst_size, new_size), 0);
 489         if (dst_size > 0) {
 490             for (i = 0; i < dst->ch_count; i++) {
 491                 c->dsp.resample_one(dst->ch[i], src->ch[i], dst_size, index2, incr);
 492                 if (i+1 == dst->ch_count) {
 493                     c->index += dst_size * c->dst_incr_div;
 494                     c->index += (c->frac + dst_size * (int64_t)c->dst_incr_mod) / c->src_incr;
 495                     av_assert2(c->index >= 0);
 496                     *consumed = c->index;
 497                     c->frac   = (c->frac + dst_size * (int64_t)c->dst_incr_mod) % c->src_incr;
 498                     c->index = 0;
 499                 }
 500             }
 501         }
 502     } else {
 503         int64_t end_index = (1LL + src_size - c->filter_length) * c->phase_count;
 504         int64_t delta_frac = (end_index - c->index) * c->src_incr - c->frac;
 505         int delta_n = (delta_frac + c->dst_incr - 1) / c->dst_incr;
 506         int (*resample_func)(struct ResampleContext *c, void *dst,
 507                              const void *src, int n, int update_ctx);
 508
 509         dst_size = FFMAX(FFMIN(dst_size, delta_n), 0);
 510         if (dst_size > 0) {
 511             /* resample_linear and resample_common should have same behavior
 512              * when frac and dst_incr_mod are zero */
 513             resample_func = (c->linear && (c->frac || c->dst_incr_mod)) ?
 514                             c->dsp.resample_linear : c->dsp.resample_common;
 515             for (i = 0; i < dst->ch_count; i++)
 516                 *consumed = resample_func(c, dst->ch[i], src->ch[i], dst_size, i+1 == dst->ch_count);
 517         }
 518     }
 519
 520     if(need_emms)
 521         emms_c();
 522
 523     if (c->compensation_distance) {
 524         c->compensation_distance -= dst_size;
 525         if (!c->compensation_distance) {
 526             c->dst_incr     = c->ideal_dst_incr;
 527             c->dst_incr_div = c->dst_incr / c->src_incr;
 528             c->dst_incr_mod = c->dst_incr % c->src_incr;
 529         }
 530     }
 531
 532     return dst_size;
 533 }
 534
 535 static int64_t get_delay(struct SwrContext *s, int64_t base){
 536     ResampleContext *c = s->resample;
 537     int64_t num = s->in_buffer_count - (c->filter_length-1)/2;
 538     num *= c->phase_count;
 539     num -= c->index;
 540     num *= c->src_incr;
 541     num -= c->frac;
 542     return av_rescale(num, base, s->in_sample_rate*(int64_t)c->src_incr * c->phase_count);
 543 }
 544
 545 static int64_t get_out_samples(struct SwrContext *s, int in_samples) {
 546     ResampleContext *c = s->resample;
 547     // The + 2 are added to allow implementations to be slightly inaccurate, they should not be needed currently.
 548     // They also make it easier to proof that changes and optimizations do not
 549     // break the upper bound.
 550     int64_t num = s->in_buffer_count + 2LL + in_samples;
 551     num *= c->phase_count;
 552     num -= c->index;
 553     num = av_rescale_rnd(num, s->out_sample_rate, ((int64_t)s->in_sample_rate) * c->phase_count, AV_ROUND_UP) + 2;
 554
 555     if (c->compensation_distance) {
 556         if (num > INT_MAX)
 557             return AVERROR(EINVAL);
 558
 559         num = FFMAX(num, (num * c->ideal_dst_incr - 1) / c->dst_incr + 1);
 560     }
 561     return num;
 562 }
 563
 564 static int resample_flush(struct SwrContext *s) {
 565     AudioData *a= &s->in_buffer;
 566     int i, j, ret;
 567     if((ret = swri_realloc_audio(a, s->in_buffer_index + 2*s->in_buffer_count)) < 0)
 568         return ret;
 569     av_assert0(a->planar);
 570     for(i=0; i<a->ch_count; i++){
 571         for(j=0; j<s->in_buffer_count; j++){
 572             memcpy(a->ch[i] + (s->in_buffer_index+s->in_buffer_count+j  )*a->bps,
 573                 a->ch[i] + (s->in_buffer_index+s->in_buffer_count-j-1)*a->bps, a->bps);
 574         }
 575     }
 576     s->in_buffer_count += (s->in_buffer_count+1)/2;
 577     return 0;
 578 }
 579
 580 // in fact the whole handle multiple ridiculously small buffers might need more thinking...
 581 static int invert_initial_buffer(ResampleContext *c, AudioData *dst, const AudioData *src,
 582                                  int in_count, int *out_idx, int *out_sz)
 583 {
 584     int n, ch, num = FFMIN(in_count + *out_sz, c->filter_length + 1), res;
 585
 586     if (c->index >= 0)
 587         return 0;
 588
 589     if ((res = swri_realloc_audio(dst, c->filter_length * 2 + 1)) < 0)
 590         return res;
 591
 592     // copy
 593     for (n = *out_sz; n < num; n++) {
 594         for (ch = 0; ch < src->ch_count; ch++) {
 595             memcpy(dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
 596                    src->ch[ch] + ((n - *out_sz) * c->felem_size), c->felem_size);
 597         }
 598     }
 599
 600     // if not enough data is in, return and wait for more
 601     if (num < c->filter_length + 1) {
 602         *out_sz = num;
 603         *out_idx = c->filter_length;
 604         return INT_MAX;
 605     }
 606
 607     // else invert
 608     for (n = 1; n <= c->filter_length; n++) {
 609         for (ch = 0; ch < src->ch_count; ch++) {
 610             memcpy(dst->ch[ch] + ((c->filter_length - n) * c->felem_size),
 611                    dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
 612                    c->felem_size);
 613         }
 614     }
 615
 616     res = num - *out_sz;
 617     *out_idx = c->filter_length;
 618     while (c->index < 0) {
 619         --*out_idx;
 620         c->index += c->phase_count;
 621     }
 622     *out_sz = FFMAX(*out_sz + c->filter_length,
 623                     1 + c->filter_length * 2) - *out_idx;
 624
 625     return FFMAX(res, 0);
 626 }
 627
 628 struct Resampler const swri_resampler={
 629   resample_init,
 630   resample_free,
 631   multiple_resample,
 632   resample_flush,
 633   set_compensation,
 634   get_delay,
 635   invert_initial_buffer,
 636   get_out_samples,
 637 };