git.sesse.net Git - mlt/blob - src/modules/plus/ebur128/ebur128.c

   1 /* See COPYING file for copyright and license details. */
   2
   3 #include "ebur128.h"
   4
   5 #include <float.h>
   6 #include <limits.h>
   7 #include <math.h> /* You may have to define _USE_MATH_DEFINES if you use MSVC */
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10
  11 /* This can be replaced by any BSD-like queue implementation. */
  12 #include "queue.h"
  13
  14 #ifdef USE_SPEEX_RESAMPLER
  15   #include <speex/speex_resampler.h>
  16 #endif
  17
  18 #define CHECK_ERROR(condition, errorcode, goto_point)                          \
  19   if ((condition)) {                                                           \
  20     errcode = (errorcode);                                                     \
  21     goto goto_point;                                                           \
  22   }
  23
  24 SLIST_HEAD(ebur128_double_queue, ebur128_dq_entry);
  25 struct ebur128_dq_entry {
  26   double z;
  27   SLIST_ENTRY(ebur128_dq_entry) entries;
  28 };
  29
  30 struct ebur128_state_internal {
  31   /** Filtered audio data (used as ring buffer). */
  32   double* audio_data;
  33   /** Size of audio_data array. */
  34   size_t audio_data_frames;
  35   /** Current index for audio_data. */
  36   size_t audio_data_index;
  37   /** How many frames are needed for a gating block. Will correspond to 400ms
  38    *  of audio at initialization, and 100ms after the first block (75% overlap
  39    *  as specified in the 2011 revision of BS1770). */
  40   unsigned long needed_frames;
  41   /** The channel map. Has as many elements as there are channels. */
  42   int* channel_map;
  43   /** How many samples fit in 100ms (rounded). */
  44   unsigned long samples_in_100ms;
  45   /** BS.1770 filter coefficients (nominator). */
  46   double b[5];
  47   /** BS.1770 filter coefficients (denominator). */
  48   double a[5];
  49   /** BS.1770 filter state. */
  50   double v[5][5];
  51   /** Linked list of block energies. */
  52   struct ebur128_double_queue block_list;
  53   /** Linked list of 3s-block energies, used to calculate LRA. */
  54   struct ebur128_double_queue short_term_block_list;
  55   int use_histogram;
  56   unsigned long *block_energy_histogram;
  57   unsigned long *short_term_block_energy_histogram;
  58   /** Keeps track of when a new short term block is needed. */
  59   size_t short_term_frame_counter;
  60   /** Maximum sample peak, one per channel */
  61   double* sample_peak;
  62   /** Maximum true peak, one per channel */
  63   double* true_peak;
  64 #ifdef USE_SPEEX_RESAMPLER
  65   SpeexResamplerState* resampler;
  66 #endif
  67   size_t oversample_factor;
  68   float* resampler_buffer_input;
  69   size_t resampler_buffer_input_frames;
  70   float* resampler_buffer_output;
  71   size_t resampler_buffer_output_frames;
  72 };
  73
  74 static double relative_gate = -10.0;
  75
  76 /* Those will be calculated when initializing the library */
  77 static double relative_gate_factor;
  78 static double minus_twenty_decibels;
  79 static double histogram_energies[1000];
  80 static double histogram_energy_boundaries[1001];
  81
  82 static void ebur128_init_filter(ebur128_state* st) {
  83   int i, j;
  84
  85   double f0 = 1681.974450955533;
  86   double G  =    3.999843853973347;
  87   double Q  =    0.7071752369554196;
  88
  89   double K  = tan(M_PI * f0 / (double) st->samplerate);
  90   double Vh = pow(10.0, G / 20.0);
  91   double Vb = pow(Vh, 0.4996667741545416);
  92
  93   double pb[3] = {0.0,  0.0, 0.0};
  94   double pa[3] = {1.0,  0.0, 0.0};
  95   double rb[3] = {1.0, -2.0, 1.0};
  96   double ra[3] = {1.0,  0.0, 0.0};
  97
  98   double a0 =      1.0 + K / Q + K * K      ;
  99   pb[0] =     (Vh + Vb * K / Q + K * K) / a0;
 100   pb[1] =           2.0 * (K * K -  Vh) / a0;
 101   pb[2] =     (Vh - Vb * K / Q + K * K) / a0;
 102   pa[1] =           2.0 * (K * K - 1.0) / a0;
 103   pa[2] =         (1.0 - K / Q + K * K) / a0;
 104
 105   /* fprintf(stderr, "%.14f %.14f %.14f %.14f %.14f\n",
 106                      b1[0], b1[1], b1[2], a1[1], a1[2]); */
 107
 108   f0 = 38.13547087602444;
 109   Q  =  0.5003270373238773;
 110   K  = tan(M_PI * f0 / (double) st->samplerate);
 111
 112   ra[1] =   2.0 * (K * K - 1.0) / (1.0 + K / Q + K * K);
 113   ra[2] = (1.0 - K / Q + K * K) / (1.0 + K / Q + K * K);
 114
 115   /* fprintf(stderr, "%.14f %.14f\n", a2[1], a2[2]); */
 116
 117   st->d->b[0] = pb[0] * rb[0];
 118   st->d->b[1] = pb[0] * rb[1] + pb[1] * rb[0];
 119   st->d->b[2] = pb[0] * rb[2] + pb[1] * rb[1] + pb[2] * rb[0];
 120   st->d->b[3] = pb[1] * rb[2] + pb[2] * rb[1];
 121   st->d->b[4] = pb[2] * rb[2];
 122
 123   st->d->a[0] = pa[0] * ra[0];
 124   st->d->a[1] = pa[0] * ra[1] + pa[1] * ra[0];
 125   st->d->a[2] = pa[0] * ra[2] + pa[1] * ra[1] + pa[2] * ra[0];
 126   st->d->a[3] = pa[1] * ra[2] + pa[2] * ra[1];
 127   st->d->a[4] = pa[2] * ra[2];
 128
 129   for (i = 0; i < 5; ++i) {
 130     for (j = 0; j < 5; ++j) {
 131       st->d->v[i][j] = 0.0;
 132     }
 133   }
 134 }
 135
 136 static int ebur128_init_channel_map(ebur128_state* st) {
 137   size_t i;
 138   st->d->channel_map = (int*) malloc(st->channels * sizeof(int));
 139   if (!st->d->channel_map) return EBUR128_ERROR_NOMEM;
 140   if (st->channels == 4) {
 141     st->d->channel_map[0] = EBUR128_LEFT;
 142     st->d->channel_map[1] = EBUR128_RIGHT;
 143     st->d->channel_map[2] = EBUR128_LEFT_SURROUND;
 144     st->d->channel_map[3] = EBUR128_RIGHT_SURROUND;
 145   } else if (st->channels == 5) {
 146     st->d->channel_map[0] = EBUR128_LEFT;
 147     st->d->channel_map[1] = EBUR128_RIGHT;
 148     st->d->channel_map[2] = EBUR128_CENTER;
 149     st->d->channel_map[3] = EBUR128_LEFT_SURROUND;
 150     st->d->channel_map[4] = EBUR128_RIGHT_SURROUND;
 151   } else {
 152     for (i = 0; i < st->channels; ++i) {
 153       switch (i) {
 154         case 0:  st->d->channel_map[i] = EBUR128_LEFT;           break;
 155         case 1:  st->d->channel_map[i] = EBUR128_RIGHT;          break;
 156         case 2:  st->d->channel_map[i] = EBUR128_CENTER;         break;
 157         case 3:  st->d->channel_map[i] = EBUR128_UNUSED;         break;
 158         case 4:  st->d->channel_map[i] = EBUR128_LEFT_SURROUND;  break;
 159         case 5:  st->d->channel_map[i] = EBUR128_RIGHT_SURROUND; break;
 160         default: st->d->channel_map[i] = EBUR128_UNUSED;         break;
 161       }
 162     }
 163   }
 164   return EBUR128_SUCCESS;
 165 }
 166
 167 #ifdef USE_SPEEX_RESAMPLER
 168 static int ebur128_init_resampler(ebur128_state* st) {
 169   int errcode = EBUR128_SUCCESS;
 170
 171   if (st->samplerate < 96000) {
 172     st->d->oversample_factor = 4;
 173   } else if (st->samplerate < 192000) {
 174     st->d->oversample_factor = 2;
 175   } else {
 176     st->d->oversample_factor = 1;
 177     st->d->resampler_buffer_input = NULL;
 178     st->d->resampler_buffer_output = NULL;
 179     st->d->resampler = NULL;
 180   }
 181
 182   st->d->resampler_buffer_input_frames = st->d->samples_in_100ms * 4;
 183   st->d->resampler_buffer_input = malloc(st->d->resampler_buffer_input_frames *
 184                                       st->channels *
 185                                       sizeof(float));
 186   CHECK_ERROR(!st->d->resampler_buffer_input, EBUR128_ERROR_NOMEM, exit)
 187
 188   st->d->resampler_buffer_output_frames =
 189                                     st->d->resampler_buffer_input_frames *
 190                                     st->d->oversample_factor;
 191   st->d->resampler_buffer_output = malloc
 192                                       (st->d->resampler_buffer_output_frames *
 193                                        st->channels *
 194                                        sizeof(float));
 195   CHECK_ERROR(!st->d->resampler_buffer_output, EBUR128_ERROR_NOMEM, free_input)
 196
 197   st->d->resampler = speex_resampler_init
 198                  ((spx_uint32_t) st->channels,
 199                   (spx_uint32_t) st->samplerate,
 200                   (spx_uint32_t) (st->samplerate * st->d->oversample_factor),
 201                   8, NULL);
 202   CHECK_ERROR(!st->d->resampler, EBUR128_ERROR_NOMEM, free_output)
 203
 204   return errcode;
 205
 206 free_output:
 207   free(st->d->resampler_buffer_output);
 208   st->d->resampler_buffer_output = NULL;
 209 free_input:
 210   free(st->d->resampler_buffer_input);
 211   st->d->resampler_buffer_input = NULL;
 212 exit:
 213   return errcode;
 214 }
 215
 216 static void ebur128_destroy_resampler(ebur128_state* st) {
 217   free(st->d->resampler_buffer_input);
 218   st->d->resampler_buffer_input = NULL;
 219   free(st->d->resampler_buffer_output);
 220   st->d->resampler_buffer_output = NULL;
 221   speex_resampler_destroy(st->d->resampler);
 222   st->d->resampler = NULL;
 223 }
 224 #endif
 225
 226 void ebur128_get_version(int* major, int* minor, int* patch) {
 227   *major = EBUR128_VERSION_MAJOR;
 228   *minor = EBUR128_VERSION_MINOR;
 229   *patch = EBUR128_VERSION_PATCH;
 230 }
 231
 232 ebur128_state* ebur128_init(unsigned int channels,
 233                             unsigned long samplerate,
 234                             int mode) {
 235   int errcode, result;
 236   ebur128_state* st;
 237   unsigned int i;
 238
 239   st = (ebur128_state*) malloc(sizeof(ebur128_state));
 240   CHECK_ERROR(!st, 0, exit)
 241   st->d = (struct ebur128_state_internal*)
 242           malloc(sizeof(struct ebur128_state_internal));
 243   CHECK_ERROR(!st->d, 0, free_state)
 244   st->channels = channels;
 245   errcode = ebur128_init_channel_map(st);
 246   CHECK_ERROR(errcode, 0, free_internal)
 247
 248   st->d->sample_peak = (double*) malloc(channels * sizeof(double));
 249   CHECK_ERROR(!st->d->sample_peak, 0, free_channel_map)
 250   st->d->true_peak = (double*) malloc(channels * sizeof(double));
 251   CHECK_ERROR(!st->d->true_peak, 0, free_sample_peak)
 252   for (i = 0; i < channels; ++i) {
 253     st->d->sample_peak[i] = 0.0;
 254     st->d->true_peak[i] = 0.0;
 255   }
 256
 257   st->d->use_histogram = mode & EBUR128_MODE_HISTOGRAM ? 1 : 0;
 258
 259   st->samplerate = samplerate;
 260   st->d->samples_in_100ms = (st->samplerate + 5) / 10;
 261   st->mode = mode;
 262   if ((mode & EBUR128_MODE_S) == EBUR128_MODE_S) {
 263     st->d->audio_data_frames = st->d->samples_in_100ms * 30;
 264   } else if ((mode & EBUR128_MODE_M) == EBUR128_MODE_M) {
 265     st->d->audio_data_frames = st->d->samples_in_100ms * 4;
 266   } else {
 267     return NULL;
 268   }
 269   st->d->audio_data = (double*) malloc(st->d->audio_data_frames *
 270                                        st->channels *
 271                                        sizeof(double));
 272   CHECK_ERROR(!st->d->audio_data, 0, free_true_peak)
 273   ebur128_init_filter(st);
 274
 275   if (st->d->use_histogram) {
 276     st->d->block_energy_histogram = malloc(1000 * sizeof(unsigned long));
 277     CHECK_ERROR(!st->d->block_energy_histogram, 0, free_audio_data)
 278     for (i = 0; i < 1000; ++i) {
 279       st->d->block_energy_histogram[i] = 0;
 280     }
 281   } else {
 282     st->d->block_energy_histogram = NULL;
 283   }
 284   if (st->d->use_histogram) {
 285     st->d->short_term_block_energy_histogram = malloc(1000 * sizeof(unsigned long));
 286     CHECK_ERROR(!st->d->short_term_block_energy_histogram, 0, free_block_energy_histogram)
 287     for (i = 0; i < 1000; ++i) {
 288       st->d->short_term_block_energy_histogram[i] = 0;
 289     }
 290   } else {
 291     st->d->short_term_block_energy_histogram = NULL;
 292   }
 293   SLIST_INIT(&st->d->block_list);
 294   SLIST_INIT(&st->d->short_term_block_list);
 295   st->d->short_term_frame_counter = 0;
 296
 297 #ifdef USE_SPEEX_RESAMPLER
 298   result = ebur128_init_resampler(st);
 299   CHECK_ERROR(result, 0, free_short_term_block_energy_histogram)
 300 #endif
 301
 302   /* the first block needs 400ms of audio data */
 303   st->d->needed_frames = st->d->samples_in_100ms * 4;
 304   /* start at the beginning of the buffer */
 305   st->d->audio_data_index = 0;
 306
 307   /* initialize static constants */
 308   relative_gate_factor = pow(10.0, relative_gate / 10.0);
 309   minus_twenty_decibels = pow(10.0, -20.0 / 10.0);
 310   histogram_energy_boundaries[0] = pow(10.0, (-70.0 + 0.691) / 10.0);
 311   if (st->d->use_histogram) {
 312     for (i = 0; i < 1000; ++i) {
 313       histogram_energies[i] = pow(10.0, ((double) i / 10.0 - 69.95 + 0.691) / 10.0);
 314     }
 315     for (i = 1; i < 1001; ++i) {
 316       histogram_energy_boundaries[i] = pow(10.0, ((double) i / 10.0 - 70.0 + 0.691) / 10.0);
 317     }
 318   }
 319
 320   return st;
 321
 322 free_short_term_block_energy_histogram:
 323   free(st->d->short_term_block_energy_histogram);
 324 free_block_energy_histogram:
 325   free(st->d->block_energy_histogram);
 326 free_audio_data:
 327   free(st->d->audio_data);
 328 free_true_peak:
 329   free(st->d->true_peak);
 330 free_sample_peak:
 331   free(st->d->sample_peak);
 332 free_channel_map:
 333   free(st->d->channel_map);
 334 free_internal:
 335   free(st->d);
 336 free_state:
 337   free(st);
 338 exit:
 339   return NULL;
 340 }
 341
 342 void ebur128_destroy(ebur128_state** st) {
 343   struct ebur128_dq_entry* entry;
 344   free((*st)->d->block_energy_histogram);
 345   free((*st)->d->short_term_block_energy_histogram);
 346   free((*st)->d->audio_data);
 347   free((*st)->d->channel_map);
 348   free((*st)->d->sample_peak);
 349   free((*st)->d->true_peak);
 350   while (!SLIST_EMPTY(&(*st)->d->block_list)) {
 351     entry = SLIST_FIRST(&(*st)->d->block_list);
 352     SLIST_REMOVE_HEAD(&(*st)->d->block_list, entries);
 353     free(entry);
 354   }
 355   while (!SLIST_EMPTY(&(*st)->d->short_term_block_list)) {
 356     entry = SLIST_FIRST(&(*st)->d->short_term_block_list);
 357     SLIST_REMOVE_HEAD(&(*st)->d->short_term_block_list, entries);
 358     free(entry);
 359   }
 360 #ifdef USE_SPEEX_RESAMPLER
 361   ebur128_destroy_resampler(*st);
 362 #endif
 363
 364   free((*st)->d);
 365   free(*st);
 366   *st = NULL;
 367 }
 368
 369 static int ebur128_use_speex_resampler(ebur128_state* st) {
 370 #ifdef USE_SPEEX_RESAMPLER
 371   return ((st->mode & EBUR128_MODE_TRUE_PEAK) == EBUR128_MODE_TRUE_PEAK);
 372 #else
 373   (void) st;
 374   return 0;
 375 #endif
 376 }
 377
 378 static void ebur128_check_true_peak(ebur128_state* st, size_t frames) {
 379 #ifdef USE_SPEEX_RESAMPLER
 380   size_t c, i;
 381   spx_uint32_t in_len = (spx_uint32_t) frames;
 382   spx_uint32_t out_len = (spx_uint32_t) st->d->resampler_buffer_output_frames;
 383   speex_resampler_process_interleaved_float(
 384                       st->d->resampler,
 385                       st->d->resampler_buffer_input,  &in_len,
 386                       st->d->resampler_buffer_output, &out_len);
 387   for (c = 0; c < st->channels; ++c) {
 388     for (i = 0; i < out_len; ++i) {
 389       if (st->d->resampler_buffer_output[i * st->channels + c] >
 390                                                          st->d->true_peak[c]) {
 391         st->d->true_peak[c] =
 392             st->d->resampler_buffer_output[i * st->channels + c];
 393       } else if (-st->d->resampler_buffer_output[i * st->channels + c] >
 394                                                          st->d->true_peak[c]) {
 395         st->d->true_peak[c] =
 396            -st->d->resampler_buffer_output[i * st->channels + c];
 397       }
 398     }
 399   }
 400 #else
 401   (void) st; (void) frames;
 402 #endif
 403 }
 404
 405 #ifdef __SSE2_MATH__
 406 #include <xmmintrin.h>
 407 #define TURN_ON_FTZ \
 408         unsigned int mxcsr = _mm_getcsr(); \
 409         _mm_setcsr(mxcsr | _MM_FLUSH_ZERO_ON);
 410 #define TURN_OFF_FTZ _mm_setcsr(mxcsr);
 411 #define FLUSH_MANUALLY
 412 #else
 413 #warning "manual FTZ is being used, please enable SSE2 (-msse2 -mfpmath=sse)"
 414 #define TURN_ON_FTZ
 415 #define TURN_OFF_FTZ
 416 #define FLUSH_MANUALLY \
 417     st->d->v[ci][4] = fabs(st->d->v[ci][4]) < DBL_MIN ? 0.0 : st->d->v[ci][4]; \
 418     st->d->v[ci][3] = fabs(st->d->v[ci][3]) < DBL_MIN ? 0.0 : st->d->v[ci][3]; \
 419     st->d->v[ci][2] = fabs(st->d->v[ci][2]) < DBL_MIN ? 0.0 : st->d->v[ci][2]; \
 420     st->d->v[ci][1] = fabs(st->d->v[ci][1]) < DBL_MIN ? 0.0 : st->d->v[ci][1];
 421 #endif
 422
 423 #define EBUR128_FILTER(type, min_scale, max_scale)                             \
 424 static void ebur128_filter_##type(ebur128_state* st, const type* src,          \
 425                                   size_t frames) {                             \
 426   static double scaling_factor = -((double) min_scale) > (double) max_scale ?  \
 427                                  -((double) min_scale) : (double) max_scale;   \
 428   double* audio_data = st->d->audio_data + st->d->audio_data_index;            \
 429   size_t i, c;                                                                 \
 430                                                                                \
 431   TURN_ON_FTZ                                                                  \
 432                                                                                \
 433   if ((st->mode & EBUR128_MODE_SAMPLE_PEAK) == EBUR128_MODE_SAMPLE_PEAK) {     \
 434     for (c = 0; c < st->channels; ++c) {                                       \
 435       double max = 0.0;                                                        \
 436       for (i = 0; i < frames; ++i) {                                           \
 437         if (src[i * st->channels + c] > max) {                                 \
 438           max =        src[i * st->channels + c];                              \
 439         } else if (-src[i * st->channels + c] > max) {                         \
 440           max = -1.0 * src[i * st->channels + c];                              \
 441         }                                                                      \
 442       }                                                                        \
 443       max /= scaling_factor;                                                   \
 444       if (max > st->d->sample_peak[c]) st->d->sample_peak[c] = max;            \
 445     }                                                                          \
 446   }                                                                            \
 447   if (ebur128_use_speex_resampler(st)) {                                       \
 448     for (c = 0; c < st->channels; ++c) {                                       \
 449       for (i = 0; i < frames; ++i) {                                           \
 450         st->d->resampler_buffer_input[i * st->channels + c] =                  \
 451                       (float) (src[i * st->channels + c] / scaling_factor);    \
 452       }                                                                        \
 453     }                                                                          \
 454     ebur128_check_true_peak(st, frames);                                       \
 455   }                                                                            \
 456   for (c = 0; c < st->channels; ++c) {                                         \
 457     int ci = st->d->channel_map[c] - 1;                                        \
 458     if (ci < 0) continue;                                                      \
 459     else if (ci > 4) ci = 0; /* dual mono */                                   \
 460     for (i = 0; i < frames; ++i) {                                             \
 461       st->d->v[ci][0] = (double) (src[i * st->channels + c] / scaling_factor)  \
 462                    - st->d->a[1] * st->d->v[ci][1]                             \
 463                    - st->d->a[2] * st->d->v[ci][2]                             \
 464                    - st->d->a[3] * st->d->v[ci][3]                             \
 465                    - st->d->a[4] * st->d->v[ci][4];                            \
 466       audio_data[i * st->channels + c] =                                       \
 467                      st->d->b[0] * st->d->v[ci][0]                             \
 468                    + st->d->b[1] * st->d->v[ci][1]                             \
 469                    + st->d->b[2] * st->d->v[ci][2]                             \
 470                    + st->d->b[3] * st->d->v[ci][3]                             \
 471                    + st->d->b[4] * st->d->v[ci][4];                            \
 472       st->d->v[ci][4] = st->d->v[ci][3];                                       \
 473       st->d->v[ci][3] = st->d->v[ci][2];                                       \
 474       st->d->v[ci][2] = st->d->v[ci][1];                                       \
 475       st->d->v[ci][1] = st->d->v[ci][0];                                       \
 476     }                                                                          \
 477     FLUSH_MANUALLY                                                             \
 478   }                                                                            \
 479   TURN_OFF_FTZ                                                                 \
 480 }
 481 EBUR128_FILTER(short, SHRT_MIN, SHRT_MAX)
 482 EBUR128_FILTER(int, INT_MIN, INT_MAX)
 483 EBUR128_FILTER(float, -1.0f, 1.0f)
 484 EBUR128_FILTER(double, -1.0, 1.0)
 485
 486 static double ebur128_energy_to_loudness(double energy) {
 487   return 10 * (log(energy) / log(10.0)) - 0.691;
 488 }
 489
 490 static size_t find_histogram_index(double energy) {
 491   size_t index_min = 0;
 492   size_t index_max = 1000;
 493   size_t index_mid;
 494
 495   do {
 496     index_mid = (index_min + index_max) / 2;
 497     if (energy >= histogram_energy_boundaries[index_mid]) {
 498       index_min = index_mid;
 499     } else {
 500       index_max = index_mid;
 501     }
 502   } while (index_max - index_min != 1);
 503
 504   return index_min;
 505 }
 506
 507 static int ebur128_calc_gating_block(ebur128_state* st, size_t frames_per_block,
 508                                      double* optional_output) {
 509   size_t i, c;
 510   double sum = 0.0;
 511   double channel_sum;
 512   for (c = 0; c < st->channels; ++c) {
 513     if (st->d->channel_map[c] == EBUR128_UNUSED) continue;
 514     channel_sum = 0.0;
 515     if (st->d->audio_data_index < frames_per_block * st->channels) {
 516       for (i = 0; i < st->d->audio_data_index / st->channels; ++i) {
 517         channel_sum += st->d->audio_data[i * st->channels + c] *
 518                        st->d->audio_data[i * st->channels + c];
 519       }
 520       for (i = st->d->audio_data_frames -
 521               (frames_per_block -
 522                st->d->audio_data_index / st->channels);
 523            i < st->d->audio_data_frames; ++i) {
 524         channel_sum += st->d->audio_data[i * st->channels + c] *
 525                        st->d->audio_data[i * st->channels + c];
 526       }
 527     } else {
 528       for (i = st->d->audio_data_index / st->channels - frames_per_block;
 529            i < st->d->audio_data_index / st->channels;
 530            ++i) {
 531         channel_sum += st->d->audio_data[i * st->channels + c] *
 532                        st->d->audio_data[i * st->channels + c];
 533       }
 534     }
 535     if (st->d->channel_map[c] == EBUR128_LEFT_SURROUND ||
 536         st->d->channel_map[c] == EBUR128_RIGHT_SURROUND) {
 537       channel_sum *= 1.41;
 538     } else if (st->d->channel_map[c] == EBUR128_DUAL_MONO) {
 539       channel_sum *= 2.0;
 540     }
 541     sum += channel_sum;
 542   }
 543   sum /= (double) frames_per_block;
 544   if (optional_output) {
 545     *optional_output = sum;
 546     return EBUR128_SUCCESS;
 547   } else if (sum >= histogram_energy_boundaries[0]) {
 548     if (st->d->use_histogram) {
 549       ++st->d->block_energy_histogram[find_histogram_index(sum)];
 550     } else {
 551       struct ebur128_dq_entry* block;
 552       block = (struct ebur128_dq_entry*) malloc(sizeof(struct ebur128_dq_entry));
 553       if (!block) return EBUR128_ERROR_NOMEM;
 554       block->z = sum;
 555       SLIST_INSERT_HEAD(&st->d->block_list, block, entries);
 556     }
 557     return EBUR128_SUCCESS;
 558   } else {
 559     return EBUR128_SUCCESS;
 560   }
 561 }
 562
 563 int ebur128_set_channel(ebur128_state* st,
 564                         unsigned int channel_number,
 565                         int value) {
 566   if (channel_number >= st->channels) {
 567     return 1;
 568   }
 569   if (value == EBUR128_DUAL_MONO &&
 570       (st->channels != 1 || channel_number != 0)) {
 571     fprintf(stderr, "EBUR128_DUAL_MONO only works with mono files!\n");
 572     return 1;
 573   }
 574   st->d->channel_map[channel_number] = value;
 575   return 0;
 576 }
 577
 578 int ebur128_change_parameters(ebur128_state* st,
 579                               unsigned int channels,
 580                               unsigned long samplerate) {
 581   int errcode;
 582   if (channels == st->channels &&
 583       samplerate == st->samplerate) {
 584     return 2;
 585   }
 586   free(st->d->audio_data);
 587   st->d->audio_data = NULL;
 588
 589   if (channels != st->channels) {
 590     unsigned int i;
 591
 592     free(st->d->channel_map); st->d->channel_map = NULL;
 593     free(st->d->sample_peak); st->d->sample_peak = NULL;
 594     free(st->d->true_peak);   st->d->true_peak = NULL;
 595     st->channels = channels;
 596
 597 #ifdef USE_SPEEX_RESAMPLER
 598     ebur128_destroy_resampler(st);
 599     ebur128_init_resampler(st);
 600 #endif
 601
 602     errcode = ebur128_init_channel_map(st);
 603     CHECK_ERROR(errcode, EBUR128_ERROR_NOMEM, exit)
 604
 605     st->d->sample_peak = (double*) malloc(channels * sizeof(double));
 606     CHECK_ERROR(!st->d->sample_peak, EBUR128_ERROR_NOMEM, exit)
 607     st->d->true_peak = (double*) malloc(channels * sizeof(double));
 608     CHECK_ERROR(!st->d->true_peak, EBUR128_ERROR_NOMEM, exit)
 609     for (i = 0; i < channels; ++i) {
 610       st->d->sample_peak[i] = 0.0;
 611       st->d->true_peak[i] = 0.0;
 612     }
 613   }
 614   if (samplerate != st->samplerate) {
 615     st->samplerate = samplerate;
 616     ebur128_init_filter(st);
 617   }
 618   if ((st->mode & EBUR128_MODE_S) == EBUR128_MODE_S) {
 619     st->d->audio_data_frames = st->d->samples_in_100ms * 30;
 620   } else if ((st->mode & EBUR128_MODE_M) == EBUR128_MODE_M) {
 621     st->d->audio_data_frames = st->d->samples_in_100ms * 4;
 622   } else {
 623     return 1;
 624   }
 625   st->d->audio_data = (double*) malloc(st->d->audio_data_frames *
 626                                        st->channels *
 627                                        sizeof(double));
 628   CHECK_ERROR(!st->d->audio_data, EBUR128_ERROR_NOMEM, exit)
 629
 630   /* the first block needs 400ms of audio data */
 631   st->d->needed_frames = st->d->samples_in_100ms * 4;
 632   /* start at the beginning of the buffer */
 633   st->d->audio_data_index = 0;
 634   /* reset short term frame counter */
 635   st->d->short_term_frame_counter = 0;
 636
 637   return 0;
 638
 639 exit:
 640   return 1;
 641 }
 642
 643
 644 static int ebur128_energy_shortterm(ebur128_state* st, double* out);
 645 #define EBUR128_ADD_FRAMES(type)                                               \
 646 int ebur128_add_frames_##type(ebur128_state* st,                               \
 647                               const type* src, size_t frames) {                \
 648   size_t src_index = 0;                                                        \
 649   while (frames > 0) {                                                         \
 650     if (frames >= st->d->needed_frames) {                                      \
 651       ebur128_filter_##type(st, src + src_index, st->d->needed_frames);        \
 652       src_index += st->d->needed_frames * st->channels;                        \
 653       frames -= st->d->needed_frames;                                          \
 654       st->d->audio_data_index += st->d->needed_frames * st->channels;          \
 655       /* calculate the new gating block */                                     \
 656       if ((st->mode & EBUR128_MODE_I) == EBUR128_MODE_I) {                     \
 657         if (ebur128_calc_gating_block(st, st->d->samples_in_100ms * 4, NULL)) {\
 658           return EBUR128_ERROR_NOMEM;                                          \
 659         }                                                                      \
 660       }                                                                        \
 661       if ((st->mode & EBUR128_MODE_LRA) == EBUR128_MODE_LRA) {                 \
 662         st->d->short_term_frame_counter += st->d->needed_frames;               \
 663         if (st->d->short_term_frame_counter == st->d->samples_in_100ms * 30) { \
 664           struct ebur128_dq_entry* block;                                      \
 665           double st_energy;                                                    \
 666           ebur128_energy_shortterm(st, &st_energy);                            \
 667           if (st_energy >= histogram_energy_boundaries[0]) {                   \
 668             if (st->d->use_histogram) {                                        \
 669               ++st->d->short_term_block_energy_histogram[                      \
 670                                               find_histogram_index(st_energy)];\
 671             } else {                                                           \
 672               block = (struct ebur128_dq_entry*)                               \
 673                       malloc(sizeof(struct ebur128_dq_entry));                 \
 674               if (!block) return EBUR128_ERROR_NOMEM;                          \
 675               block->z = st_energy;                                            \
 676               SLIST_INSERT_HEAD(&st->d->short_term_block_list, block, entries);\
 677             }                                                                  \
 678           }                                                                    \
 679           st->d->short_term_frame_counter = st->d->samples_in_100ms * 20;      \
 680         }                                                                      \
 681       }                                                                        \
 682       /* 100ms are needed for all blocks besides the first one */              \
 683       st->d->needed_frames = st->d->samples_in_100ms;                          \
 684       /* reset audio_data_index when buffer full */                            \
 685       if (st->d->audio_data_index == st->d->audio_data_frames * st->channels) {\
 686         st->d->audio_data_index = 0;                                           \
 687       }                                                                        \
 688     } else {                                                                   \
 689       ebur128_filter_##type(st, src + src_index, frames);                      \
 690       st->d->audio_data_index += frames * st->channels;                        \
 691       if ((st->mode & EBUR128_MODE_LRA) == EBUR128_MODE_LRA) {                 \
 692         st->d->short_term_frame_counter += frames;                             \
 693       }                                                                        \
 694       st->d->needed_frames -= frames;                                          \
 695       frames = 0;                                                              \
 696     }                                                                          \
 697   }                                                                            \
 698   return EBUR128_SUCCESS;                                                      \
 699 }
 700 EBUR128_ADD_FRAMES(short)
 701 EBUR128_ADD_FRAMES(int)
 702 EBUR128_ADD_FRAMES(float)
 703 EBUR128_ADD_FRAMES(double)
 704
 705 static int ebur128_gated_loudness(ebur128_state** sts, size_t size,
 706                                   double* out) {
 707   struct ebur128_dq_entry* it;
 708   double relative_threshold = 0.0;
 709   double gated_loudness = 0.0;
 710   size_t above_thresh_counter = 0;
 711   size_t i, j, start_index;
 712
 713   for (i = 0; i < size; i++) {
 714     if (sts[i] && (sts[i]->mode & EBUR128_MODE_I) != EBUR128_MODE_I) {
 715       return EBUR128_ERROR_INVALID_MODE;
 716     }
 717   }
 718
 719   for (i = 0; i < size; i++) {
 720     if (!sts[i]) continue;
 721     if (sts[i]->d->use_histogram) {
 722       for (j = 0; j < 1000; ++j) {
 723         relative_threshold += sts[i]->d->block_energy_histogram[j] *
 724                               histogram_energies[j];
 725         above_thresh_counter += sts[i]->d->block_energy_histogram[j];
 726       }
 727     } else {
 728       SLIST_FOREACH(it, &sts[i]->d->block_list, entries) {
 729         ++above_thresh_counter;
 730         relative_threshold += it->z;
 731       }
 732     }
 733   }
 734   if (!above_thresh_counter) {
 735     *out = -HUGE_VAL;
 736     return EBUR128_SUCCESS;
 737   }
 738   relative_threshold /= (double) above_thresh_counter;
 739   relative_threshold *= relative_gate_factor;
 740   above_thresh_counter = 0;
 741   if (relative_threshold < histogram_energy_boundaries[0]) {
 742     start_index = 0;
 743   } else {
 744     start_index = find_histogram_index(relative_threshold);
 745     if (relative_threshold > histogram_energies[start_index]) {
 746       ++start_index;
 747     }
 748   }
 749   for (i = 0; i < size; i++) {
 750     if (!sts[i]) continue;
 751     if (sts[i]->d->use_histogram) {
 752       for (j = start_index; j < 1000; ++j) {
 753         gated_loudness += sts[i]->d->block_energy_histogram[j] *
 754                           histogram_energies[j];
 755         above_thresh_counter += sts[i]->d->block_energy_histogram[j];
 756       }
 757     } else {
 758       SLIST_FOREACH(it, &sts[i]->d->block_list, entries) {
 759         if (it->z >= relative_threshold) {
 760           ++above_thresh_counter;
 761           gated_loudness += it->z;
 762         }
 763       }
 764     }
 765   }
 766   if (!above_thresh_counter) {
 767     *out = -HUGE_VAL;
 768     return EBUR128_SUCCESS;
 769   }
 770   gated_loudness /= (double) above_thresh_counter;
 771   *out = ebur128_energy_to_loudness(gated_loudness);
 772   return EBUR128_SUCCESS;
 773 }
 774
 775 int ebur128_loudness_global(ebur128_state* st, double* out) {
 776   return ebur128_gated_loudness(&st, 1, out);
 777 }
 778
 779 int ebur128_loudness_global_multiple(ebur128_state** sts, size_t size,
 780                                      double* out) {
 781   return ebur128_gated_loudness(sts, size, out);
 782 }
 783
 784 static int ebur128_energy_in_interval(ebur128_state* st,
 785                                       size_t interval_frames,
 786                                       double* out) {
 787   if (interval_frames > st->d->audio_data_frames) {
 788     return EBUR128_ERROR_INVALID_MODE;
 789   }
 790   ebur128_calc_gating_block(st, interval_frames, out);
 791   return EBUR128_SUCCESS;
 792 }
 793
 794 static int ebur128_energy_shortterm(ebur128_state* st, double* out) {
 795   return ebur128_energy_in_interval(st, st->d->samples_in_100ms * 30, out);
 796 }
 797
 798 int ebur128_loudness_momentary(ebur128_state* st, double* out) {
 799   double energy;
 800   int error = ebur128_energy_in_interval(st, st->d->samples_in_100ms * 4,
 801                                          &energy);
 802   if (error) {
 803     return error;
 804   } else if (energy <= 0.0) {
 805     *out = -HUGE_VAL;
 806     return EBUR128_SUCCESS;
 807   }
 808   *out = ebur128_energy_to_loudness(energy);
 809   return EBUR128_SUCCESS;
 810 }
 811
 812 int ebur128_loudness_shortterm(ebur128_state* st, double* out) {
 813   double energy;
 814   int error = ebur128_energy_shortterm(st, &energy);
 815   if (error) {
 816     return error;
 817   } else if (energy <= 0.0) {
 818     *out = -HUGE_VAL;
 819     return EBUR128_SUCCESS;
 820   }
 821   *out = ebur128_energy_to_loudness(energy);
 822   return EBUR128_SUCCESS;
 823 }
 824
 825 static int ebur128_double_cmp(const void *p1, const void *p2) {
 826   const double* d1 = (const double*) p1;
 827   const double* d2 = (const double*) p2;
 828   return (*d1 > *d2) - (*d1 < *d2);
 829 }
 830
 831 /* EBU - TECH 3342 */
 832 int ebur128_loudness_range_multiple(ebur128_state** sts, size_t size,
 833                                     double* out) {
 834   size_t i, j;
 835   struct ebur128_dq_entry* it;
 836   double* stl_vector;
 837   size_t stl_size;
 838   double* stl_relgated;
 839   size_t stl_relgated_size;
 840   double stl_power, stl_integrated;
 841   /* High and low percentile energy */
 842   double h_en, l_en;
 843   int use_histogram = 0;
 844
 845   for (i = 0; i < size; ++i) {
 846     if (sts[i]) {
 847       if ((sts[i]->mode & EBUR128_MODE_LRA) != EBUR128_MODE_LRA) {
 848         return EBUR128_ERROR_INVALID_MODE;
 849       }
 850       if (i == 0 && sts[i]->mode & EBUR128_MODE_HISTOGRAM) {
 851         use_histogram = 1;
 852       } else if (use_histogram != !!(sts[i]->mode & EBUR128_MODE_HISTOGRAM)) {
 853         return EBUR128_ERROR_INVALID_MODE;
 854       }
 855     }
 856   }
 857
 858   if (use_histogram) {
 859     unsigned long hist[1000] = { 0 };
 860     size_t percentile_low, percentile_high;
 861     size_t index;
 862
 863     stl_size = 0;
 864     stl_power = 0.0;
 865     for (i = 0; i < size; ++i) {
 866       if (!sts[i]) continue;
 867       for (j = 0; j < 1000; ++j) {
 868         hist[j]   += sts[i]->d->short_term_block_energy_histogram[j];
 869         stl_size  += sts[i]->d->short_term_block_energy_histogram[j];
 870         stl_power += sts[i]->d->short_term_block_energy_histogram[j]
 871                      * histogram_energies[j];
 872       }
 873     }
 874     if (!stl_size) {
 875       *out = 0.0;
 876       return EBUR128_SUCCESS;
 877     }
 878
 879     stl_power /= stl_size;
 880     stl_integrated = minus_twenty_decibels * stl_power;
 881
 882     if (stl_integrated < histogram_energy_boundaries[0]) {
 883       index = 0;
 884     } else {
 885       index = find_histogram_index(stl_integrated);
 886       if (stl_integrated > histogram_energies[index]) {
 887         ++index;
 888       }
 889     }
 890     stl_size = 0;
 891     for (j = index; j < 1000; ++j) {
 892       stl_size += hist[j];
 893     }
 894     if (!stl_size) {
 895       *out = 0.0;
 896       return EBUR128_SUCCESS;
 897     }
 898
 899     percentile_low  = (size_t) ((stl_size - 1) * 0.1 + 0.5);
 900     percentile_high = (size_t) ((stl_size - 1) * 0.95 + 0.5);
 901
 902     stl_size = 0;
 903     j = index;
 904     while (stl_size <= percentile_low) {
 905       stl_size += hist[j++];
 906     }
 907     l_en = histogram_energies[j - 1];
 908     while (stl_size <= percentile_high) {
 909       stl_size += hist[j++];
 910     }
 911     h_en = histogram_energies[j - 1];
 912     *out = ebur128_energy_to_loudness(h_en) - ebur128_energy_to_loudness(l_en);
 913     return EBUR128_SUCCESS;
 914
 915   } else {
 916     stl_size = 0;
 917     for (i = 0; i < size; ++i) {
 918       if (!sts[i]) continue;
 919       SLIST_FOREACH(it, &sts[i]->d->short_term_block_list, entries) {
 920         ++stl_size;
 921       }
 922     }
 923     if (!stl_size) {
 924       *out = 0.0;
 925       return EBUR128_SUCCESS;
 926     }
 927     stl_vector = (double*) malloc(stl_size * sizeof(double));
 928     if (!stl_vector)
 929       return EBUR128_ERROR_NOMEM;
 930
 931     for (j = 0, i = 0; i < size; ++i) {
 932       if (!sts[i]) continue;
 933       SLIST_FOREACH(it, &sts[i]->d->short_term_block_list, entries) {
 934         stl_vector[j] = it->z;
 935         ++j;
 936       }
 937     }
 938     qsort(stl_vector, stl_size, sizeof(double), ebur128_double_cmp);
 939     stl_power = 0.0;
 940     for (i = 0; i < stl_size; ++i) {
 941       stl_power += stl_vector[i];
 942     }
 943     stl_power /= (double) stl_size;
 944     stl_integrated = minus_twenty_decibels * stl_power;
 945
 946     stl_relgated = stl_vector;
 947     stl_relgated_size = stl_size;
 948     while (stl_relgated_size > 0 && *stl_relgated < stl_integrated) {
 949       ++stl_relgated;
 950       --stl_relgated_size;
 951     }
 952
 953     if (stl_relgated_size) {
 954       h_en = stl_relgated[(size_t) ((stl_relgated_size - 1) * 0.95 + 0.5)];
 955       l_en = stl_relgated[(size_t) ((stl_relgated_size - 1) * 0.1 + 0.5)];
 956       free(stl_vector);
 957       *out = ebur128_energy_to_loudness(h_en) - ebur128_energy_to_loudness(l_en);
 958       return EBUR128_SUCCESS;
 959     } else {
 960       free(stl_vector);
 961       *out = 0.0;
 962       return EBUR128_SUCCESS;
 963     }
 964   }
 965 }
 966
 967 int ebur128_loudness_range(ebur128_state* st, double* out) {
 968   return ebur128_loudness_range_multiple(&st, 1, out);
 969 }
 970
 971 int ebur128_sample_peak(ebur128_state* st,
 972                         unsigned int channel_number,
 973                         double* out) {
 974   if ((st->mode & EBUR128_MODE_SAMPLE_PEAK) != EBUR128_MODE_SAMPLE_PEAK) {
 975     return EBUR128_ERROR_INVALID_MODE;
 976   } else if (channel_number >= st->channels) {
 977     return EBUR128_ERROR_INVALID_CHANNEL_INDEX;
 978   }
 979   *out = st->d->sample_peak[channel_number];
 980   return EBUR128_SUCCESS;
 981 }
 982
 983 #ifdef USE_SPEEX_RESAMPLER
 984 int ebur128_true_peak(ebur128_state* st,
 985                       unsigned int channel_number,
 986                       double* out) {
 987   if ((st->mode & EBUR128_MODE_TRUE_PEAK) != EBUR128_MODE_TRUE_PEAK) {
 988     return EBUR128_ERROR_INVALID_MODE;
 989   } else if (channel_number >= st->channels) {
 990     return EBUR128_ERROR_INVALID_CHANNEL_INDEX;
 991   }
 992   *out = st->d->true_peak[channel_number] > st->d->sample_peak[channel_number]
 993        ? st->d->true_peak[channel_number]
 994        : st->d->sample_peak[channel_number];
 995   return EBUR128_SUCCESS;
 996 }
 997 #endif