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(double, 1.0)
 366
 367 static double ebur128_energy_to_loudness(double energy)
 368 {
 369     return 10 * log10(energy) - 0.691;
 370 }
 371
 372 static size_t find_histogram_index(double energy)
 373 {
 374     size_t index_min = 0;
 375     size_t index_max = 1000;
 376     size_t index_mid;
 377
 378     do {
 379         index_mid = (index_min + index_max) / 2;
 380         if (energy >= histogram_energy_boundaries[index_mid]) {
 381             index_min = index_mid;
 382         } else {
 383             index_max = index_mid;
 384         }
 385     } while (index_max - index_min != 1);
 386
 387     return index_min;
 388 }
 389
 390 static void ebur128_calc_gating_block(FFEBUR128State * st,
 391                                       size_t frames_per_block,
 392                                       double *optional_output)
 393 {
 394     size_t i, c;
 395     double sum = 0.0;
 396     double channel_sum;
 397     for (c = 0; c < st->channels; ++c) {
 398         if (st->d->channel_map[c] == FF_EBUR128_UNUSED)
 399             continue;
 400         channel_sum = 0.0;
 401         if (st->d->audio_data_index < frames_per_block * st->channels) {
 402             for (i = 0; i < st->d->audio_data_index / st->channels; ++i) {
 403                 channel_sum += st->d->audio_data[i * st->channels + c] *
 404                     st->d->audio_data[i * st->channels + c];
 405             }
 406             for (i = st->d->audio_data_frames -
 407                  (frames_per_block -
 408                   st->d->audio_data_index / st->channels);
 409                  i < st->d->audio_data_frames; ++i) {
 410                 channel_sum += st->d->audio_data[i * st->channels + c] *
 411                     st->d->audio_data[i * st->channels + c];
 412             }
 413         } else {
 414             for (i =
 415                  st->d->audio_data_index / st->channels - frames_per_block;
 416                  i < st->d->audio_data_index / st->channels; ++i) {
 417                 channel_sum +=
 418                     st->d->audio_data[i * st->channels +
 419                                       c] * st->d->audio_data[i *
 420                                                              st->channels +
 421                                                              c];
 422             }
 423         }
 424         if (st->d->channel_map[c] == FF_EBUR128_Mp110 ||
 425             st->d->channel_map[c] == FF_EBUR128_Mm110 ||
 426             st->d->channel_map[c] == FF_EBUR128_Mp060 ||
 427             st->d->channel_map[c] == FF_EBUR128_Mm060 ||
 428             st->d->channel_map[c] == FF_EBUR128_Mp090 ||
 429             st->d->channel_map[c] == FF_EBUR128_Mm090) {
 430             channel_sum *= 1.41;
 431         } else if (st->d->channel_map[c] == FF_EBUR128_DUAL_MONO) {
 432             channel_sum *= 2.0;
 433         }
 434         sum += channel_sum;
 435     }
 436     sum /= (double) frames_per_block;
 437     if (optional_output) {
 438         *optional_output = sum;
 439     } else if (sum >= histogram_energy_boundaries[0]) {
 440         ++st->d->block_energy_histogram[find_histogram_index(sum)];
 441     }
 442 }
 443
 444 int ff_ebur128_set_channel(FFEBUR128State * st,
 445                            unsigned int channel_number, int value)
 446 {
 447     if (channel_number >= st->channels) {
 448         return 1;
 449     }
 450     if (value == FF_EBUR128_DUAL_MONO &&
 451         (st->channels != 1 || channel_number != 0)) {
 452         return 1;
 453     }
 454     st->d->channel_map[channel_number] = value;
 455     return 0;
 456 }
 457
 458 static int ebur128_energy_shortterm(FFEBUR128State * st, double *out);
 459 #define EBUR128_ADD_FRAMES_PLANAR(type)                                                \
 460 static void ebur128_add_frames_planar_##type(FFEBUR128State* st, const type** srcs,    \
 461                                  size_t frames, int stride) {                          \
 462     size_t src_index = 0;                                                              \
 463     while (frames > 0) {                                                               \
 464         if (frames >= st->d->needed_frames) {                                          \
 465             ebur128_filter_##type(st, srcs, src_index, st->d->needed_frames, stride);  \
 466             src_index += st->d->needed_frames * stride;                                \
 467             frames -= st->d->needed_frames;                                            \
 468             st->d->audio_data_index += st->d->needed_frames * st->channels;            \
 469             /* calculate the new gating block */                                       \
 470             if ((st->mode & FF_EBUR128_MODE_I) == FF_EBUR128_MODE_I) {                 \
 471                 ebur128_calc_gating_block(st, st->d->samples_in_100ms * 4, NULL);      \
 472             }                                                                          \
 473             if ((st->mode & FF_EBUR128_MODE_LRA) == FF_EBUR128_MODE_LRA) {             \
 474                 st->d->short_term_frame_counter += st->d->needed_frames;               \
 475                 if (st->d->short_term_frame_counter == st->d->samples_in_100ms * 30) { \
 476                     double st_energy;                                                  \
 477                     ebur128_energy_shortterm(st, &st_energy);                          \
 478                     if (st_energy >= histogram_energy_boundaries[0]) {                 \
 479                         ++st->d->short_term_block_energy_histogram[                    \
 480                                                     find_histogram_index(st_energy)];  \
 481                     }                                                                  \
 482                     st->d->short_term_frame_counter = st->d->samples_in_100ms * 20;    \
 483                 }                                                                      \
 484             }                                                                          \
 485             /* 100ms are needed for all blocks besides the first one */                \
 486             st->d->needed_frames = st->d->samples_in_100ms;                            \
 487             /* reset audio_data_index when buffer full */                              \
 488             if (st->d->audio_data_index == st->d->audio_data_frames * st->channels) {  \
 489                 st->d->audio_data_index = 0;                                           \
 490             }                                                                          \
 491         } else {                                                                       \
 492             ebur128_filter_##type(st, srcs, src_index, frames, stride);                \
 493             st->d->audio_data_index += frames * st->channels;                          \
 494             if ((st->mode & FF_EBUR128_MODE_LRA) == FF_EBUR128_MODE_LRA) {             \
 495                 st->d->short_term_frame_counter += frames;                             \
 496             }                                                                          \
 497             st->d->needed_frames -= frames;                                            \
 498             frames = 0;                                                                \
 499         }                                                                              \
 500     }                                                                                  \
 501 }
 502 EBUR128_ADD_FRAMES_PLANAR(double)
 503 #define FF_EBUR128_ADD_FRAMES(type)                                            \
 504 void ff_ebur128_add_frames_##type(FFEBUR128State* st, const type* src,         \
 505                                     size_t frames) {                           \
 506   int i;                                                                       \
 507   const type **buf = (const type**)st->d->data_ptrs;                           \
 508   for (i = 0; i < st->channels; i++)                                           \
 509     buf[i] = src + i;                                                          \
 510   ebur128_add_frames_planar_##type(st, buf, frames, st->channels);             \
 511 }
 512 FF_EBUR128_ADD_FRAMES(double)
 513
 514 static int ebur128_calc_relative_threshold(FFEBUR128State **sts, size_t size,
 515                                            double *relative_threshold)
 516 {
 517     size_t i, j;
 518     int above_thresh_counter = 0;
 519     *relative_threshold = 0.0;
 520
 521     for (i = 0; i < size; i++) {
 522         unsigned long *block_energy_histogram = sts[i]->d->block_energy_histogram;
 523         for (j = 0; j < 1000; ++j) {
 524             *relative_threshold += block_energy_histogram[j] * histogram_energies[j];
 525             above_thresh_counter += block_energy_histogram[j];
 526         }
 527     }
 528
 529     if (above_thresh_counter != 0) {
 530         *relative_threshold /= (double)above_thresh_counter;
 531         *relative_threshold *= RELATIVE_GATE_FACTOR;
 532     }
 533
 534     return above_thresh_counter;
 535 }
 536
 537 static int ebur128_gated_loudness(FFEBUR128State ** sts, size_t size,
 538                                   double *out)
 539 {
 540     double gated_loudness = 0.0;
 541     double relative_threshold;
 542     size_t above_thresh_counter;
 543     size_t i, j, start_index;
 544
 545     for (i = 0; i < size; i++)
 546         if ((sts[i]->mode & FF_EBUR128_MODE_I) != FF_EBUR128_MODE_I)
 547             return AVERROR(EINVAL);
 548
 549     if (!ebur128_calc_relative_threshold(sts, size, &relative_threshold)) {
 550         *out = -HUGE_VAL;
 551         return 0;
 552     }
 553
 554     above_thresh_counter = 0;
 555     if (relative_threshold < histogram_energy_boundaries[0]) {
 556         start_index = 0;
 557     } else {
 558         start_index = find_histogram_index(relative_threshold);
 559         if (relative_threshold > histogram_energies[start_index]) {
 560             ++start_index;
 561         }
 562     }
 563     for (i = 0; i < size; i++) {
 564         for (j = start_index; j < 1000; ++j) {
 565             gated_loudness += sts[i]->d->block_energy_histogram[j] *
 566                 histogram_energies[j];
 567             above_thresh_counter += sts[i]->d->block_energy_histogram[j];
 568         }
 569     }
 570     if (!above_thresh_counter) {
 571         *out = -HUGE_VAL;
 572         return 0;
 573     }
 574     gated_loudness /= (double) above_thresh_counter;
 575     *out = ebur128_energy_to_loudness(gated_loudness);
 576     return 0;
 577 }
 578
 579 int ff_ebur128_relative_threshold(FFEBUR128State * st, double *out)
 580 {
 581     double relative_threshold;
 582
 583     if ((st->mode & FF_EBUR128_MODE_I) != FF_EBUR128_MODE_I)
 584         return AVERROR(EINVAL);
 585
 586     if (!ebur128_calc_relative_threshold(&st, 1, &relative_threshold)) {
 587         *out = -70.0;
 588         return 0;
 589     }
 590
 591     *out = ebur128_energy_to_loudness(relative_threshold);
 592     return 0;
 593 }
 594
 595 int ff_ebur128_loudness_global(FFEBUR128State * st, double *out)
 596 {
 597     return ebur128_gated_loudness(&st, 1, out);
 598 }
 599
 600 static int ebur128_energy_in_interval(FFEBUR128State * st,
 601                                       size_t interval_frames, double *out)
 602 {
 603     if (interval_frames > st->d->audio_data_frames) {
 604         return AVERROR(EINVAL);
 605     }
 606     ebur128_calc_gating_block(st, interval_frames, out);
 607     return 0;
 608 }
 609
 610 static int ebur128_energy_shortterm(FFEBUR128State * st, double *out)
 611 {
 612     return ebur128_energy_in_interval(st, st->d->samples_in_100ms * 30,
 613                                       out);
 614 }
 615
 616 int ff_ebur128_loudness_shortterm(FFEBUR128State * st, double *out)
 617 {
 618     double energy;
 619     int error = ebur128_energy_shortterm(st, &energy);
 620     if (error) {
 621         return error;
 622     } else if (energy <= 0.0) {
 623         *out = -HUGE_VAL;
 624         return 0;
 625     }
 626     *out = ebur128_energy_to_loudness(energy);
 627     return 0;
 628 }
 629
 630 /* EBU - TECH 3342 */
 631 int ff_ebur128_loudness_range_multiple(FFEBUR128State ** sts, size_t size,
 632                                        double *out)
 633 {
 634     size_t i, j;
 635     size_t stl_size;
 636     double stl_power, stl_integrated;
 637     /* High and low percentile energy */
 638     double h_en, l_en;
 639     unsigned long hist[1000] = { 0 };
 640     size_t percentile_low, percentile_high;
 641     size_t index;
 642
 643     for (i = 0; i < size; ++i) {
 644         if (sts[i]) {
 645             if ((sts[i]->mode & FF_EBUR128_MODE_LRA) !=
 646                 FF_EBUR128_MODE_LRA) {
 647                 return AVERROR(EINVAL);
 648             }
 649         }
 650     }
 651
 652     stl_size = 0;
 653     stl_power = 0.0;
 654     for (i = 0; i < size; ++i) {
 655         if (!sts[i])
 656             continue;
 657         for (j = 0; j < 1000; ++j) {
 658             hist[j] += sts[i]->d->short_term_block_energy_histogram[j];
 659             stl_size += sts[i]->d->short_term_block_energy_histogram[j];
 660             stl_power += sts[i]->d->short_term_block_energy_histogram[j]
 661                 * histogram_energies[j];
 662         }
 663     }
 664     if (!stl_size) {
 665         *out = 0.0;
 666         return 0;
 667     }
 668
 669     stl_power /= stl_size;
 670     stl_integrated = MINUS_20DB * stl_power;
 671
 672     if (stl_integrated < histogram_energy_boundaries[0]) {
 673         index = 0;
 674     } else {
 675         index = find_histogram_index(stl_integrated);
 676         if (stl_integrated > histogram_energies[index]) {
 677             ++index;
 678         }
 679     }
 680     stl_size = 0;
 681     for (j = index; j < 1000; ++j) {
 682         stl_size += hist[j];
 683     }
 684     if (!stl_size) {
 685         *out = 0.0;
 686         return 0;
 687     }
 688
 689     percentile_low = (size_t) ((stl_size - 1) * 0.1 + 0.5);
 690     percentile_high = (size_t) ((stl_size - 1) * 0.95 + 0.5);
 691
 692     stl_size = 0;
 693     j = index;
 694     while (stl_size <= percentile_low) {
 695         stl_size += hist[j++];
 696     }
 697     l_en = histogram_energies[j - 1];
 698     while (stl_size <= percentile_high) {
 699         stl_size += hist[j++];
 700     }
 701     h_en = histogram_energies[j - 1];
 702     *out =
 703         ebur128_energy_to_loudness(h_en) -
 704         ebur128_energy_to_loudness(l_en);
 705     return 0;
 706 }
 707
 708 int ff_ebur128_loudness_range(FFEBUR128State * st, double *out)
 709 {
 710     return ff_ebur128_loudness_range_multiple(&st, 1, out);
 711 }
 712
 713 int ff_ebur128_sample_peak(FFEBUR128State * st,
 714                            unsigned int channel_number, double *out)
 715 {
 716     if ((st->mode & FF_EBUR128_MODE_SAMPLE_PEAK) !=
 717         FF_EBUR128_MODE_SAMPLE_PEAK) {
 718         return AVERROR(EINVAL);
 719     } else if (channel_number >= st->channels) {
 720         return AVERROR(EINVAL);
 721     }
 722     *out = st->d->sample_peak[channel_number];
 723     return 0;
 724 }