git.sesse.net Git - ffmpeg/blob - libavfilter/af_atempo.c

   1 /*
   2  * Copyright (c) 2012 Pavel Koshevoy <pkoshevoy at gmail dot com>
   3  *
   4  * This file is part of FFmpeg.
   5  *
   6  * FFmpeg is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Lesser General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2.1 of the License, or (at your option) any later version.
  10  *
  11  * FFmpeg is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public
  17  * License along with FFmpeg; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  19  */
  20
  21 /**
  22  * @file
  23  * tempo scaling audio filter -- an implementation of WSOLA algorithm
  24  *
  25  * Based on MIT licensed yaeAudioTempoFilter.h and yaeAudioFragment.h
  26  * from Apprentice Video player by Pavel Koshevoy.
  27  * https://sourceforge.net/projects/apprenticevideo/
  28  *
  29  * An explanation of SOLA algorithm is available at
  30  * http://www.surina.net/article/time-and-pitch-scaling.html
  31  *
  32  * WSOLA is very similar to SOLA, only one major difference exists between
  33  * these algorithms.  SOLA shifts audio fragments along the output stream,
  34  * where as WSOLA shifts audio fragments along the input stream.
  35  *
  36  * The advantage of WSOLA algorithm is that the overlap region size is
  37  * always the same, therefore the blending function is constant and
  38  * can be precomputed.
  39  */
  40
  41 #include <float.h>
  42 #include "libavcodec/avfft.h"
  43 #include "libavutil/avassert.h"
  44 #include "libavutil/avstring.h"
  45 #include "libavutil/channel_layout.h"
  46 #include "libavutil/eval.h"
  47 #include "libavutil/opt.h"
  48 #include "libavutil/samplefmt.h"
  49 #include "avfilter.h"
  50 #include "audio.h"
  51 #include "internal.h"
  52
  53 /**
  54  * A fragment of audio waveform
  55  */
  56 typedef struct {
  57     // index of the first sample of this fragment in the overall waveform;
  58     // 0: input sample position
  59     // 1: output sample position
  60     int64_t position[2];
  61
  62     // original packed multi-channel samples:
  63     uint8_t *data;
  64
  65     // number of samples in this fragment:
  66     int nsamples;
  67
  68     // rDFT transform of the down-mixed mono fragment, used for
  69     // fast waveform alignment via correlation in frequency domain:
  70     FFTSample *xdat;
  71 } AudioFragment;
  72
  73 /**
  74  * Filter state machine states
  75  */
  76 typedef enum {
  77     YAE_LOAD_FRAGMENT,
  78     YAE_ADJUST_POSITION,
  79     YAE_RELOAD_FRAGMENT,
  80     YAE_OUTPUT_OVERLAP_ADD,
  81     YAE_FLUSH_OUTPUT,
  82 } FilterState;
  83
  84 /**
  85  * Filter state machine
  86  */
  87 typedef struct {
  88     const AVClass *class;
  89
  90     // ring-buffer of input samples, necessary because some times
  91     // input fragment position may be adjusted backwards:
  92     uint8_t *buffer;
  93
  94     // ring-buffer maximum capacity, expressed in sample rate time base:
  95     int ring;
  96
  97     // ring-buffer house keeping:
  98     int size;
  99     int head;
 100     int tail;
 101
 102     // 0: input sample position corresponding to the ring buffer tail
 103     // 1: output sample position
 104     int64_t position[2];
 105
 106     // sample format:
 107     enum AVSampleFormat format;
 108
 109     // number of channels:
 110     int channels;
 111
 112     // row of bytes to skip from one sample to next, across multple channels;
 113     // stride = (number-of-channels * bits-per-sample-per-channel) / 8
 114     int stride;
 115
 116     // fragment window size, power-of-two integer:
 117     int window;
 118
 119     // Hann window coefficients, for feathering
 120     // (blending) the overlapping fragment region:
 121     float *hann;
 122
 123     // tempo scaling factor:
 124     double tempo;
 125
 126     // cumulative alignment drift:
 127     int drift;
 128
 129     // current/previous fragment ring-buffer:
 130     AudioFragment frag[2];
 131
 132     // current fragment index:
 133     uint64_t nfrag;
 134
 135     // current state:
 136     FilterState state;
 137
 138     // for fast correlation calculation in frequency domain:
 139     RDFTContext *real_to_complex;
 140     RDFTContext *complex_to_real;
 141     FFTSample *correlation;
 142
 143     // for managing AVFilterPad.request_frame and AVFilterPad.filter_frame
 144     AVFrame *dst_buffer;
 145     uint8_t *dst;
 146     uint8_t *dst_end;
 147     uint64_t nsamples_in;
 148     uint64_t nsamples_out;
 149 } ATempoContext;
 150
 151 #define OFFSET(x) offsetof(ATempoContext, x)
 152
 153 static const AVOption atempo_options[] = {
 154     { "tempo", "set tempo scale factor",
 155       OFFSET(tempo), AV_OPT_TYPE_DOUBLE, { .dbl = 1.0 }, 0.5, 2.0,
 156       AV_OPT_FLAG_AUDIO_PARAM | AV_OPT_FLAG_FILTERING_PARAM },
 157     { NULL }
 158 };
 159
 160 AVFILTER_DEFINE_CLASS(atempo);
 161
 162 /**
 163  * Reset filter to initial state, do not deallocate existing local buffers.
 164  */
 165 static void yae_clear(ATempoContext *atempo)
 166 {
 167     atempo->size = 0;
 168     atempo->head = 0;
 169     atempo->tail = 0;
 170
 171     atempo->drift = 0;
 172     atempo->nfrag = 0;
 173     atempo->state = YAE_LOAD_FRAGMENT;
 174
 175     atempo->position[0] = 0;
 176     atempo->position[1] = 0;
 177
 178     atempo->frag[0].position[0] = 0;
 179     atempo->frag[0].position[1] = 0;
 180     atempo->frag[0].nsamples    = 0;
 181
 182     atempo->frag[1].position[0] = 0;
 183     atempo->frag[1].position[1] = 0;
 184     atempo->frag[1].nsamples    = 0;
 185
 186     // shift left position of 1st fragment by half a window
 187     // so that no re-normalization would be required for
 188     // the left half of the 1st fragment:
 189     atempo->frag[0].position[0] = -(int64_t)(atempo->window / 2);
 190     atempo->frag[0].position[1] = -(int64_t)(atempo->window / 2);
 191
 192     av_frame_free(&atempo->dst_buffer);
 193     atempo->dst     = NULL;
 194     atempo->dst_end = NULL;
 195
 196     atempo->nsamples_in       = 0;
 197     atempo->nsamples_out      = 0;
 198 }
 199
 200 /**
 201  * Reset filter to initial state and deallocate all buffers.
 202  */
 203 static void yae_release_buffers(ATempoContext *atempo)
 204 {
 205     yae_clear(atempo);
 206
 207     av_freep(&atempo->frag[0].data);
 208     av_freep(&atempo->frag[1].data);
 209     av_freep(&atempo->frag[0].xdat);
 210     av_freep(&atempo->frag[1].xdat);
 211
 212     av_freep(&atempo->buffer);
 213     av_freep(&atempo->hann);
 214     av_freep(&atempo->correlation);
 215
 216     av_rdft_end(atempo->real_to_complex);
 217     atempo->real_to_complex = NULL;
 218
 219     av_rdft_end(atempo->complex_to_real);
 220     atempo->complex_to_real = NULL;
 221 }
 222
 223 /* av_realloc is not aligned enough; fortunately, the data does not need to
 224  * be preserved */
 225 #define RE_MALLOC_OR_FAIL(field, field_size)                    \
 226     do {                                                        \
 227         av_freep(&field);                                       \
 228         field = av_malloc(field_size);                          \
 229         if (!field) {                                           \
 230             yae_release_buffers(atempo);                        \
 231             return AVERROR(ENOMEM);                             \
 232         }                                                       \
 233     } while (0)
 234
 235 /**
 236  * Prepare filter for processing audio data of given format,
 237  * sample rate and number of channels.
 238  */
 239 static int yae_reset(ATempoContext *atempo,
 240                      enum AVSampleFormat format,
 241                      int sample_rate,
 242                      int channels)
 243 {
 244     const int sample_size = av_get_bytes_per_sample(format);
 245     uint32_t nlevels  = 0;
 246     uint32_t pot;
 247     int i;
 248
 249     atempo->format   = format;
 250     atempo->channels = channels;
 251     atempo->stride   = sample_size * channels;
 252
 253     // pick a segment window size:
 254     atempo->window = sample_rate / 24;
 255
 256     // adjust window size to be a power-of-two integer:
 257     nlevels = av_log2(atempo->window);
 258     pot = 1 << nlevels;
 259     av_assert0(pot <= atempo->window);
 260
 261     if (pot < atempo->window) {
 262         atempo->window = pot * 2;
 263         nlevels++;
 264     }
 265
 266     // initialize audio fragment buffers:
 267     RE_MALLOC_OR_FAIL(atempo->frag[0].data, atempo->window * atempo->stride);
 268     RE_MALLOC_OR_FAIL(atempo->frag[1].data, atempo->window * atempo->stride);
 269     RE_MALLOC_OR_FAIL(atempo->frag[0].xdat, atempo->window * sizeof(FFTComplex));
 270     RE_MALLOC_OR_FAIL(atempo->frag[1].xdat, atempo->window * sizeof(FFTComplex));
 271
 272     // initialize rDFT contexts:
 273     av_rdft_end(atempo->real_to_complex);
 274     atempo->real_to_complex = NULL;
 275
 276     av_rdft_end(atempo->complex_to_real);
 277     atempo->complex_to_real = NULL;
 278
 279     atempo->real_to_complex = av_rdft_init(nlevels + 1, DFT_R2C);
 280     if (!atempo->real_to_complex) {
 281         yae_release_buffers(atempo);
 282         return AVERROR(ENOMEM);
 283     }
 284
 285     atempo->complex_to_real = av_rdft_init(nlevels + 1, IDFT_C2R);
 286     if (!atempo->complex_to_real) {
 287         yae_release_buffers(atempo);
 288         return AVERROR(ENOMEM);
 289     }
 290
 291     RE_MALLOC_OR_FAIL(atempo->correlation, atempo->window * sizeof(FFTComplex));
 292
 293     atempo->ring = atempo->window * 3;
 294     RE_MALLOC_OR_FAIL(atempo->buffer, atempo->ring * atempo->stride);
 295
 296     // initialize the Hann window function:
 297     RE_MALLOC_OR_FAIL(atempo->hann, atempo->window * sizeof(float));
 298
 299     for (i = 0; i < atempo->window; i++) {
 300         double t = (double)i / (double)(atempo->window - 1);
 301         double h = 0.5 * (1.0 - cos(2.0 * M_PI * t));
 302         atempo->hann[i] = (float)h;
 303     }
 304
 305     yae_clear(atempo);
 306     return 0;
 307 }
 308
 309 static int yae_set_tempo(AVFilterContext *ctx, const char *arg_tempo)
 310 {
 311     ATempoContext *atempo = ctx->priv;
 312     char   *tail = NULL;
 313     double tempo = av_strtod(arg_tempo, &tail);
 314
 315     if (tail && *tail) {
 316         av_log(ctx, AV_LOG_ERROR, "Invalid tempo value '%s'\n", arg_tempo);
 317         return AVERROR(EINVAL);
 318     }
 319
 320     if (tempo < 0.5 || tempo > 2.0) {
 321         av_log(ctx, AV_LOG_ERROR, "Tempo value %f exceeds [0.5, 2.0] range\n",
 322                tempo);
 323         return AVERROR(EINVAL);
 324     }
 325
 326     atempo->tempo = tempo;
 327     return 0;
 328 }
 329
 330 inline static AudioFragment *yae_curr_frag(ATempoContext *atempo)
 331 {
 332     return &atempo->frag[atempo->nfrag % 2];
 333 }
 334
 335 inline static AudioFragment *yae_prev_frag(ATempoContext *atempo)
 336 {
 337     return &atempo->frag[(atempo->nfrag + 1) % 2];
 338 }
 339
 340 /**
 341  * A helper macro for initializing complex data buffer with scalar data
 342  * of a given type.
 343  */
 344 #define yae_init_xdat(scalar_type, scalar_max)                          \
 345     do {                                                                \
 346         const uint8_t *src_end = src +                                  \
 347             frag->nsamples * atempo->channels * sizeof(scalar_type);    \
 348                                                                         \
 349         FFTSample *xdat = frag->xdat;                                   \
 350         scalar_type tmp;                                                \
 351                                                                         \
 352         if (atempo->channels == 1) {                                    \
 353             for (; src < src_end; xdat++) {                             \
 354                 tmp = *(const scalar_type *)src;                        \
 355                 src += sizeof(scalar_type);                             \
 356                                                                         \
 357                 *xdat = (FFTSample)tmp;                                 \
 358             }                                                           \
 359         } else {                                                        \
 360             FFTSample s, max, ti, si;                                   \
 361             int i;                                                      \
 362                                                                         \
 363             for (; src < src_end; xdat++) {                             \
 364                 tmp = *(const scalar_type *)src;                        \
 365                 src += sizeof(scalar_type);                             \
 366                                                                         \
 367                 max = (FFTSample)tmp;                                   \
 368                 s = FFMIN((FFTSample)scalar_max,                        \
 369                           (FFTSample)fabsf(max));                       \
 370                                                                         \
 371                 for (i = 1; i < atempo->channels; i++) {                \
 372                     tmp = *(const scalar_type *)src;                    \
 373                     src += sizeof(scalar_type);                         \
 374                                                                         \
 375                     ti = (FFTSample)tmp;                                \
 376                     si = FFMIN((FFTSample)scalar_max,                   \
 377                                (FFTSample)fabsf(ti));                   \
 378                                                                         \
 379                     if (s < si) {                                       \
 380                         s   = si;                                       \
 381                         max = ti;                                       \
 382                     }                                                   \
 383                 }                                                       \
 384                                                                         \
 385                 *xdat = max;                                            \
 386             }                                                           \
 387         }                                                               \
 388     } while (0)
 389
 390 /**
 391  * Initialize complex data buffer of a given audio fragment
 392  * with down-mixed mono data of appropriate scalar type.
 393  */
 394 static void yae_downmix(ATempoContext *atempo, AudioFragment *frag)
 395 {
 396     // shortcuts:
 397     const uint8_t *src = frag->data;
 398
 399     // init complex data buffer used for FFT and Correlation:
 400     memset(frag->xdat, 0, sizeof(FFTComplex) * atempo->window);
 401
 402     if (atempo->format == AV_SAMPLE_FMT_U8) {
 403         yae_init_xdat(uint8_t, 127);
 404     } else if (atempo->format == AV_SAMPLE_FMT_S16) {
 405         yae_init_xdat(int16_t, 32767);
 406     } else if (atempo->format == AV_SAMPLE_FMT_S32) {
 407         yae_init_xdat(int, 2147483647);
 408     } else if (atempo->format == AV_SAMPLE_FMT_FLT) {
 409         yae_init_xdat(float, 1);
 410     } else if (atempo->format == AV_SAMPLE_FMT_DBL) {
 411         yae_init_xdat(double, 1);
 412     }
 413 }
 414
 415 /**
 416  * Populate the internal data buffer on as-needed basis.
 417  *
 418  * @return
 419  *   0 if requested data was already available or was successfully loaded,
 420  *   AVERROR(EAGAIN) if more input data is required.
 421  */
 422 static int yae_load_data(ATempoContext *atempo,
 423                          const uint8_t **src_ref,
 424                          const uint8_t *src_end,
 425                          int64_t stop_here)
 426 {
 427     // shortcut:
 428     const uint8_t *src = *src_ref;
 429     const int read_size = stop_here - atempo->position[0];
 430
 431     if (stop_here <= atempo->position[0]) {
 432         return 0;
 433     }
 434
 435     // samples are not expected to be skipped:
 436     av_assert0(read_size <= atempo->ring);
 437
 438     while (atempo->position[0] < stop_here && src < src_end) {
 439         int src_samples = (src_end - src) / atempo->stride;
 440
 441         // load data piece-wise, in order to avoid complicating the logic:
 442         int nsamples = FFMIN(read_size, src_samples);
 443         int na;
 444         int nb;
 445
 446         nsamples = FFMIN(nsamples, atempo->ring);
 447         na = FFMIN(nsamples, atempo->ring - atempo->tail);
 448         nb = FFMIN(nsamples - na, atempo->ring);
 449
 450         if (na) {
 451             uint8_t *a = atempo->buffer + atempo->tail * atempo->stride;
 452             memcpy(a, src, na * atempo->stride);
 453
 454             src += na * atempo->stride;
 455             atempo->position[0] += na;
 456
 457             atempo->size = FFMIN(atempo->size + na, atempo->ring);
 458             atempo->tail = (atempo->tail + na) % atempo->ring;
 459             atempo->head =
 460                 atempo->size < atempo->ring ?
 461                 atempo->tail - atempo->size :
 462                 atempo->tail;
 463         }
 464
 465         if (nb) {
 466             uint8_t *b = atempo->buffer;
 467             memcpy(b, src, nb * atempo->stride);
 468
 469             src += nb * atempo->stride;
 470             atempo->position[0] += nb;
 471
 472             atempo->size = FFMIN(atempo->size + nb, atempo->ring);
 473             atempo->tail = (atempo->tail + nb) % atempo->ring;
 474             atempo->head =
 475                 atempo->size < atempo->ring ?
 476                 atempo->tail - atempo->size :
 477                 atempo->tail;
 478         }
 479     }
 480
 481     // pass back the updated source buffer pointer:
 482     *src_ref = src;
 483
 484     // sanity check:
 485     av_assert0(atempo->position[0] <= stop_here);
 486
 487     return atempo->position[0] == stop_here ? 0 : AVERROR(EAGAIN);
 488 }
 489
 490 /**
 491  * Populate current audio fragment data buffer.
 492  *
 493  * @return
 494  *   0 when the fragment is ready,
 495  *   AVERROR(EAGAIN) if more input data is required.
 496  */
 497 static int yae_load_frag(ATempoContext *atempo,
 498                          const uint8_t **src_ref,
 499                          const uint8_t *src_end)
 500 {
 501     // shortcuts:
 502     AudioFragment *frag = yae_curr_frag(atempo);
 503     uint8_t *dst;
 504     int64_t missing, start, zeros;
 505     uint32_t nsamples;
 506     const uint8_t *a, *b;
 507     int i0, i1, n0, n1, na, nb;
 508
 509     int64_t stop_here = frag->position[0] + atempo->window;
 510     if (src_ref && yae_load_data(atempo, src_ref, src_end, stop_here) != 0) {
 511         return AVERROR(EAGAIN);
 512     }
 513
 514     // calculate the number of samples we don't have:
 515     missing =
 516         stop_here > atempo->position[0] ?
 517         stop_here - atempo->position[0] : 0;
 518
 519     nsamples =
 520         missing < (int64_t)atempo->window ?
 521         (uint32_t)(atempo->window - missing) : 0;
 522
 523     // setup the output buffer:
 524     frag->nsamples = nsamples;
 525     dst = frag->data;
 526
 527     start = atempo->position[0] - atempo->size;
 528     zeros = 0;
 529
 530     if (frag->position[0] < start) {
 531         // what we don't have we substitute with zeros:
 532         zeros = FFMIN(start - frag->position[0], (int64_t)nsamples);
 533         av_assert0(zeros != nsamples);
 534
 535         memset(dst, 0, zeros * atempo->stride);
 536         dst += zeros * atempo->stride;
 537     }
 538
 539     if (zeros == nsamples) {
 540         return 0;
 541     }
 542
 543     // get the remaining data from the ring buffer:
 544     na = (atempo->head < atempo->tail ?
 545           atempo->tail - atempo->head :
 546           atempo->ring - atempo->head);
 547
 548     nb = atempo->head < atempo->tail ? 0 : atempo->tail;
 549
 550     // sanity check:
 551     av_assert0(nsamples <= zeros + na + nb);
 552
 553     a = atempo->buffer + atempo->head * atempo->stride;
 554     b = atempo->buffer;
 555
 556     i0 = frag->position[0] + zeros - start;
 557     i1 = i0 < na ? 0 : i0 - na;
 558
 559     n0 = i0 < na ? FFMIN(na - i0, (int)(nsamples - zeros)) : 0;
 560     n1 = nsamples - zeros - n0;
 561
 562     if (n0) {
 563         memcpy(dst, a + i0 * atempo->stride, n0 * atempo->stride);
 564         dst += n0 * atempo->stride;
 565     }
 566
 567     if (n1) {
 568         memcpy(dst, b + i1 * atempo->stride, n1 * atempo->stride);
 569     }
 570
 571     return 0;
 572 }
 573
 574 /**
 575  * Prepare for loading next audio fragment.
 576  */
 577 static void yae_advance_to_next_frag(ATempoContext *atempo)
 578 {
 579     const double fragment_step = atempo->tempo * (double)(atempo->window / 2);
 580
 581     const AudioFragment *prev;
 582     AudioFragment       *frag;
 583
 584     atempo->nfrag++;
 585     prev = yae_prev_frag(atempo);
 586     frag = yae_curr_frag(atempo);
 587
 588     frag->position[0] = prev->position[0] + (int64_t)fragment_step;
 589     frag->position[1] = prev->position[1] + atempo->window / 2;
 590     frag->nsamples    = 0;
 591 }
 592
 593 /**
 594  * Calculate cross-correlation via rDFT.
 595  *
 596  * Multiply two vectors of complex numbers (result of real_to_complex rDFT)
 597  * and transform back via complex_to_real rDFT.
 598  */
 599 static void yae_xcorr_via_rdft(FFTSample *xcorr,
 600                                RDFTContext *complex_to_real,
 601                                const FFTComplex *xa,
 602                                const FFTComplex *xb,
 603                                const int window)
 604 {
 605     FFTComplex *xc = (FFTComplex *)xcorr;
 606     int i;
 607
 608     // NOTE: first element requires special care -- Given Y = rDFT(X),
 609     // Im(Y[0]) and Im(Y[N/2]) are always zero, therefore av_rdft_calc
 610     // stores Re(Y[N/2]) in place of Im(Y[0]).
 611
 612     xc->re = xa->re * xb->re;
 613     xc->im = xa->im * xb->im;
 614     xa++;
 615     xb++;
 616     xc++;
 617
 618     for (i = 1; i < window; i++, xa++, xb++, xc++) {
 619         xc->re = (xa->re * xb->re + xa->im * xb->im);
 620         xc->im = (xa->im * xb->re - xa->re * xb->im);
 621     }
 622
 623     // apply inverse rDFT:
 624     av_rdft_calc(complex_to_real, xcorr);
 625 }
 626
 627 /**
 628  * Calculate alignment offset for given fragment
 629  * relative to the previous fragment.
 630  *
 631  * @return alignment offset of current fragment relative to previous.
 632  */
 633 static int yae_align(AudioFragment *frag,
 634                      const AudioFragment *prev,
 635                      const int window,
 636                      const int delta_max,
 637                      const int drift,
 638                      FFTSample *correlation,
 639                      RDFTContext *complex_to_real)
 640 {
 641     int       best_offset = -drift;
 642     FFTSample best_metric = -FLT_MAX;
 643     FFTSample *xcorr;
 644
 645     int i0;
 646     int i1;
 647     int i;
 648
 649     yae_xcorr_via_rdft(correlation,
 650                        complex_to_real,
 651                        (const FFTComplex *)prev->xdat,
 652                        (const FFTComplex *)frag->xdat,
 653                        window);
 654
 655     // identify search window boundaries:
 656     i0 = FFMAX(window / 2 - delta_max - drift, 0);
 657     i0 = FFMIN(i0, window);
 658
 659     i1 = FFMIN(window / 2 + delta_max - drift, window - window / 16);
 660     i1 = FFMAX(i1, 0);
 661
 662     // identify cross-correlation peaks within search window:
 663     xcorr = correlation + i0;
 664
 665     for (i = i0; i < i1; i++, xcorr++) {
 666         FFTSample metric = *xcorr;
 667
 668         // normalize:
 669         FFTSample drifti = (FFTSample)(drift + i);
 670         metric *= drifti * (FFTSample)(i - i0) * (FFTSample)(i1 - i);
 671
 672         if (metric > best_metric) {
 673             best_metric = metric;
 674             best_offset = i - window / 2;
 675         }
 676     }
 677
 678     return best_offset;
 679 }
 680
 681 /**
 682  * Adjust current fragment position for better alignment
 683  * with previous fragment.
 684  *
 685  * @return alignment correction.
 686  */
 687 static int yae_adjust_position(ATempoContext *atempo)
 688 {
 689     const AudioFragment *prev = yae_prev_frag(atempo);
 690     AudioFragment       *frag = yae_curr_frag(atempo);
 691
 692     const int delta_max  = atempo->window / 2;
 693     const int correction = yae_align(frag,
 694                                      prev,
 695                                      atempo->window,
 696                                      delta_max,
 697                                      atempo->drift,
 698                                      atempo->correlation,
 699                                      atempo->complex_to_real);
 700
 701     if (correction) {
 702         // adjust fragment position:
 703         frag->position[0] -= correction;
 704
 705         // clear so that the fragment can be reloaded:
 706         frag->nsamples = 0;
 707
 708         // update cumulative correction drift counter:
 709         atempo->drift += correction;
 710     }
 711
 712     return correction;
 713 }
 714
 715 /**
 716  * A helper macro for blending the overlap region of previous
 717  * and current audio fragment.
 718  */
 719 #define yae_blend(scalar_type)                                          \
 720     do {                                                                \
 721         const scalar_type *aaa = (const scalar_type *)a;                \
 722         const scalar_type *bbb = (const scalar_type *)b;                \
 723                                                                         \
 724         scalar_type *out     = (scalar_type *)dst;                      \
 725         scalar_type *out_end = (scalar_type *)dst_end;                  \
 726         int64_t i;                                                      \
 727                                                                         \
 728         for (i = 0; i < overlap && out < out_end;                       \
 729              i++, atempo->position[1]++, wa++, wb++) {                  \
 730             float w0 = *wa;                                             \
 731             float w1 = *wb;                                             \
 732             int j;                                                      \
 733                                                                         \
 734             for (j = 0; j < atempo->channels;                           \
 735                  j++, aaa++, bbb++, out++) {                            \
 736                 float t0 = (float)*aaa;                                 \
 737                 float t1 = (float)*bbb;                                 \
 738                                                                         \
 739                 *out =                                                  \
 740                     frag->position[0] + i < 0 ?                         \
 741                     *aaa :                                              \
 742                     (scalar_type)(t0 * w0 + t1 * w1);                   \
 743             }                                                           \
 744         }                                                               \
 745         dst = (uint8_t *)out;                                           \
 746     } while (0)
 747
 748 /**
 749  * Blend the overlap region of previous and current audio fragment
 750  * and output the results to the given destination buffer.
 751  *
 752  * @return
 753  *   0 if the overlap region was completely stored in the dst buffer,
 754  *   AVERROR(EAGAIN) if more destination buffer space is required.
 755  */
 756 static int yae_overlap_add(ATempoContext *atempo,
 757                            uint8_t **dst_ref,
 758                            uint8_t *dst_end)
 759 {
 760     // shortcuts:
 761     const AudioFragment *prev = yae_prev_frag(atempo);
 762     const AudioFragment *frag = yae_curr_frag(atempo);
 763
 764     const int64_t start_here = FFMAX(atempo->position[1],
 765                                      frag->position[1]);
 766
 767     const int64_t stop_here = FFMIN(prev->position[1] + prev->nsamples,
 768                                     frag->position[1] + frag->nsamples);
 769
 770     const int64_t overlap = stop_here - start_here;
 771
 772     const int64_t ia = start_here - prev->position[1];
 773     const int64_t ib = start_here - frag->position[1];
 774
 775     const float *wa = atempo->hann + ia;
 776     const float *wb = atempo->hann + ib;
 777
 778     const uint8_t *a = prev->data + ia * atempo->stride;
 779     const uint8_t *b = frag->data + ib * atempo->stride;
 780
 781     uint8_t *dst = *dst_ref;
 782
 783     av_assert0(start_here <= stop_here &&
 784                frag->position[1] <= start_here &&
 785                overlap <= frag->nsamples);
 786
 787     if (atempo->format == AV_SAMPLE_FMT_U8) {
 788         yae_blend(uint8_t);
 789     } else if (atempo->format == AV_SAMPLE_FMT_S16) {
 790         yae_blend(int16_t);
 791     } else if (atempo->format == AV_SAMPLE_FMT_S32) {
 792         yae_blend(int);
 793     } else if (atempo->format == AV_SAMPLE_FMT_FLT) {
 794         yae_blend(float);
 795     } else if (atempo->format == AV_SAMPLE_FMT_DBL) {
 796         yae_blend(double);
 797     }
 798
 799     // pass-back the updated destination buffer pointer:
 800     *dst_ref = dst;
 801
 802     return atempo->position[1] == stop_here ? 0 : AVERROR(EAGAIN);
 803 }
 804
 805 /**
 806  * Feed as much data to the filter as it is able to consume
 807  * and receive as much processed data in the destination buffer
 808  * as it is able to produce or store.
 809  */
 810 static void
 811 yae_apply(ATempoContext *atempo,
 812           const uint8_t **src_ref,
 813           const uint8_t *src_end,
 814           uint8_t **dst_ref,
 815           uint8_t *dst_end)
 816 {
 817     while (1) {
 818         if (atempo->state == YAE_LOAD_FRAGMENT) {
 819             // load additional data for the current fragment:
 820             if (yae_load_frag(atempo, src_ref, src_end) != 0) {
 821                 break;
 822             }
 823
 824             // down-mix to mono:
 825             yae_downmix(atempo, yae_curr_frag(atempo));
 826
 827             // apply rDFT:
 828             av_rdft_calc(atempo->real_to_complex, yae_curr_frag(atempo)->xdat);
 829
 830             // must load the second fragment before alignment can start:
 831             if (!atempo->nfrag) {
 832                 yae_advance_to_next_frag(atempo);
 833                 continue;
 834             }
 835
 836             atempo->state = YAE_ADJUST_POSITION;
 837         }
 838
 839         if (atempo->state == YAE_ADJUST_POSITION) {
 840             // adjust position for better alignment:
 841             if (yae_adjust_position(atempo)) {
 842                 // reload the fragment at the corrected position, so that the
 843                 // Hann window blending would not require normalization:
 844                 atempo->state = YAE_RELOAD_FRAGMENT;
 845             } else {
 846                 atempo->state = YAE_OUTPUT_OVERLAP_ADD;
 847             }
 848         }
 849
 850         if (atempo->state == YAE_RELOAD_FRAGMENT) {
 851             // load additional data if necessary due to position adjustment:
 852             if (yae_load_frag(atempo, src_ref, src_end) != 0) {
 853                 break;
 854             }
 855
 856             // down-mix to mono:
 857             yae_downmix(atempo, yae_curr_frag(atempo));
 858
 859             // apply rDFT:
 860             av_rdft_calc(atempo->real_to_complex, yae_curr_frag(atempo)->xdat);
 861
 862             atempo->state = YAE_OUTPUT_OVERLAP_ADD;
 863         }
 864
 865         if (atempo->state == YAE_OUTPUT_OVERLAP_ADD) {
 866             // overlap-add and output the result:
 867             if (yae_overlap_add(atempo, dst_ref, dst_end) != 0) {
 868                 break;
 869             }
 870
 871             // advance to the next fragment, repeat:
 872             yae_advance_to_next_frag(atempo);
 873             atempo->state = YAE_LOAD_FRAGMENT;
 874         }
 875     }
 876 }
 877
 878 /**
 879  * Flush any buffered data from the filter.
 880  *
 881  * @return
 882  *   0 if all data was completely stored in the dst buffer,
 883  *   AVERROR(EAGAIN) if more destination buffer space is required.
 884  */
 885 static int yae_flush(ATempoContext *atempo,
 886                      uint8_t **dst_ref,
 887                      uint8_t *dst_end)
 888 {
 889     AudioFragment *frag = yae_curr_frag(atempo);
 890     int64_t overlap_end;
 891     int64_t start_here;
 892     int64_t stop_here;
 893     int64_t offset;
 894
 895     const uint8_t *src;
 896     uint8_t *dst;
 897
 898     int src_size;
 899     int dst_size;
 900     int nbytes;
 901
 902     atempo->state = YAE_FLUSH_OUTPUT;
 903
 904     if (atempo->position[0] == frag->position[0] + frag->nsamples &&
 905         atempo->position[1] == frag->position[1] + frag->nsamples) {
 906         // the current fragment is already flushed:
 907         return 0;
 908     }
 909
 910     if (frag->position[0] + frag->nsamples < atempo->position[0]) {
 911         // finish loading the current (possibly partial) fragment:
 912         yae_load_frag(atempo, NULL, NULL);
 913
 914         if (atempo->nfrag) {
 915             // down-mix to mono:
 916             yae_downmix(atempo, frag);
 917
 918             // apply rDFT:
 919             av_rdft_calc(atempo->real_to_complex, frag->xdat);
 920
 921             // align current fragment to previous fragment:
 922             if (yae_adjust_position(atempo)) {
 923                 // reload the current fragment due to adjusted position:
 924                 yae_load_frag(atempo, NULL, NULL);
 925             }
 926         }
 927     }
 928
 929     // flush the overlap region:
 930     overlap_end = frag->position[1] + FFMIN(atempo->window / 2,
 931                                             frag->nsamples);
 932
 933     while (atempo->position[1] < overlap_end) {
 934         if (yae_overlap_add(atempo, dst_ref, dst_end) != 0) {
 935             return AVERROR(EAGAIN);
 936         }
 937     }
 938
 939     // flush the remaininder of the current fragment:
 940     start_here = FFMAX(atempo->position[1], overlap_end);
 941     stop_here  = frag->position[1] + frag->nsamples;
 942     offset     = start_here - frag->position[1];
 943     av_assert0(start_here <= stop_here && frag->position[1] <= start_here);
 944
 945     src = frag->data + offset * atempo->stride;
 946     dst = (uint8_t *)*dst_ref;
 947
 948     src_size = (int)(stop_here - start_here) * atempo->stride;
 949     dst_size = dst_end - dst;
 950     nbytes = FFMIN(src_size, dst_size);
 951
 952     memcpy(dst, src, nbytes);
 953     dst += nbytes;
 954
 955     atempo->position[1] += (nbytes / atempo->stride);
 956
 957     // pass-back the updated destination buffer pointer:
 958     *dst_ref = (uint8_t *)dst;
 959
 960     return atempo->position[1] == stop_here ? 0 : AVERROR(EAGAIN);
 961 }
 962
 963 static av_cold int init(AVFilterContext *ctx)
 964 {
 965     ATempoContext *atempo = ctx->priv;
 966     atempo->format = AV_SAMPLE_FMT_NONE;
 967     atempo->state  = YAE_LOAD_FRAGMENT;
 968     return 0;
 969 }
 970
 971 static av_cold void uninit(AVFilterContext *ctx)
 972 {
 973     ATempoContext *atempo = ctx->priv;
 974     yae_release_buffers(atempo);
 975 }
 976
 977 static int query_formats(AVFilterContext *ctx)
 978 {
 979     AVFilterChannelLayouts *layouts = NULL;
 980     AVFilterFormats        *formats = NULL;
 981
 982     // WSOLA necessitates an internal sliding window ring buffer
 983     // for incoming audio stream.
 984     //
 985     // Planar sample formats are too cumbersome to store in a ring buffer,
 986     // therefore planar sample formats are not supported.
 987     //
 988     static const enum AVSampleFormat sample_fmts[] = {
 989         AV_SAMPLE_FMT_U8,
 990         AV_SAMPLE_FMT_S16,
 991         AV_SAMPLE_FMT_S32,
 992         AV_SAMPLE_FMT_FLT,
 993         AV_SAMPLE_FMT_DBL,
 994         AV_SAMPLE_FMT_NONE
 995     };
 996
 997     layouts = ff_all_channel_layouts();
 998     if (!layouts) {
 999         return AVERROR(ENOMEM);
1000     }
1001     ff_set_common_channel_layouts(ctx, layouts);
1002
1003     formats = ff_make_format_list(sample_fmts);
1004     if (!formats) {
1005         return AVERROR(ENOMEM);
1006     }
1007     ff_set_common_formats(ctx, formats);
1008
1009     formats = ff_all_samplerates();
1010     if (!formats) {
1011         return AVERROR(ENOMEM);
1012     }
1013     ff_set_common_samplerates(ctx, formats);
1014
1015     return 0;
1016 }
1017
1018 static int config_props(AVFilterLink *inlink)
1019 {
1020     AVFilterContext  *ctx = inlink->dst;
1021     ATempoContext *atempo = ctx->priv;
1022
1023     enum AVSampleFormat format = inlink->format;
1024     int sample_rate = (int)inlink->sample_rate;
1025     int channels = av_get_channel_layout_nb_channels(inlink->channel_layout);
1026
1027     ctx->outputs[0]->flags |= FF_LINK_FLAG_REQUEST_LOOP;
1028
1029     return yae_reset(atempo, format, sample_rate, channels);
1030 }
1031
1032 static int push_samples(ATempoContext *atempo,
1033                         AVFilterLink *outlink,
1034                         int n_out)
1035 {
1036     int ret;
1037
1038     atempo->dst_buffer->sample_rate = outlink->sample_rate;
1039     atempo->dst_buffer->nb_samples  = n_out;
1040
1041     // adjust the PTS:
1042     atempo->dst_buffer->pts =
1043         av_rescale_q(atempo->nsamples_out,
1044                      (AVRational){ 1, outlink->sample_rate },
1045                      outlink->time_base);
1046
1047     ret = ff_filter_frame(outlink, atempo->dst_buffer);
1048     if (ret < 0)
1049         return ret;
1050     atempo->dst_buffer = NULL;
1051     atempo->dst        = NULL;
1052     atempo->dst_end    = NULL;
1053
1054     atempo->nsamples_out += n_out;
1055     return 0;
1056 }
1057
1058 static int filter_frame(AVFilterLink *inlink, AVFrame *src_buffer)
1059 {
1060     AVFilterContext  *ctx = inlink->dst;
1061     ATempoContext *atempo = ctx->priv;
1062     AVFilterLink *outlink = ctx->outputs[0];
1063
1064     int ret = 0;
1065     int n_in = src_buffer->nb_samples;
1066     int n_out = (int)(0.5 + ((double)n_in) / atempo->tempo);
1067
1068     const uint8_t *src = src_buffer->data[0];
1069     const uint8_t *src_end = src + n_in * atempo->stride;
1070
1071     while (src < src_end) {
1072         if (!atempo->dst_buffer) {
1073             atempo->dst_buffer = ff_get_audio_buffer(outlink, n_out);
1074             if (!atempo->dst_buffer)
1075                 return AVERROR(ENOMEM);
1076             av_frame_copy_props(atempo->dst_buffer, src_buffer);
1077
1078             atempo->dst = atempo->dst_buffer->data[0];
1079             atempo->dst_end = atempo->dst + n_out * atempo->stride;
1080         }
1081
1082         yae_apply(atempo, &src, src_end, &atempo->dst, atempo->dst_end);
1083
1084         if (atempo->dst == atempo->dst_end) {
1085             ret = push_samples(atempo, outlink, n_out);
1086             if (ret < 0)
1087                 goto end;
1088         }
1089     }
1090
1091     atempo->nsamples_in += n_in;
1092 end:
1093     av_frame_free(&src_buffer);
1094     return ret;
1095 }
1096
1097 static int request_frame(AVFilterLink *outlink)
1098 {
1099     AVFilterContext  *ctx = outlink->src;
1100     ATempoContext *atempo = ctx->priv;
1101     int ret;
1102
1103     ret = ff_request_frame(ctx->inputs[0]);
1104
1105     if (ret == AVERROR_EOF) {
1106         // flush the filter:
1107         int n_max = atempo->ring;
1108         int n_out;
1109         int err = AVERROR(EAGAIN);
1110
1111         while (err == AVERROR(EAGAIN)) {
1112             if (!atempo->dst_buffer) {
1113                 atempo->dst_buffer = ff_get_audio_buffer(outlink, n_max);
1114                 if (!atempo->dst_buffer)
1115                     return AVERROR(ENOMEM);
1116
1117                 atempo->dst = atempo->dst_buffer->data[0];
1118                 atempo->dst_end = atempo->dst + n_max * atempo->stride;
1119             }
1120
1121             err = yae_flush(atempo, &atempo->dst, atempo->dst_end);
1122
1123             n_out = ((atempo->dst - atempo->dst_buffer->data[0]) /
1124                      atempo->stride);
1125
1126             if (n_out) {
1127                 ret = push_samples(atempo, outlink, n_out);
1128             }
1129         }
1130
1131         av_frame_free(&atempo->dst_buffer);
1132         atempo->dst     = NULL;
1133         atempo->dst_end = NULL;
1134
1135         return AVERROR_EOF;
1136     }
1137
1138     return ret;
1139 }
1140
1141 static int process_command(AVFilterContext *ctx,
1142                            const char *cmd,
1143                            const char *arg,
1144                            char *res,
1145                            int res_len,
1146                            int flags)
1147 {
1148     return !strcmp(cmd, "tempo") ? yae_set_tempo(ctx, arg) : AVERROR(ENOSYS);
1149 }
1150
1151 static const AVFilterPad atempo_inputs[] = {
1152     {
1153         .name         = "default",
1154         .type         = AVMEDIA_TYPE_AUDIO,
1155         .filter_frame = filter_frame,
1156         .config_props = config_props,
1157     },
1158     { NULL }
1159 };
1160
1161 static const AVFilterPad atempo_outputs[] = {
1162     {
1163         .name          = "default",
1164         .request_frame = request_frame,
1165         .type          = AVMEDIA_TYPE_AUDIO,
1166     },
1167     { NULL }
1168 };
1169
1170 AVFilter avfilter_af_atempo = {
1171     .name            = "atempo",
1172     .description     = NULL_IF_CONFIG_SMALL("Adjust audio tempo."),
1173     .init            = init,
1174     .uninit          = uninit,
1175     .query_formats   = query_formats,
1176     .process_command = process_command,
1177     .priv_size       = sizeof(ATempoContext),
1178     .priv_class      = &atempo_class,
1179     .inputs          = atempo_inputs,
1180     .outputs         = atempo_outputs,
1181 };