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

   1 /*
   2  * Copyright (c) 2011 Jan Kokemüller
   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  * This file is based on libebur128 which is available at
  21  * https://github.com/jiixyj/libebur128/
  22  *
  23  * Libebur128 has the following copyright:
  24  *
  25  * Permission is hereby granted, free of charge, to any person obtaining a copy
  26  * of this software and associated documentation files (the "Software"), to deal
  27  * in the Software without restriction, including without limitation the rights
  28  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  29  * copies of the Software, and to permit persons to whom the Software is
  30  * furnished to do so, subject to the following conditions:
  31  *
  32  * The above copyright notice and this permission notice shall be included in
  33  * all copies or substantial portions of the Software.
  34  *
  35  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  36  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  37  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  38  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  39  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  40  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  41  * THE SOFTWARE.
  42 */
  43
  44 #include "ebur128.h"
  45
  46 #include <float.h>
  47 #include <limits.h>
  48 #include <math.h>               /* You may have to define _USE_MATH_DEFINES if you use MSVC */
  49
  50 #include "libavutil/common.h"
  51 #include "libavutil/mem.h"
  52 #include "libavutil/mem_internal.h"
  53 #include "libavutil/thread.h"
  54
  55 #define CHECK_ERROR(condition, errorcode, goto_point)                          \
  56     if ((condition)) {                                                         \
  57         errcode = (errorcode);                                                 \
  58         goto goto_point;                                                       \
  59     }
  60
  61 #define ALMOST_ZERO 0.000001
  62
  63 #define RELATIVE_GATE         (-10.0)
  64 #define RELATIVE_GATE_FACTOR  pow(10.0, RELATIVE_GATE / 10.0)
  65 #define MINUS_20DB            pow(10.0, -20.0 / 10.0)
  66
  67 struct FFEBUR128StateInternal {
  68     /** Filtered audio data (used as ring buffer). */
  69     double *audio_data;
  70     /** Size of audio_data array. */
  71     size_t audio_data_frames;
  72     /** Current index for audio_data. */
  73     size_t audio_data_index;
  74     /** How many frames are needed for a gating block. Will correspond to 400ms
  75      *  of audio at initialization, and 100ms after the first block (75% overlap
  76      *  as specified in the 2011 revision of BS1770). */
  77     unsigned long needed_frames;
  78     /** The channel map. Has as many elements as there are channels. */
  79     int *channel_map;
  80     /** How many samples fit in 100ms (rounded). */
  81     unsigned long samples_in_100ms;
  82     /** BS.1770 filter coefficients (nominator). */
  83     double b[5];
  84     /** BS.1770 filter coefficients (denominator). */
  85     double a[5];
  86     /** BS.1770 filter state. */
  87     double v[5][5];
  88     /** Histograms, used to calculate LRA. */
  89     unsigned long *block_energy_histogram;
  90     unsigned long *short_term_block_energy_histogram;
  91     /** Keeps track of when a new short term block is needed. */
  92     size_t short_term_frame_counter;
  93     /** Maximum sample peak, one per channel */
  94     double *sample_peak;
  95     /** The maximum window duration in ms. */
  96     unsigned long window;
  97     /** Data pointer array for interleaved data */
  98     void **data_ptrs;
  99 };
 100
 101 static AVOnce histogram_init = AV_ONCE_INIT;
 102 static DECLARE_ALIGNED(32, double, histogram_energies)[1000];
 103 static DECLARE_ALIGNED(32, double, histogram_energy_boundaries)[1001];
 104
 105 static void ebur128_init_filter(FFEBUR128State * st)
 106 {
 107     int i, j;
 108
 109     double f0 = 1681.974450955533;
 110     double G = 3.999843853973347;
 111     double Q = 0.7071752369554196;
 112
 113     double K = tan(M_PI * f0 / (double) st->samplerate);
 114     double Vh = pow(10.0, G / 20.0);
 115     double Vb = pow(Vh, 0.4996667741545416);
 116
 117     double pb[3] = { 0.0, 0.0, 0.0 };
 118     double pa[3] = { 1.0, 0.0, 0.0 };
 119     double rb[3] = { 1.0, -2.0, 1.0 };
 120     double ra[3] = { 1.0, 0.0, 0.0 };
 121
 122     double a0 = 1.0 + K / Q + K * K;
 123     pb[0] = (Vh + Vb * K / Q + K * K) / a0;
 124     pb[1] = 2.0 * (K * K - Vh) / a0;
 125     pb[2] = (Vh - Vb * K / Q + K * K) / a0;
 126     pa[1] = 2.0 * (K * K - 1.0) / a0;
 127     pa[2] = (1.0 - K / Q + K * K) / a0;
 128
 129     f0 = 38.13547087602444;
 130     Q = 0.5003270373238773;
 131     K = tan(M_PI * f0 / (double) st->samplerate);
 132
 133     ra[1] = 2.0 * (K * K - 1.0) / (1.0 + K / Q + K * K);
 134     ra[2] = (1.0 - K / Q + K * K) / (1.0 + K / Q + K * K);
 135
 136     st->d->b[0] = pb[0] * rb[0];
 137     st->d->b[1] = pb[0] * rb[1] + pb[1] * rb[0];
 138     st->d->b[2] = pb[0] * rb[2] + pb[1] * rb[1] + pb[2] * rb[0];
 139     st->d->b[3] = pb[1] * rb[2] + pb[2] * rb[1];
 140     st->d->b[4] = pb[2] * rb[2];
 141
 142     st->d->a[0] = pa[0] * ra[0];
 143     st->d->a[1] = pa[0] * ra[1] + pa[1] * ra[0];
 144     st->d->a[2] = pa[0] * ra[2] + pa[1] * ra[1] + pa[2] * ra[0];
 145     st->d->a[3] = pa[1] * ra[2] + pa[2] * ra[1];
 146     st->d->a[4] = pa[2] * ra[2];
 147
 148     for (i = 0; i < 5; ++i) {
 149         for (j = 0; j < 5; ++j) {
 150             st->d->v[i][j] = 0.0;
 151         }
 152     }
 153 }
 154
 155 static int ebur128_init_channel_map(FFEBUR128State * st)
 156 {
 157     size_t i;
 158     st->d->channel_map =
 159         (int *) av_malloc_array(st->channels, sizeof(*st->d->channel_map));
 160     if (!st->d->channel_map)
 161         return AVERROR(ENOMEM);
 162     if (st->channels == 4) {
 163         st->d->channel_map[0] = FF_EBUR128_LEFT;
 164         st->d->channel_map[1] = FF_EBUR128_RIGHT;
 165         st->d->channel_map[2] = FF_EBUR128_LEFT_SURROUND;
 166         st->d->channel_map[3] = FF_EBUR128_RIGHT_SURROUND;
 167     } else if (st->channels == 5) {
 168         st->d->channel_map[0] = FF_EBUR128_LEFT;
 169         st->d->channel_map[1] = FF_EBUR128_RIGHT;
 170         st->d->channel_map[2] = FF_EBUR128_CENTER;
 171         st->d->channel_map[3] = FF_EBUR128_LEFT_SURROUND;
 172         st->d->channel_map[4] = FF_EBUR128_RIGHT_SURROUND;
 173     } else {
 174         for (i = 0; i < st->channels; ++i) {
 175             switch (i) {
 176             case 0:
 177                 st->d->channel_map[i] = FF_EBUR128_LEFT;
 178                 break;
 179             case 1:
 180                 st->d->channel_map[i] = FF_EBUR128_RIGHT;
 181                 break;
 182             case 2:
 183                 st->d->channel_map[i] = FF_EBUR128_CENTER;
 184                 break;
 185             case 3:
 186                 st->d->channel_map[i] = FF_EBUR128_UNUSED;
 187                 break;
 188             case 4:
 189                 st->d->channel_map[i] = FF_EBUR128_LEFT_SURROUND;
 190                 break;
 191             case 5:
 192                 st->d->channel_map[i] = FF_EBUR128_RIGHT_SURROUND;
 193                 break;
 194             default:
 195                 st->d->channel_map[i] = FF_EBUR128_UNUSED;
 196                 break;
 197             }
 198         }
 199     }
 200     return 0;
 201 }
 202
 203 static inline void init_histogram(void)
 204 {
 205     int i;
 206     /* initialize static constants */
 207     histogram_energy_boundaries[0] = pow(10.0, (-70.0 + 0.691) / 10.0);
 208     for (i = 0; i < 1000; ++i) {
 209         histogram_energies[i] =
 210             pow(10.0, ((double) i / 10.0 - 69.95 + 0.691) / 10.0);
 211     }
 212     for (i = 1; i < 1001; ++i) {
 213         histogram_energy_boundaries[i] =
 214             pow(10.0, ((double) i / 10.0 - 70.0 + 0.691) / 10.0);
 215     }
 216 }
 217
 218 FFEBUR128State *ff_ebur128_init(unsigned int channels,
 219                                 unsigned long samplerate,
 220                                 unsigned long window, int mode)
 221 {
 222     int errcode;
 223     FFEBUR128State *st;
 224
 225     st = (FFEBUR128State *) av_malloc(sizeof(*st));
 226     CHECK_ERROR(!st, 0, exit)
 227     st->d = (struct FFEBUR128StateInternal *)
 228         av_malloc(sizeof(*st->d));
 229     CHECK_ERROR(!st->d, 0, free_state)
 230     st->channels = channels;
 231     errcode = ebur128_init_channel_map(st);
 232     CHECK_ERROR(errcode, 0, free_internal)
 233
 234     st->d->sample_peak =
 235         (double *) av_mallocz_array(channels, sizeof(*st->d->sample_peak));
 236     CHECK_ERROR(!st->d->sample_peak, 0, free_channel_map)
 237
 238     st->samplerate = samplerate;
 239     st->d->samples_in_100ms = (st->samplerate + 5) / 10;
 240     st->mode = mode;
 241     if ((mode & FF_EBUR128_MODE_S) == FF_EBUR128_MODE_S) {
 242         st->d->window = FFMAX(window, 3000);
 243     } else if ((mode & FF_EBUR128_MODE_M) == FF_EBUR128_MODE_M) {
 244         st->d->window = FFMAX(window, 400);
 245     } else {
 246         goto free_sample_peak;
 247     }
 248     st->d->audio_data_frames = st->samplerate * st->d->window / 1000;
 249     if (st->d->audio_data_frames % st->d->samples_in_100ms) {
 250         /* round up to multiple of samples_in_100ms */
 251         st->d->audio_data_frames = st->d->audio_data_frames
 252             + st->d->samples_in_100ms
 253             - (st->d->audio_data_frames % st->d->samples_in_100ms);
 254     }
 255     st->d->audio_data =
 256         (double *) av_mallocz_array(st->d->audio_data_frames,
 257                                     st->channels * sizeof(*st->d->audio_data));
 258     CHECK_ERROR(!st->d->audio_data, 0, free_sample_peak)
 259
 260     ebur128_init_filter(st);
 261
 262     st->d->block_energy_histogram =
 263         av_mallocz(1000 * sizeof(*st->d->block_energy_histogram));
 264     CHECK_ERROR(!st->d->block_energy_histogram, 0, free_audio_data)
 265     st->d->short_term_block_energy_histogram =
 266         av_mallocz(1000 * sizeof(*st->d->short_term_block_energy_histogram));
 267     CHECK_ERROR(!st->d->short_term_block_energy_histogram, 0,
 268                 free_block_energy_histogram)
 269     st->d->short_term_frame_counter = 0;
 270
 271     /* the first block needs 400ms of audio data */
 272     st->d->needed_frames = st->d->samples_in_100ms * 4;
 273     /* start at the beginning of the buffer */
 274     st->d->audio_data_index = 0;
 275
 276     if (ff_thread_once(&histogram_init, &init_histogram) != 0)
 277         goto free_short_term_block_energy_histogram;
 278
 279     st->d->data_ptrs = av_malloc_array(channels, sizeof(*st->d->data_ptrs));
 280     CHECK_ERROR(!st->d->data_ptrs, 0,
 281                 free_short_term_block_energy_histogram);
 282
 283     return st;
 284
 285 free_short_term_block_energy_histogram:
 286     av_free(st->d->short_term_block_energy_histogram);
 287 free_block_energy_histogram:
 288     av_free(st->d->block_energy_histogram);
 289 free_audio_data:
 290     av_free(st->d->audio_data);
 291 free_sample_peak:
 292     av_free(st->d->sample_peak);
 293 free_channel_map:
 294     av_free(st->d->channel_map);
 295 free_internal:
 296     av_free(st->d);
 297 free_state:
 298     av_free(st);
 299 exit:
 300     return NULL;
 301 }
 302
 303 void ff_ebur128_destroy(FFEBUR128State ** st)
 304 {
 305     av_free((*st)->d->block_energy_histogram);
 306     av_free((*st)->d->short_term_block_energy_histogram);
 307     av_free((*st)->d->audio_data);
 308     av_free((*st)->d->channel_map);
 309     av_free((*st)->d->sample_peak);
 310     av_free((*st)->d->data_ptrs);
 311     av_free((*st)->d);
 312     av_free(*st);
 313     *st = NULL;
 314 }
 315
 316 #define EBUR128_FILTER(type, scaling_factor)                                       \
 317 static void ebur128_filter_##type(FFEBUR128State* st, const type** srcs,           \
 318                                   size_t src_index, size_t frames,                 \
 319                                   int stride) {                                    \
 320     double* audio_data = st->d->audio_data + st->d->audio_data_index;              \
 321     size_t i, c;                                                                   \
 322                                                                                    \
 323     if ((st->mode & FF_EBUR128_MODE_SAMPLE_PEAK) == FF_EBUR128_MODE_SAMPLE_PEAK) { \
 324         for (c = 0; c < st->channels; ++c) {                                       \
 325             double max = 0.0;                                                      \
 326             for (i = 0; i < frames; ++i) {                                         \
 327                 type v = srcs[c][src_index + i * stride];                          \
 328                 if (v > max) {                                                     \
 329                     max =        v;                                                \
 330                 } else if (-v > max) {                                             \
 331                     max = -1.0 * v;                                                \
 332                 }                                                                  \
 333             }                                                                      \
 334             max /= scaling_factor;                                                 \
 335             if (max > st->d->sample_peak[c]) st->d->sample_peak[c] = max;          \
 336         }                                                                          \
 337     }                                                                              \
 338     for (c = 0; c < st->channels; ++c) {                                           \
 339         int ci = st->d->channel_map[c] - 1;                                        \
 340         if (ci < 0) continue;                                                      \
 341         else if (ci == FF_EBUR128_DUAL_MONO - 1) ci = 0; /*dual mono */            \
 342         for (i = 0; i < frames; ++i) {                                             \
 343             st->d->v[ci][0] = (double) (srcs[c][src_index + i * stride] / scaling_factor) \
 344                          - st->d->a[1] * st->d->v[ci][1]                           \
 345                          - st->d->a[2] * st->d->v[ci][2]                           \
 346                          - st->d->a[3] * st->d->v[ci][3]                           \
 347                          - st->d->a[4] * st->d->v[ci][4];                          \
 348             audio_data[i * st->channels + c] =                                     \
 349                            st->d->b[0] * st->d->v[ci][0]                           \
 350                          + st->d->b[1] * st->d->v[ci][1]                           \
 351                          + st->d->b[2] * st->d->v[ci][2]                           \
 352                          + st->d->b[3] * st->d->v[ci][3]                           \
 353                          + st->d->b[4] * st->d->v[ci][4];                          \
 354             st->d->v[ci][4] = st->d->v[ci][3];                                     \
 355             st->d->v[ci][3] = st->d->v[ci][2];                                     \
 356             st->d->v[ci][2] = st->d->v[ci][1];                                     \
 357             st->d->v[ci][1] = st->d->v[ci][0];                                     \
 358         }                                                                          \
 359         st->d->v[ci][4] = fabs(st->d->v[ci][4]) < DBL_MIN ? 0.0 : st->d->v[ci][4]; \
 360         st->d->v[ci][3] = fabs(st->d->v[ci][3]) < DBL_MIN ? 0.0 : st->d->v[ci][3]; \
 361         st->d->v[ci][2] = fabs(st->d->v[ci][2]) < DBL_MIN ? 0.0 : st->d->v[ci][2]; \
 362         st->d->v[ci][1] = fabs(st->d->v[ci][1]) < DBL_MIN ? 0.0 : st->d->v[ci][1]; \
 363     }                                                                              \
 364 }
 365 EBUR128_FILTER(short, -((double)SHRT_MIN))
 366 EBUR128_FILTER(int, -((double)INT_MIN))
 367 EBUR128_FILTER(float,  1.0)
 368 EBUR128_FILTER(double, 1.0)
 369
 370 static double ebur128_energy_to_loudness(double energy)
 371 {
 372     return 10 * log10(energy) - 0.691;
 373 }
 374
 375 static size_t find_histogram_index(double energy)
 376 {
 377     size_t index_min = 0;
 378     size_t index_max = 1000;
 379     size_t index_mid;
 380
 381     do {
 382         index_mid = (index_min + index_max) / 2;
 383         if (energy >= histogram_energy_boundaries[index_mid]) {
 384             index_min = index_mid;
 385         } else {
 386             index_max = index_mid;
 387         }
 388     } while (index_max - index_min != 1);
 389
 390     return index_min;
 391 }
 392
 393 static void ebur128_calc_gating_block(FFEBUR128State * st,
 394                                       size_t frames_per_block,
 395                                       double *optional_output)
 396 {
 397     size_t i, c;
 398     double sum = 0.0;
 399     double channel_sum;
 400     for (c = 0; c < st->channels; ++c) {
 401         if (st->d->channel_map[c] == FF_EBUR128_UNUSED)
 402             continue;
 403         channel_sum = 0.0;
 404         if (st->d->audio_data_index < frames_per_block * st->channels) {
 405             for (i = 0; i < st->d->audio_data_index / st->channels; ++i) {
 406                 channel_sum += st->d->audio_data[i * st->channels + c] *
 407                     st->d->audio_data[i * st->channels + c];
 408             }
 409             for (i = st->d->audio_data_frames -
 410                  (frames_per_block -
 411                   st->d->audio_data_index / st->channels);
 412                  i < st->d->audio_data_frames; ++i) {
 413                 channel_sum += st->d->audio_data[i * st->channels + c] *
 414                     st->d->audio_data[i * st->channels + c];
 415             }
 416         } else {
 417             for (i =
 418                  st->d->audio_data_index / st->channels - frames_per_block;
 419                  i < st->d->audio_data_index / st->channels; ++i) {
 420                 channel_sum +=
 421                     st->d->audio_data[i * st->channels +
 422                                       c] * st->d->audio_data[i *
 423                                                              st->channels +
 424                                                              c];
 425             }
 426         }
 427         if (st->d->channel_map[c] == FF_EBUR128_Mp110 ||
 428             st->d->channel_map[c] == FF_EBUR128_Mm110 ||
 429             st->d->channel_map[c] == FF_EBUR128_Mp060 ||
 430             st->d->channel_map[c] == FF_EBUR128_Mm060 ||
 431             st->d->channel_map[c] == FF_EBUR128_Mp090 ||
 432             st->d->channel_map[c] == FF_EBUR128_Mm090) {
 433             channel_sum *= 1.41;
 434         } else if (st->d->channel_map[c] == FF_EBUR128_DUAL_MONO) {
 435             channel_sum *= 2.0;
 436         }
 437         sum += channel_sum;
 438     }
 439     sum /= (double) frames_per_block;
 440     if (optional_output) {
 441         *optional_output = sum;
 442     } else if (sum >= histogram_energy_boundaries[0]) {
 443         ++st->d->block_energy_histogram[find_histogram_index(sum)];
 444     }
 445 }
 446
 447 int ff_ebur128_set_channel(FFEBUR128State * st,
 448                            unsigned int channel_number, int value)
 449 {
 450     if (channel_number >= st->channels) {
 451         return 1;
 452     }
 453     if (value == FF_EBUR128_DUAL_MONO &&
 454         (st->channels != 1 || channel_number != 0)) {
 455         return 1;
 456     }
 457     st->d->channel_map[channel_number] = value;
 458     return 0;
 459 }
 460
 461 static int ebur128_energy_shortterm(FFEBUR128State * st, double *out);
 462 #define FF_EBUR128_ADD_FRAMES_PLANAR(type)                                             \
 463 void ff_ebur128_add_frames_planar_##type(FFEBUR128State* st, const type** srcs,        \
 464                                  size_t frames, int stride) {                          \
 465     size_t src_index = 0;                                                              \
 466     while (frames > 0) {                                                               \
 467         if (frames >= st->d->needed_frames) {                                          \
 468             ebur128_filter_##type(st, srcs, src_index, st->d->needed_frames, stride);  \
 469             src_index += st->d->needed_frames * stride;                                \
 470             frames -= st->d->needed_frames;                                            \
 471             st->d->audio_data_index += st->d->needed_frames * st->channels;            \
 472             /* calculate the new gating block */                                       \
 473             if ((st->mode & FF_EBUR128_MODE_I) == FF_EBUR128_MODE_I) {                 \
 474                 ebur128_calc_gating_block(st, st->d->samples_in_100ms * 4, NULL);      \
 475             }                                                                          \
 476             if ((st->mode & FF_EBUR128_MODE_LRA) == FF_EBUR128_MODE_LRA) {             \
 477                 st->d->short_term_frame_counter += st->d->needed_frames;               \
 478                 if (st->d->short_term_frame_counter == st->d->samples_in_100ms * 30) { \
 479                     double st_energy;                                                  \
 480                     ebur128_energy_shortterm(st, &st_energy);                          \
 481                     if (st_energy >= histogram_energy_boundaries[0]) {                 \
 482                         ++st->d->short_term_block_energy_histogram[                    \
 483                                                     find_histogram_index(st_energy)];  \
 484                     }                                                                  \
 485                     st->d->short_term_frame_counter = st->d->samples_in_100ms * 20;    \
 486                 }                                                                      \
 487             }                                                                          \
 488             /* 100ms are needed for all blocks besides the first one */                \
 489             st->d->needed_frames = st->d->samples_in_100ms;                            \
 490             /* reset audio_data_index when buffer full */                              \
 491             if (st->d->audio_data_index == st->d->audio_data_frames * st->channels) {  \
 492                 st->d->audio_data_index = 0;                                           \
 493             }                                                                          \
 494         } else {                                                                       \
 495             ebur128_filter_##type(st, srcs, src_index, frames, stride);                \
 496             st->d->audio_data_index += frames * st->channels;                          \
 497             if ((st->mode & FF_EBUR128_MODE_LRA) == FF_EBUR128_MODE_LRA) {             \
 498                 st->d->short_term_frame_counter += frames;                             \
 499             }                                                                          \
 500             st->d->needed_frames -= frames;                                            \
 501             frames = 0;                                                                \
 502         }                                                                              \
 503     }                                                                                  \
 504 }
 505 FF_EBUR128_ADD_FRAMES_PLANAR(short)
 506 FF_EBUR128_ADD_FRAMES_PLANAR(int)
 507 FF_EBUR128_ADD_FRAMES_PLANAR(float)
 508 FF_EBUR128_ADD_FRAMES_PLANAR(double)
 509 #define FF_EBUR128_ADD_FRAMES(type)                                            \
 510 void ff_ebur128_add_frames_##type(FFEBUR128State* st, const type* src,         \
 511                                     size_t frames) {                           \
 512   int i;                                                                       \
 513   const type **buf = (const type**)st->d->data_ptrs;                           \
 514   for (i = 0; i < st->channels; i++)                                           \
 515     buf[i] = src + i;                                                          \
 516   ff_ebur128_add_frames_planar_##type(st, buf, frames, st->channels);          \
 517 }
 518 FF_EBUR128_ADD_FRAMES(short)
 519 FF_EBUR128_ADD_FRAMES(int)
 520 FF_EBUR128_ADD_FRAMES(float)
 521 FF_EBUR128_ADD_FRAMES(double)
 522
 523 static int ebur128_calc_relative_threshold(FFEBUR128State **sts, size_t size,
 524                                            double *relative_threshold)
 525 {
 526     size_t i, j;
 527     int above_thresh_counter = 0;
 528     *relative_threshold = 0.0;
 529
 530     for (i = 0; i < size; i++) {
 531         unsigned long *block_energy_histogram = sts[i]->d->block_energy_histogram;
 532         for (j = 0; j < 1000; ++j) {
 533             *relative_threshold += block_energy_histogram[j] * histogram_energies[j];
 534             above_thresh_counter += block_energy_histogram[j];
 535         }
 536     }
 537
 538     if (above_thresh_counter != 0) {
 539         *relative_threshold /= (double)above_thresh_counter;
 540         *relative_threshold *= RELATIVE_GATE_FACTOR;
 541     }
 542
 543     return above_thresh_counter;
 544 }
 545
 546 static int ebur128_gated_loudness(FFEBUR128State ** sts, size_t size,
 547                                   double *out)
 548 {
 549     double gated_loudness = 0.0;
 550     double relative_threshold;
 551     size_t above_thresh_counter;
 552     size_t i, j, start_index;
 553
 554     for (i = 0; i < size; i++)
 555         if ((sts[i]->mode & FF_EBUR128_MODE_I) != FF_EBUR128_MODE_I)
 556             return AVERROR(EINVAL);
 557
 558     if (!ebur128_calc_relative_threshold(sts, size, &relative_threshold)) {
 559         *out = -HUGE_VAL;
 560         return 0;
 561     }
 562
 563     above_thresh_counter = 0;
 564     if (relative_threshold < histogram_energy_boundaries[0]) {
 565         start_index = 0;
 566     } else {
 567         start_index = find_histogram_index(relative_threshold);
 568         if (relative_threshold > histogram_energies[start_index]) {
 569             ++start_index;
 570         }
 571     }
 572     for (i = 0; i < size; i++) {
 573         for (j = start_index; j < 1000; ++j) {
 574             gated_loudness += sts[i]->d->block_energy_histogram[j] *
 575                 histogram_energies[j];
 576             above_thresh_counter += sts[i]->d->block_energy_histogram[j];
 577         }
 578     }
 579     if (!above_thresh_counter) {
 580         *out = -HUGE_VAL;
 581         return 0;
 582     }
 583     gated_loudness /= (double) above_thresh_counter;
 584     *out = ebur128_energy_to_loudness(gated_loudness);
 585     return 0;
 586 }
 587
 588 int ff_ebur128_relative_threshold(FFEBUR128State * st, double *out)
 589 {
 590     double relative_threshold;
 591
 592     if ((st->mode & FF_EBUR128_MODE_I) != FF_EBUR128_MODE_I)
 593         return AVERROR(EINVAL);
 594
 595     if (!ebur128_calc_relative_threshold(&st, 1, &relative_threshold)) {
 596         *out = -70.0;
 597         return 0;
 598     }
 599
 600     *out = ebur128_energy_to_loudness(relative_threshold);
 601     return 0;
 602 }
 603
 604 int ff_ebur128_loudness_global(FFEBUR128State * st, double *out)
 605 {
 606     return ebur128_gated_loudness(&st, 1, out);
 607 }
 608
 609 int ff_ebur128_loudness_global_multiple(FFEBUR128State ** sts, size_t size,
 610                                         double *out)
 611 {
 612     return ebur128_gated_loudness(sts, size, out);
 613 }
 614
 615 static int ebur128_energy_in_interval(FFEBUR128State * st,
 616                                       size_t interval_frames, double *out)
 617 {
 618     if (interval_frames > st->d->audio_data_frames) {
 619         return AVERROR(EINVAL);
 620     }
 621     ebur128_calc_gating_block(st, interval_frames, out);
 622     return 0;
 623 }
 624
 625 static int ebur128_energy_shortterm(FFEBUR128State * st, double *out)
 626 {
 627     return ebur128_energy_in_interval(st, st->d->samples_in_100ms * 30,
 628                                       out);
 629 }
 630
 631 int ff_ebur128_loudness_momentary(FFEBUR128State * st, double *out)
 632 {
 633     double energy;
 634     int error = ebur128_energy_in_interval(st, st->d->samples_in_100ms * 4,
 635                                            &energy);
 636     if (error) {
 637         return error;
 638     } else if (energy <= 0.0) {
 639         *out = -HUGE_VAL;
 640         return 0;
 641     }
 642     *out = ebur128_energy_to_loudness(energy);
 643     return 0;
 644 }
 645
 646 int ff_ebur128_loudness_shortterm(FFEBUR128State * st, double *out)
 647 {
 648     double energy;
 649     int error = ebur128_energy_shortterm(st, &energy);
 650     if (error) {
 651         return error;
 652     } else if (energy <= 0.0) {
 653         *out = -HUGE_VAL;
 654         return 0;
 655     }
 656     *out = ebur128_energy_to_loudness(energy);
 657     return 0;
 658 }
 659
 660 int ff_ebur128_loudness_window(FFEBUR128State * st,
 661                                unsigned long window, double *out)
 662 {
 663     double energy;
 664     size_t interval_frames = st->samplerate * window / 1000;
 665     int error = ebur128_energy_in_interval(st, interval_frames, &energy);
 666     if (error) {
 667         return error;
 668     } else if (energy <= 0.0) {
 669         *out = -HUGE_VAL;
 670         return 0;
 671     }
 672     *out = ebur128_energy_to_loudness(energy);
 673     return 0;
 674 }
 675
 676 /* EBU - TECH 3342 */
 677 int ff_ebur128_loudness_range_multiple(FFEBUR128State ** sts, size_t size,
 678                                        double *out)
 679 {
 680     size_t i, j;
 681     size_t stl_size;
 682     double stl_power, stl_integrated;
 683     /* High and low percentile energy */
 684     double h_en, l_en;
 685     unsigned long hist[1000] = { 0 };
 686     size_t percentile_low, percentile_high;
 687     size_t index;
 688
 689     for (i = 0; i < size; ++i) {
 690         if (sts[i]) {
 691             if ((sts[i]->mode & FF_EBUR128_MODE_LRA) !=
 692                 FF_EBUR128_MODE_LRA) {
 693                 return AVERROR(EINVAL);
 694             }
 695         }
 696     }
 697
 698     stl_size = 0;
 699     stl_power = 0.0;
 700     for (i = 0; i < size; ++i) {
 701         if (!sts[i])
 702             continue;
 703         for (j = 0; j < 1000; ++j) {
 704             hist[j] += sts[i]->d->short_term_block_energy_histogram[j];
 705             stl_size += sts[i]->d->short_term_block_energy_histogram[j];
 706             stl_power += sts[i]->d->short_term_block_energy_histogram[j]
 707                 * histogram_energies[j];
 708         }
 709     }
 710     if (!stl_size) {
 711         *out = 0.0;
 712         return 0;
 713     }
 714
 715     stl_power /= stl_size;
 716     stl_integrated = MINUS_20DB * stl_power;
 717
 718     if (stl_integrated < histogram_energy_boundaries[0]) {
 719         index = 0;
 720     } else {
 721         index = find_histogram_index(stl_integrated);
 722         if (stl_integrated > histogram_energies[index]) {
 723             ++index;
 724         }
 725     }
 726     stl_size = 0;
 727     for (j = index; j < 1000; ++j) {
 728         stl_size += hist[j];
 729     }
 730     if (!stl_size) {
 731         *out = 0.0;
 732         return 0;
 733     }
 734
 735     percentile_low = (size_t) ((stl_size - 1) * 0.1 + 0.5);
 736     percentile_high = (size_t) ((stl_size - 1) * 0.95 + 0.5);
 737
 738     stl_size = 0;
 739     j = index;
 740     while (stl_size <= percentile_low) {
 741         stl_size += hist[j++];
 742     }
 743     l_en = histogram_energies[j - 1];
 744     while (stl_size <= percentile_high) {
 745         stl_size += hist[j++];
 746     }
 747     h_en = histogram_energies[j - 1];
 748     *out =
 749         ebur128_energy_to_loudness(h_en) -
 750         ebur128_energy_to_loudness(l_en);
 751     return 0;
 752 }
 753
 754 int ff_ebur128_loudness_range(FFEBUR128State * st, double *out)
 755 {
 756     return ff_ebur128_loudness_range_multiple(&st, 1, out);
 757 }
 758
 759 int ff_ebur128_sample_peak(FFEBUR128State * st,
 760                            unsigned int channel_number, double *out)
 761 {
 762     if ((st->mode & FF_EBUR128_MODE_SAMPLE_PEAK) !=
 763         FF_EBUR128_MODE_SAMPLE_PEAK) {
 764         return AVERROR(EINVAL);
 765     } else if (channel_number >= st->channels) {
 766         return AVERROR(EINVAL);
 767     }
 768     *out = st->d->sample_peak[channel_number];
 769     return 0;
 770 }