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