resample.c 23.2 KB
Newer Older
M
Michael Niedermayer 已提交
1 2
/*
 * audio resampling
3
 * Copyright (c) 2004-2012 Michael Niedermayer <michaelni@gmx.at>
4
 * bessel function: Copyright (c) 2006 Xiaogang Zhang
M
Michael Niedermayer 已提交
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
 *
 * This file is part of FFmpeg.
 *
 * FFmpeg is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * FFmpeg is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with FFmpeg; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

/**
 * @file
 * audio resampling
 * @author Michael Niedermayer <michaelni@gmx.at>
 */

29
#include "libavutil/avassert.h"
R
Ronald S. Bultje 已提交
30
#include "resample.h"
M
Michael Niedermayer 已提交
31

32 33 34 35 36 37 38 39 40 41
static inline double eval_poly(const double *coeff, int size, double x) {
    double sum = coeff[size-1];
    int i;
    for (i = size-2; i >= 0; --i) {
        sum *= x;
        sum += coeff[i];
    }
    return sum;
}

M
Michael Niedermayer 已提交
42 43
/**
 * 0th order modified bessel function of the first kind.
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
 * Algorithm taken from the Boost project, source:
 * https://searchcode.com/codesearch/view/14918379/
 * Use, modification and distribution are subject to the
 * Boost Software License, Version 1.0 (see notice below).
 * Boost Software License - Version 1.0 - August 17th, 2003
Permission is hereby granted, free of charge, to any person or organization
obtaining a copy of the software and accompanying documentation covered by
this license (the "Software") to use, reproduce, display, distribute,
execute, and transmit the Software, and to prepare derivative works of the
Software, and to permit third-parties to whom the Software is furnished to
do so, all subject to the following:

The copyright notices in the Software and this entire statement, including
the above license grant, this restriction and the following disclaimer,
must be included in all copies of the Software, in whole or in part, and
all derivative works of the Software, unless such copies or derivative
works are solely in the form of machine-executable object code generated by
a source language processor.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
M
Michael Niedermayer 已提交
70 71
 */

72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
static double bessel(double x) {
// Modified Bessel function of the first kind of order zero
// minimax rational approximations on intervals, see
// Blair and Edwards, Chalk River Report AECL-4928, 1974
    static const double p1[] = {
        -2.2335582639474375249e+15,
        -5.5050369673018427753e+14,
        -3.2940087627407749166e+13,
        -8.4925101247114157499e+11,
        -1.1912746104985237192e+10,
        -1.0313066708737980747e+08,
        -5.9545626019847898221e+05,
        -2.4125195876041896775e+03,
        -7.0935347449210549190e+00,
        -1.5453977791786851041e-02,
        -2.5172644670688975051e-05,
        -3.0517226450451067446e-08,
        -2.6843448573468483278e-11,
        -1.5982226675653184646e-14,
        -5.2487866627945699800e-18,
    };
    static const double q1[] = {
        -2.2335582639474375245e+15,
         7.8858692566751002988e+12,
        -1.2207067397808979846e+10,
         1.0377081058062166144e+07,
        -4.8527560179962773045e+03,
99
         1.0,
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
    };
    static const double p2[] = {
        -2.2210262233306573296e-04,
         1.3067392038106924055e-02,
        -4.4700805721174453923e-01,
         5.5674518371240761397e+00,
        -2.3517945679239481621e+01,
         3.1611322818701131207e+01,
        -9.6090021968656180000e+00,
    };
    static const double q2[] = {
        -5.5194330231005480228e-04,
         3.2547697594819615062e-02,
        -1.1151759188741312645e+00,
         1.3982595353892851542e+01,
        -6.0228002066743340583e+01,
         8.5539563258012929600e+01,
        -3.1446690275135491500e+01,
118
        1.0,
119 120 121 122 123 124 125 126 127 128 129 130 131 132
    };
    double y, r, factor;
    if (x == 0)
        return 1.0;
    x = fabs(x);
    if (x <= 15) {
        y = x * x;
        return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
    }
    else {
        y = 1 / x - 1.0 / 15;
        r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
        factor = exp(x) / sqrt(x);
        return factor * r;
M
Michael Niedermayer 已提交
133 134 135 136 137 138 139
    }
}

/**
 * builds a polyphase filterbank.
 * @param factor resampling factor
 * @param scale wanted sum of coefficients for each filter
140 141
 * @param filter_type  filter type
 * @param kaiser_beta  kaiser window beta
M
Michael Niedermayer 已提交
142 143
 * @return 0 on success, negative on error
 */
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){
M
Michael Niedermayer 已提交
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));
M
Michael Niedermayer 已提交
151
    const int center= (tap_count-1)/2;
152
    int ret = AVERROR(ENOMEM);
M
Michael Niedermayer 已提交
153

154 155
    if (!tab || !sin_lut)
        goto fail;
M
Michael Niedermayer 已提交
156 157 158 159 160

    /* if upsampling, only need to interpolate, no filter */
    if (factor > 1.0)
        factor = 1.0;

161
    if (factor == 1.0) {
162
        for (ph = 0; ph < ph_nb; ph++)
163 164
            sin_lut[ph] = sin(M_PI * ph / phase_count);
    }
165
    for(ph = 0; ph < ph_nb; ph++) {
M
Michael Niedermayer 已提交
166
        double norm = 0;
167
        s = sin_lut[ph];
168
        for(i=0;i<=tap_count;i++) {
M
Michael Niedermayer 已提交
169 170
            x = M_PI * ((double)(i - center) - (double)ph / phase_count) * factor;
            if (x == 0) y = 1.0;
171 172 173 174
            else if (factor == 1.0)
                y = s / x;
            else
                y = sin(x) / x;
175 176
            switch(filter_type){
            case SWR_FILTER_TYPE_CUBIC:{
M
Michael Niedermayer 已提交
177 178 179 180 181
                const float d= -0.5; //first order derivative = -0.5
                x = fabs(((double)(i - center) - (double)ph / phase_count) * factor);
                if(x<1.0) y= 1 - 3*x*x + 2*x*x*x + d*(            -x*x + x*x*x);
                else      y=                       d*(-4 + 8*x - 5*x*x + x*x*x);
                break;}
182
            case SWR_FILTER_TYPE_BLACKMAN_NUTTALL:
183 184
                w = 2.0*x / (factor*tap_count);
                t = -cos(w);
185
                y *= 0.3635819 - 0.4891775 * t + 0.1365995 * (2*t*t-1) - 0.0106411 * (4*t*t*t - 3*t);
M
Michael Niedermayer 已提交
186
                break;
187
            case SWR_FILTER_TYPE_KAISER:
M
Michael Niedermayer 已提交
188
                w = 2.0*x / (factor*tap_count*M_PI);
189
                y *= bessel(kaiser_beta*sqrt(FFMAX(1-w*w, 0)));
M
Michael Niedermayer 已提交
190
                break;
191 192
            default:
                av_assert0(0);
M
Michael Niedermayer 已提交
193 194 195
            }

            tab[i] = y;
196
            s = -s;
197 198
            if (i < tap_count)
                norm += y;
M
Michael Niedermayer 已提交
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 209 210 211 212 213
                for (i = 0; i < tap_count; i++)
                    ((int16_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int16_t*)filter)[ph * alloc + i];
            }
            else {
                for (i = 1; i <= tap_count; i++)
                    ((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 223 224 225 226 227 228 229
                for (i = 0; i < tap_count; i++)
                    ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((int32_t*)filter)[ph * alloc + i];
            }
            else {
                for (i = 1; i <= tap_count; i++)
                    ((int32_t*)filter)[(phase_count-ph) * alloc + tap_count-i] =
                        av_clipl_int32(llrint(tab[i] * scale / (norm - tab[0] + tab[tap_count])));
            }
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 237 238 239 240 241 242
                for (i = 0; i < tap_count; i++)
                    ((float*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((float*)filter)[ph * alloc + i];
            }
            else {
                for (i = 1; i <= tap_count; i++)
                    ((float*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
            }
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 250 251 252 253 254 255
                for (i = 0; i < tap_count; i++)
                    ((double*)filter)[(phase_count-ph) * alloc + tap_count-1-i] = ((double*)filter)[ph * alloc + i];
            }
            else {
                for (i = 1; i <= tap_count; i++)
                    ((double*)filter)[(phase_count-ph) * alloc + tap_count-i] = tab[i] * scale / (norm - tab[0] + tab[tap_count]);
            }
256
            break;
M
Michael Niedermayer 已提交
257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295
        }
    }
#if 0
    {
#define LEN 1024
        int j,k;
        double sine[LEN + tap_count];
        double filtered[LEN];
        double maxff=-2, minff=2, maxsf=-2, minsf=2;
        for(i=0; i<LEN; i++){
            double ss=0, sf=0, ff=0;
            for(j=0; j<LEN+tap_count; j++)
                sine[j]= cos(i*j*M_PI/LEN);
            for(j=0; j<LEN; j++){
                double sum=0;
                ph=0;
                for(k=0; k<tap_count; k++)
                    sum += filter[ph * tap_count + k] * sine[k+j];
                filtered[j]= sum / (1<<FILTER_SHIFT);
                ss+= sine[j + center] * sine[j + center];
                ff+= filtered[j] * filtered[j];
                sf+= sine[j + center] * filtered[j];
            }
            ss= sqrt(2*ss/LEN);
            ff= sqrt(2*ff/LEN);
            sf= 2*sf/LEN;
            maxff= FFMAX(maxff, ff);
            minff= FFMIN(minff, ff);
            maxsf= FFMAX(maxsf, sf);
            minsf= FFMIN(minsf, sf);
            if(i%11==0){
                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);
                minff=minsf= 2;
                maxff=maxsf= -2;
            }
        }
    }
#endif

296
    ret = 0;
297
fail:
M
Michael Niedermayer 已提交
298
    av_free(tab);
299
    av_free(sin_lut);
300
    return ret;
M
Michael Niedermayer 已提交
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)
R
Ronald S. Bultje 已提交
306
{
307
    double cutoff = cutoff0? cutoff0 : 0.97;
M
Michael Niedermayer 已提交
308 309
    double factor= FFMIN(out_rate * cutoff / in_rate, 1.0);
    int phase_count= 1<<phase_shift;
310
    int phase_count_compensation = phase_count;
M
Michael Niedermayer 已提交
311

312 313 314 315 316
    if (exact_rational) {
        int phase_count_exact, phase_count_exact_den;

        av_reduce(&phase_count_exact, &phase_count_exact_den, out_rate, in_rate, INT_MAX);
        if (phase_count_exact <= phase_count) {
317
            phase_count_compensation = phase_count_exact * (phase_count / phase_count_exact);
318 319 320 321 322
            phase_count = phase_count_exact;
        }
    }

    if (!c || c->phase_count != phase_count || c->linear!=linear || c->factor != factor
323 324
           || c->filter_length != FFMAX((int)ceil(filter_size/factor), 1) || c->format != format
           || c->filter_type != filter_type || c->kaiser_beta != kaiser_beta) {
325
        c = av_mallocz(sizeof(*c));
326 327 328
        if (!c)
            return NULL;

329 330
        c->format= format;

331 332
        c->felem_size= av_get_bytes_per_sample(c->format);

333
        switch(c->format){
334
        case AV_SAMPLE_FMT_S16P:
335 336
            c->filter_shift = 15;
            break;
337
        case AV_SAMPLE_FMT_S32P:
338 339
            c->filter_shift = 30;
            break;
340 341
        case AV_SAMPLE_FMT_FLTP:
        case AV_SAMPLE_FMT_DBLP:
342 343
            c->filter_shift = 0;
            break;
344 345
        default:
            av_log(NULL, AV_LOG_ERROR, "Unsupported sample format\n");
346
            av_assert0(0);
347 348
        }

349 350 351 352 353
        if (filter_size/factor > INT32_MAX/256) {
            av_log(NULL, AV_LOG_ERROR, "Filter length too large\n");
            goto error;
        }

354
        c->phase_count   = phase_count;
355 356 357
        c->linear        = linear;
        c->factor        = factor;
        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 361
        c->filter_type   = filter_type;
        c->kaiser_beta   = kaiser_beta;
362
        c->phase_count_compensation = phase_count_compensation;
363 364
        if (!c->filter_bank)
            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 368
        memcpy(c->filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, c->filter_bank, (c->filter_alloc-1)*c->felem_size);
        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);
M
Michael Niedermayer 已提交
369 370 371
    }

    c->compensation_distance= 0;
M
Michael Niedermayer 已提交
372 373
    if(!av_reduce(&c->src_incr, &c->dst_incr, out_rate, in_rate * (int64_t)phase_count, INT32_MAX/2))
        goto error;
374 375 376 377
    while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
        c->dst_incr *= 2;
        c->src_incr *= 2;
    }
378 379 380
    c->ideal_dst_incr = c->dst_incr;
    c->dst_incr_div   = c->dst_incr / c->src_incr;
    c->dst_incr_mod   = c->dst_incr % c->src_incr;
M
Michael Niedermayer 已提交
381

M
Michael Niedermayer 已提交
382 383 384
    c->index= -phase_count*((c->filter_length-1)/2);
    c->frac= 0;

385
    swri_resample_dsp_init(c);
R
Ronald S. Bultje 已提交
386

M
Michael Niedermayer 已提交
387 388
    return c;
error:
389
    av_freep(&c->filter_bank);
M
Michael Niedermayer 已提交
390 391 392 393
    av_free(c);
    return NULL;
}

394
static void resample_free(ResampleContext **c){
M
Michael Niedermayer 已提交
395 396 397 398 399 400
    if(!*c)
        return;
    av_freep(&(*c)->filter_bank);
    av_freep(c);
}

401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448
static int rebuild_filter_bank_with_compensation(ResampleContext *c)
{
    uint8_t *new_filter_bank;
    int new_src_incr, new_dst_incr;
    int phase_count = c->phase_count_compensation;
    int ret;

    if (phase_count == c->phase_count)
        return 0;

    av_assert0(!c->frac && !c->dst_incr_mod && !c->compensation_distance);

    new_filter_bank = av_calloc(c->filter_alloc, (phase_count + 1) * c->felem_size);
    if (!new_filter_bank)
        return AVERROR(ENOMEM);

    ret = build_filter(c, new_filter_bank, c->factor, c->filter_length, c->filter_alloc,
                       phase_count, 1 << c->filter_shift, c->filter_type, c->kaiser_beta);
    if (ret < 0) {
        av_freep(&new_filter_bank);
        return ret;
    }
    memcpy(new_filter_bank + (c->filter_alloc*phase_count+1)*c->felem_size, new_filter_bank, (c->filter_alloc-1)*c->felem_size);
    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);

    if (!av_reduce(&new_src_incr, &new_dst_incr, c->src_incr,
                   c->dst_incr * (int64_t)(phase_count/c->phase_count), INT32_MAX/2))
    {
        av_freep(&new_filter_bank);
        return AVERROR(EINVAL);
    }

    c->src_incr = new_src_incr;
    c->dst_incr = new_dst_incr;
    while (c->dst_incr < (1<<20) && c->src_incr < (1<<20)) {
        c->dst_incr *= 2;
        c->src_incr *= 2;
    }
    c->ideal_dst_incr = c->dst_incr;
    c->dst_incr_div   = c->dst_incr / c->src_incr;
    c->dst_incr_mod   = c->dst_incr % c->src_incr;
    c->index         *= phase_count / c->phase_count;
    c->phase_count    = phase_count;
    av_freep(&c->filter_bank);
    c->filter_bank = new_filter_bank;
    return 0;
}

449
static int set_compensation(ResampleContext *c, int sample_delta, int compensation_distance){
450 451
    int ret;

452
    if (compensation_distance && sample_delta) {
453 454 455 456 457
        ret = rebuild_filter_bank_with_compensation(c);
        if (ret < 0)
            return ret;
    }

M
Michael Niedermayer 已提交
458
    c->compensation_distance= compensation_distance;
459 460 461 462
    if (compensation_distance)
        c->dst_incr = c->ideal_dst_incr - c->ideal_dst_incr * (int64_t)sample_delta / compensation_distance;
    else
        c->dst_incr = c->ideal_dst_incr;
463 464 465 466

    c->dst_incr_div   = c->dst_incr / c->src_incr;
    c->dst_incr_mod   = c->dst_incr % c->src_incr;

467
    return 0;
M
Michael Niedermayer 已提交
468 469
}

R
Ronald S. Bultje 已提交
470 471 472 473
static int swri_resample(ResampleContext *c,
                         uint8_t *dst, const uint8_t *src, int *consumed,
                         int src_size, int dst_size, int update_ctx)
{
474
    if (c->filter_length == 1 && c->phase_count == 1) {
R
Ronald S. Bultje 已提交
475 476 477 478 479 480 481
        int index= c->index;
        int frac= c->frac;
        int64_t index2= (1LL<<32)*c->frac/c->src_incr + (1LL<<32)*index;
        int64_t incr= (1LL<<32) * c->dst_incr / c->src_incr;
        int new_size = (src_size * (int64_t)c->src_incr - frac + c->dst_incr - 1) / c->dst_incr;

        dst_size= FFMIN(dst_size, new_size);
482
        c->dsp.resample_one(dst, src, dst_size, index2, incr);
R
Ronald S. Bultje 已提交
483

484 485
        index += dst_size * c->dst_incr_div;
        index += (frac + dst_size * (int64_t)c->dst_incr_mod) / c->src_incr;
R
Ronald S. Bultje 已提交
486 487 488
        av_assert2(index >= 0);
        *consumed= index;
        if (update_ctx) {
489
            c->frac   = (frac + dst_size * (int64_t)c->dst_incr_mod) % c->src_incr;
R
Ronald S. Bultje 已提交
490 491 492
            c->index = 0;
        }
    } else {
493
        int64_t end_index = (1LL + src_size - c->filter_length) * c->phase_count;
R
Ronald S. Bultje 已提交
494 495 496 497
        int64_t delta_frac = (end_index - c->index) * c->src_incr - c->frac;
        int delta_n = (delta_frac + c->dst_incr - 1) / c->dst_incr;

        dst_size = FFMIN(dst_size, delta_n);
498
        if (dst_size > 0) {
499 500 501 502 503 504
            /* resample_linear and resample_common should have same behavior
             * when frac and dst_incr_mod are zero */
            if (c->linear && (c->frac || c->dst_incr_mod))
                *consumed = c->dsp.resample_linear(c, dst, src, dst_size, update_ctx);
            else
                *consumed = c->dsp.resample_common(c, dst, src, dst_size, update_ctx);
R
Ronald S. Bultje 已提交
505
        } else {
506
            *consumed = 0;
R
Ronald S. Bultje 已提交
507 508
        }
    }
509

R
Ronald S. Bultje 已提交
510 511
    return dst_size;
}
M
Michael Niedermayer 已提交
512

513
static int multiple_resample(ResampleContext *c, AudioData *dst, int dst_size, AudioData *src, int src_size, int *consumed){
M
Michael Niedermayer 已提交
514
    int i, ret= -1;
515
    int av_unused mm_flags = av_get_cpu_flags();
R
Ronald S. Bultje 已提交
516 517
    int need_emms = c->format == AV_SAMPLE_FMT_S16P && ARCH_X86_32 &&
                    (mm_flags & (AV_CPU_FLAG_MMX2 | AV_CPU_FLAG_SSE2)) == AV_CPU_FLAG_MMX2;
518
    int64_t max_src_size = (INT64_MAX/2 / c->phase_count) / c->src_incr;
M
Michael Niedermayer 已提交
519

520 521
    if (c->compensation_distance)
        dst_size = FFMIN(dst_size, c->compensation_distance);
522
    src_size = FFMIN(src_size, max_src_size);
523

M
Michael Niedermayer 已提交
524
    for(i=0; i<dst->ch_count; i++){
R
Ronald S. Bultje 已提交
525 526
        ret= swri_resample(c, dst->ch[i], src->ch[i],
                           consumed, src_size, dst_size, i+1==dst->ch_count);
M
Michael Niedermayer 已提交
527
    }
528 529
    if(need_emms)
        emms_c();
530 531 532

    if (c->compensation_distance) {
        c->compensation_distance -= ret;
533 534 535 536 537
        if (!c->compensation_distance) {
            c->dst_incr     = c->ideal_dst_incr;
            c->dst_incr_div = c->dst_incr / c->src_incr;
            c->dst_incr_mod = c->dst_incr % c->src_incr;
        }
538 539
    }

M
Michael Niedermayer 已提交
540 541
    return ret;
}
542

543
static int64_t get_delay(struct SwrContext *s, int64_t base){
544
    ResampleContext *c = s->resample;
545
    int64_t num = s->in_buffer_count - (c->filter_length-1)/2;
546
    num *= c->phase_count;
547 548 549
    num -= c->index;
    num *= c->src_incr;
    num -= c->frac;
550
    return av_rescale(num, base, s->in_sample_rate*(int64_t)c->src_incr * c->phase_count);
551
}
552

553 554
static int64_t get_out_samples(struct SwrContext *s, int in_samples) {
    ResampleContext *c = s->resample;
555
    // The + 2 are added to allow implementations to be slightly inaccurate, they should not be needed currently.
556
    // They also make it easier to proof that changes and optimizations do not
557
    // break the upper bound.
558
    int64_t num = s->in_buffer_count + 2LL + in_samples;
559
    num *= c->phase_count;
560
    num -= c->index;
561
    num = av_rescale_rnd(num, s->out_sample_rate, ((int64_t)s->in_sample_rate) * c->phase_count, AV_ROUND_UP) + 2;
562 563 564 565 566 567 568 569 570 571

    if (c->compensation_distance) {
        if (num > INT_MAX)
            return AVERROR(EINVAL);

        num = FFMAX(num, (num * c->ideal_dst_incr - 1) / c->dst_incr + 1);
    }
    return num;
}

572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587
static int resample_flush(struct SwrContext *s) {
    AudioData *a= &s->in_buffer;
    int i, j, ret;
    if((ret = swri_realloc_audio(a, s->in_buffer_index + 2*s->in_buffer_count)) < 0)
        return ret;
    av_assert0(a->planar);
    for(i=0; i<a->ch_count; i++){
        for(j=0; j<s->in_buffer_count; j++){
            memcpy(a->ch[i] + (s->in_buffer_index+s->in_buffer_count+j  )*a->bps,
                a->ch[i] + (s->in_buffer_index+s->in_buffer_count-j-1)*a->bps, a->bps);
        }
    }
    s->in_buffer_count += (s->in_buffer_count+1)/2;
    return 0;
}

588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624
// in fact the whole handle multiple ridiculously small buffers might need more thinking...
static int invert_initial_buffer(ResampleContext *c, AudioData *dst, const AudioData *src,
                                 int in_count, int *out_idx, int *out_sz)
{
    int n, ch, num = FFMIN(in_count + *out_sz, c->filter_length + 1), res;

    if (c->index >= 0)
        return 0;

    if ((res = swri_realloc_audio(dst, c->filter_length * 2 + 1)) < 0)
        return res;

    // copy
    for (n = *out_sz; n < num; n++) {
        for (ch = 0; ch < src->ch_count; ch++) {
            memcpy(dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
                   src->ch[ch] + ((n - *out_sz) * c->felem_size), c->felem_size);
        }
    }

    // if not enough data is in, return and wait for more
    if (num < c->filter_length + 1) {
        *out_sz = num;
        *out_idx = c->filter_length;
        return INT_MAX;
    }

    // else invert
    for (n = 1; n <= c->filter_length; n++) {
        for (ch = 0; ch < src->ch_count; ch++) {
            memcpy(dst->ch[ch] + ((c->filter_length - n) * c->felem_size),
                   dst->ch[ch] + ((c->filter_length + n) * c->felem_size),
                   c->felem_size);
        }
    }

    res = num - *out_sz;
625 626 627 628 629
    *out_idx = c->filter_length;
    while (c->index < 0) {
        --*out_idx;
        c->index += c->phase_count;
    }
630 631
    *out_sz = FFMAX(*out_sz + c->filter_length,
                    1 + c->filter_length * 2) - *out_idx;
632

633
    return FFMAX(res, 0);
634 635
}

636 637 638 639
struct Resampler const swri_resampler={
  resample_init,
  resample_free,
  multiple_resample,
640
  resample_flush,
641 642
  set_compensation,
  get_delay,
643
  invert_initial_buffer,
644
  get_out_samples,
645
};