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
 153     if (!tab || !sin_lut)
 154         goto fail;
 155
 156     /* if upsampling, only need to interpolate, no filter */
 157     if (factor > 1.0)
 158         factor = 1.0;
 159
 160     if (factor == 1.0) {
 161         for (ph = 0; ph < ph_nb; ph++)
 162             sin_lut[ph] = sin(M_PI * ph / phase_count);
 163     }
 164     for(ph = 0; ph < ph_nb; ph++) {
 165         double norm = 0;
 166         s = sin_lut[ph];
 167         for(i=0;i<=tap_count;i++) {
 168             x = M_PI * ((double)(i - center) - (double)ph / phase_count) * factor;
 169             if (x == 0) y = 1.0;
 170             else if (factor == 1.0)
 171                 y = s / x;
 172             else
 173                 y = sin(x) / x;
 174             switch(filter_type){
 175             case SWR_FILTER_TYPE_CUBIC:{
 176                 const float d= -0.5; //first order derivative = -0.5
 177                 x = fabs(((double)(i - center) - (double)ph / phase_count) * factor);
 178                 if(x<1.0) y= 1 - 3*x*x + 2*x*x*x + d*(            -x*x + x*x*x);
 179                 else      y=                       d*(-4 + 8*x - 5*x*x + x*x*x);
 180                 break;}
 181             case SWR_FILTER_TYPE_BLACKMAN_NUTTALL:
 182                 w = 2.0*x / (factor*tap_count);
 183                 t = -cos(w);
 184                 y *= 0.3635819 - 0.4891775 * t + 0.1365995 * (2*t*t-1) - 0.0106411 * (4*t*t*t - 3*t);
 185                 break;
 186             case SWR_FILTER_TYPE_KAISER:
 187                 w = 2.0*x / (factor*tap_count*M_PI);
 188                 y *= bessel(kaiser_beta*sqrt(FFMAX(1-w*w, 0)));
 189                 break;
 190             default:
 191                 av_assert0(0);
 192             }
 193
 194             tab[i] = y;
 195             s = -s;
 196             if (i < tap_count)
 197                 norm += y;
 198         }
 199
 200         /* normalize so that an uniform color remains the same */
 201         switch(c->format){
 202         case AV_SAMPLE_FMT_S16P:
 203             for(i=0;i<tap_count;i++)
 204                 ((int16_t*)filter)[ph * alloc + i] = av_clip_int16(lrintf(tab[i] * scale / norm));
 205             if (phase_count % 2) break;
 206             if (tap_count % 2 == 0 || tap_count == 1) {
 207                 for (i = 0; i < tap_count; i++)
 208                     ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int16_t*)filter)[ph * alloc + i];
 209             }
 210             else {
 211                 for (i = 1; i <= tap_count; i++)
 212                     ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
 213                         av_clip_int16(lrintf(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
 214             }
 215             break;
 216         case AV_SAMPLE_FMT_S32P:
 217             for(i=0;i<tap_count;i++)
 218                 ((int32_t*)filter)[ph * alloc + i] = av_clipl_int32(llrint(tab[i] * scale / norm));
 219             if (phase_count % 2) break;
 220             if (tap_count % 2 == 0 || tap_count == 1) {
 221                 for (i = 0; i < tap_count; i++)
 222                     ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int32_t*)filter)[ph * alloc + i];
 223             }
 224             else {
 225                 for (i = 1; i <= tap_count; i++)
 226                     ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
 227                         av_clipl_int32(llrint(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
 228             }
 229             break;
 230         case AV_SAMPLE_FMT_FLTP:
 231             for(i=0;i<tap_count;i++)
 232                 ((float*)filter)[ph * alloc + i] = tab[i] * scale / norm;
 233             if (phase_count % 2) break;
 234             if (tap_count % 2 == 0 || tap_count == 1) {
 235                 for (i = 0; i < tap_count; i++)
 236                     ((float*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((float*)filter)[ph * alloc + i];
 237             }
 238             else {
 239                 for (i = 1; i <= tap_count; i++)
 240                     ((float*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
 241             }
 242             break;
 243         case AV_SAMPLE_FMT_DBLP:
 244             for(i=0;i<tap_count;i++)
 245                 ((double*)filter)[ph * alloc + i] = tab[i] * scale / norm;
 246             if (phase_count % 2) break;
 247             if (tap_count % 2 == 0 || tap_count == 1) {
 248                 for (i = 0; i < tap_count; i++)
 249                     ((double*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((double*)filter)[ph * alloc + i];
 250             }
 251             else {
 252                 for (i = 1; i <= tap_count; i++)
 253                     ((double*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
 254             }
 255             break;
 256         }
 257     }
 258 #if 0
 259     {
 260 #define LEN 1024
 261         int j,k;
 262         double sine[LEN + tap_count];
 263         double filtered[LEN];
 264         double maxff=-2, minff=2, maxsf=-2, minsf=2;
 265         for(i=0; i<LEN; i++){
 266             double ss=0, sf=0, ff=0;
 267             for(j=0; j<LEN+tap_count; j++)
 268                 sine[j]= cos(i*j*M_PI/LEN);
 269             for(j=0; j<LEN; j++){
 270                 double sum=0;
 271                 ph=0;
 272                 for(k=0; k<tap_count; k++)
 273                     sum += filter[ph * tap_count + k] * sine[k+j];
 274                 filtered[j]= sum / (1<<FILTER_SHIFT);
 275                 ss+= sine[j + center] * sine[j + center];
 276                 ff+= filtered[j] * filtered[j];
 277                 sf+= sine[j + center] * filtered[j];
 278             }
 279             ss= sqrt(2*ss/LEN);
 280             ff= sqrt(2*ff/LEN);
 281             sf= 2*sf/LEN;
 282             maxff= FFMAX(maxff, ff);
 283             minff= FFMIN(minff, ff);
 284             maxsf= FFMAX(maxsf, sf);
 285             minsf= FFMIN(minsf, sf);
 286             if(i%11==0){
 287                 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);
 288                 minff=minsf= 2;
 289                 maxff=maxsf= -2;
 290             }
 291         }
 292     }
 293 #endif
 294
 295 fail:
 296     av_free(tab);
 297     av_free(sin_lut);
 298     return 0;
 299 }
 300
 301 static ResampleContext *resample_init(ResampleContext *c, int out_rate, int in_rate, int filter_size, int phase_shift, int linear,
 302                                     double cutoff0, enum AVSampleFormat format, enum SwrFilterType filter_type, double kaiser_beta,
 303                                     double precision, int cheby, int exact_rational)
 304 {
 305     double cutoff = cutoff0? cutoff0 : 0.97;
 306     double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
 307     int phase_count= 1<<phase_shift;
 308     int phase_count_compensation = phase_count;
 309
 310     if (exact_rational) {
 311         int phase_count_exact, phase_count_exact_den;
 312
 313         av_reduce(&phase_count_exact, &phase_count_exact_den, out_rate, in_rate, INT_MAX);
 314         if (phase_count_exact <= phase_count) {
 315             phase_count_compensation = phase_count_exact * (phase_count / phase_count_exact);
 316             phase_count = phase_count_exact;
 317         }
 318     }
 319
 320     if (!c || c->phase_count != phase_count || c->linear!=linear || c->factor != factor
 321            || c->filter_length != FFMAX((int)ceil(filter_size/factor), 1) || c->format != format
 322            || c->filter_type != filter_type || c->kaiser_beta != kaiser_beta) {
 323         c = av_mallocz(sizeof(*c));
 324         if (!c)
 325             return NULL;
 326
 327         c->format= format;
 328
 329         c->felem_size= av_get_bytes_per_sample(c->format);
 330
 331         switch(c->format){
 332         case AV_SAMPLE_FMT_S16P:
 333             c->filter_shift = 15;
 334             break;
 335         case AV_SAMPLE_FMT_S32P:
 336             c->filter_shift = 30;
 337             break;
 338         case AV_SAMPLE_FMT_FLTP:
 339         case AV_SAMPLE_FMT_DBLP:
 340             c->filter_shift = 0;
 341             break;
 342         default:
 343             av_log(NULL, AV_LOG_ERROR, "Unsupported sample format\n");
 344             av_assert0(0);
 345         }
 346
 347         if (filter_size/factor > INT32_MAX/256) {
 348             av_log(NULL, AV_LOG_ERROR, "Filter length too large\n");
 349             goto error;
 350         }
 351
 352         c->phase_count   = phase_count;
 353         c->linear        = linear;
 354         c->factor        = factor;
 355         c->filter_length = FFMAX((int)ceil(filter_size/factor), 1);
 356         c->filter_alloc  = FFALIGN(c->filter_length, 8);
 357         c->filter_bank   = av_calloc(c->filter_alloc, (phase_count+1)*c->felem_size);
 358         c->filter_type   = filter_type;
 359         c->kaiser_beta   = kaiser_beta;
 360         c->phase_count_compensation = phase_count_compensation;
 361         if (!c->filter_bank)
 362             goto error;
 363         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))
 364             goto error;
 365         memcpy(c->filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, c->filter_bank, (c->filter_alloc-1)*c->felem_size);
 366         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);
 367     }
 368
 369     c->compensation_distance= 0;
 370     if(!av_reduce(&c->src_incr, &c->dst_incr, out_rate, in_rate * (int64_t)phase_count, INT32_MAX/2))
 371         goto error;
 372     while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
 373         c->dst_incr *= 2;
 374         c->src_incr *= 2;
 375     }
 376     c->ideal_dst_incr = c->dst_incr;
 377     c->dst_incr_div   = c->dst_incr / c->src_incr;
 378     c->dst_incr_mod   = c->dst_incr % c->src_incr;
 379
 380     c->index= -phase_count*((c->filter_length-1)/2);
 381     c->frac= 0;
 382
 383     swri_resample_dsp_init(c);
 384
 385     return c;
 386 error:
 387     av_freep(&c->filter_bank);
 388     av_free(c);
 389     return NULL;
 390 }
 391
 392 static void resample_free(ResampleContext **c){
 393     if(!*c)
 394         return;
 395     av_freep(&(*c)->filter_bank);
 396     av_freep(c);
 397 }
 398
 399 static int rebuild_filter_bank_with_compensation(ResampleContext *c)
 400 {
 401     uint8_t *new_filter_bank;
 402     int new_src_incr, new_dst_incr;
 403     int phase_count = c->phase_count_compensation;
 404     int ret;
 405
 406     if (phase_count == c->phase_count)
 407         return 0;
 408
 409     av_assert0(!c->frac && !c->dst_incr_mod && !c->compensation_distance);
 410
 411     new_filter_bank = av_calloc(c->filter_alloc, (phase_count + 1) * c->felem_size);
 412     if (!new_filter_bank)
 413         return AVERROR(ENOMEM);
 414
 415     ret = build_filter(c, new_filter_bank, c->factor, c->filter_length, c->filter_alloc,
 416                        phase_count, 1 << c->filter_shift, c->filter_type, c->kaiser_beta);
 417     if (ret < 0) {
 418         av_freep(&new_filter_bank);
 419         return ret;
 420     }
 421     memcpy(new_filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, new_filter_bank, (c->filter_alloc-1)*c->felem_size);
 422     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);
 423
 424     if (!av_reduce(&new_src_incr, &new_dst_incr, c->src_incr,
 425                    c->dst_incr * (int64_t)(phase_count/c->phase_count), INT32_MAX/2))
 426     {
 427         av_freep(&new_filter_bank);
 428         return AVERROR(EINVAL);
 429     }
 430
 431     c->src_incr = new_src_incr;
 432     c->dst_incr = new_dst_incr;
 433     while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
 434         c->dst_incr *= 2;
 435         c->src_incr *= 2;
 436     }
 437     c->ideal_dst_incr = c->dst_incr;
 438     c->dst_incr_div   = c->dst_incr / c->src_incr;
 439     c->dst_incr_mod   = c->dst_incr % c->src_incr;
 440     c->index         *= phase_count / c->phase_count;
 441     c->phase_count    = phase_count;
 442     av_freep(&c->filter_bank);
 443     c->filter_bank = new_filter_bank;
 444     return 0;
 445 }
 446
 447 static int set_compensation(ResampleContext *c, int sample_delta, int compensation_distance){
 448     int ret;
 449
 450     if (compensation_distance) {
 451         ret = rebuild_filter_bank_with_compensation(c);
 452         if (ret < 0)
 453             return ret;
 454     }
 455
 456     c->compensation_distance= compensation_distance;
 457     if (compensation_distance)
 458         c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
 459     else
 460         c->dst_incr = c->ideal_dst_incr;
 461
 462     c->dst_incr_div   = c->dst_incr / c->src_incr;
 463     c->dst_incr_mod   = c->dst_incr % c->src_incr;
 464
 465     return 0;
 466 }
 467
 468 static int swri_resample(ResampleContext *c,
 469                          uint8_t *dst, const uint8_t *src, int *consumed,
 470                          int src_size, int dst_size, int update_ctx)
 471 {
 472     if (c->filter_length == 1 && c->phase_count == 1) {
 473         int index= c->index;
 474         int frac= c->frac;
 475         int64_t index2= (1LL<<32)*c->frac/c->src_incr + (1LL<<32)*index;
 476         int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
 477         int new_size = (src_size * (int64_t)c->src_incr - frac + c->dst_incr - 1) / c->dst_incr;
 478
 479         dst_size= FFMIN(dst_size, new_size);
 480         c->dsp.resample_one(dst, src, dst_size, index2, incr);
 481
 482         index += dst_size * c->dst_incr_div;
 483         index += (frac + dst_size * (int64_t)c->dst_incr_mod) / c->src_incr;
 484         av_assert2(index >= 0);
 485         *consumed= index;
 486         if (update_ctx) {
 487             c->frac   = (frac + dst_size * (int64_t)c->dst_incr_mod) % c->src_incr;
 488             c->index = 0;
 489         }
 490     } else {
 491         int64_t end_index = (1LL + src_size - c->filter_length) * c->phase_count;
 492         int64_t delta_frac = (end_index - c->index) * c->src_incr - c->frac;
 493         int delta_n = (delta_frac + c->dst_incr - 1) / c->dst_incr;
 494
 495         dst_size = FFMIN(dst_size, delta_n);
 496         if (dst_size > 0) {
 497             *consumed = c->dsp.resample(c, dst, src, dst_size, update_ctx);
 498         } else {
 499             *consumed = 0;
 500         }
 501     }
 502
 503     return dst_size;
 504 }
 505
 506 static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, AudioData *src, int src_size, int *consumed){
 507     int i, ret= -1;
 508     int av_unused mm_flags = av_get_cpu_flags();
 509     int need_emms = c->format == AV_SAMPLE_FMT_S16P && ARCH_X86_32 &&
 510                     (mm_flags & (AV_CPU_FLAG_MMX2 | AV_CPU_FLAG_SSE2)) == AV_CPU_FLAG_MMX2;
 511     int64_t max_src_size = (INT64_MAX/2 / c->phase_count) / c->src_incr;
 512
 513     if (c->compensation_distance)
 514         dst_size = FFMIN(dst_size, c->compensation_distance);
 515     src_size = FFMIN(src_size, max_src_size);
 516
 517     for(i=0; i<dst->ch_count; i++){
 518         ret= swri_resample(c, dst->ch[i], src->ch[i],
 519                            consumed, src_size, dst_size, i+1==dst->ch_count);
 520     }
 521     if(need_emms)
 522         emms_c();
 523
 524     if (c->compensation_distance) {
 525         c->compensation_distance -= ret;
 526         if (!c->compensation_distance) {
 527             c->dst_incr     = c->ideal_dst_incr;
 528             c->dst_incr_div = c->dst_incr / c->src_incr;
 529             c->dst_incr_mod = c->dst_incr % c->src_incr;
 530         }
 531     }
 532
 533     return ret;
 534 }
 535
 536 static int64_t get_delay(struct SwrContext *s, int64_t base){
 537     ResampleContext *c = s->resample;
 538     int64_t num = s->in_buffer_count - (c->filter_length-1)/2;
 539     num *= c->phase_count;
 540     num -= c->index;
 541     num *= c->src_incr;
 542     num -= c->frac;
 543     return av_rescale(num, base, s->in_sample_rate*(int64_t)c->src_incr * c->phase_count);
 544 }
 545
 546 static int64_t get_out_samples(struct SwrContext *s, int in_samples) {
 547     ResampleContext *c = s->resample;
 548     // The + 2 are added to allow implementations to be slightly inaccurate, they should not be needed currently.
 549     // They also make it easier to proof that changes and optimizations do not
 550     // break the upper bound.
 551     int64_t num = s->in_buffer_count + 2LL + in_samples;
 552     num *= c->phase_count;
 553     num -= c->index;
 554     num = av_rescale_rnd(num, s->out_sample_rate, ((int64_t)s->in_sample_rate) * c->phase_count, AV_ROUND_UP) + 2;
 555
 556     if (c->compensation_distance) {
 557         if (num > INT_MAX)
 558             return AVERROR(EINVAL);
 559
 560         num = FFMAX(num, (num * c->ideal_dst_incr - 1) / c->dst_incr + 1);
 561     }
 562     return num;
 563 }
 564
 565 static int resample_flush(struct SwrContext *s) {
 566     AudioData *a= &s->in_buffer;
 567     int i, j, ret;
 568     if((ret = swri_realloc_audio(a, s->in_buffer_index + 2*s->in_buffer_count)) < 0)
 569         return ret;
 570     av_assert0(a->planar);
 571     for(i=0; i<a->ch_count; i++){
 572         for(j=0; j<s->in_buffer_count; j++){
 573             memcpy(a->ch[i] + (s->in_buffer_index+s->in_buffer_count+j  )*a->bps,
 574                 a->ch[i] + (s->in_buffer_index+s->in_buffer_count-j-1)*a->bps, a->bps);
 575         }
 576     }
 577     s->in_buffer_count += (s->in_buffer_count+1)/2;
 578     return 0;
 579 }
 580
 581 // in fact the whole handle multiple ridiculously small buffers might need more thinking...
 582 static int invert_initial_buffer(ResampleContext *c, AudioData *dst, const AudioData *src,
 583                                  int in_count, int *out_idx, int *out_sz)
 584 {
 585     int n, ch, num = FFMIN(in_count + *out_sz, c->filter_length + 1), res;
 586
 587     if (c->index >= 0)
 588         return 0;
 589
 590     if ((res = swri_realloc_audio(dst, c->filter_length * 2 + 1)) < 0)
 591         return res;
 592
 593     // copy
 594     for (n = *out_sz; n < num; n++) {
 595         for (ch = 0; ch < src->ch_count; ch++) {
 596             memcpy(dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
 597                    src->ch[ch] + ((n - *out_sz) * c->felem_size), c->felem_size);
 598         }
 599     }
 600
 601     // if not enough data is in, return and wait for more
 602     if (num < c->filter_length + 1) {
 603         *out_sz = num;
 604         *out_idx = c->filter_length;
 605         return INT_MAX;
 606     }
 607
 608     // else invert
 609     for (n = 1; n <= c->filter_length; n++) {
 610         for (ch = 0; ch < src->ch_count; ch++) {
 611             memcpy(dst->ch[ch] + ((c->filter_length - n) * c->felem_size),
 612                    dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
 613                    c->felem_size);
 614         }
 615     }
 616
 617     res = num - *out_sz;
 618     *out_idx = c->filter_length;
 619     while (c->index < 0) {
 620         --*out_idx;
 621         c->index += c->phase_count;
 622     }
 623     *out_sz = FFMAX(*out_sz + c->filter_length,
 624                     1 + c->filter_length * 2) - *out_idx;
 625
 626     return FFMAX(res, 0);
 627 }
 628
 629 struct Resampler const swri_resampler={
 630   resample_init,
 631   resample_free,
 632   multiple_resample,
 633   resample_flush,
 634   set_compensation,
 635   get_delay,
 636   invert_initial_buffer,
 637   get_out_samples,
 638 };