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 <sys/queue.h>
13 
14 #define CHECK_ERROR(condition, errorcode, goto_point)                          \
15   if ((condition)) {                                                           \
16     errcode = (errorcode);                                                     \
17     goto goto_point;                                                           \
18   }
19 #define EBUR128_MAX(a, b) (((a) > (b)) ? (a) : (b))
20 
safe_size_mul(size_t nmemb,size_t size,size_t * result)21 static int safe_size_mul(size_t nmemb, size_t size, size_t* result) {
22   /* Adapted from OpenBSD reallocarray. */
23 #define MUL_NO_OVERFLOW (((size_t) 1) << (sizeof(size_t) * 4))
24   if ((nmemb >= MUL_NO_OVERFLOW || size >= MUL_NO_OVERFLOW) && /**/
25       nmemb > 0 && ((size_t) (-1)) / nmemb < size) {
26     return 1;
27   }
28 #undef MUL_NO_OVERFLOW
29   *result = nmemb * size;
30   return 0;
31 }
32 
33 STAILQ_HEAD(ebur128_double_queue, ebur128_dq_entry);
34 struct ebur128_dq_entry {
35   double z;
36   STAILQ_ENTRY(ebur128_dq_entry) entries;
37 };
38 
39 #define ALMOST_ZERO 0.000001
40 #define FILTER_STATE_SIZE 5
41 
42 typedef struct {
43   unsigned int count;  /* Number of coefficients in this subfilter */
44   unsigned int* index; /* Delay index of corresponding filter coeff */
45   double* coeff;       /* List of subfilter coefficients */
46 } interp_filter;
47 
48 typedef struct {         /* Data structure for polyphase FIR interpolator */
49   unsigned int factor;   /* Interpolation factor of the interpolator */
50   unsigned int taps;     /* Taps (prefer odd to increase zero coeffs) */
51   unsigned int channels; /* Number of channels */
52   unsigned int delay;    /* Size of delay buffer */
53   interp_filter* filter; /* List of subfilters (one for each factor) */
54   float** z;             /* List of delay buffers (one for each channel) */
55   unsigned int zi;       /* Current delay buffer index */
56 } interpolator;
57 
58 /** BS.1770 filter state. */
59 typedef double filter_state[FILTER_STATE_SIZE];
60 
61 struct ebur128_state_internal {
62   /** Filtered audio data (used as ring buffer). */
63   double* audio_data;
64   /** Size of audio_data array. */
65   size_t audio_data_frames;
66   /** Current index for audio_data. */
67   size_t audio_data_index;
68   /** How many frames are needed for a gating block. Will correspond to 400ms
69    *  of audio at initialization, and 100ms after the first block (75% overlap
70    *  as specified in the 2011 revision of BS1770). */
71   unsigned long needed_frames;
72   /** The channel map. Has as many elements as there are channels. */
73   int* channel_map;
74   /** How many samples fit in 100ms (rounded). */
75   unsigned long samples_in_100ms;
76   /** BS.1770 filter coefficients (nominator). */
77   double b[5];
78   /** BS.1770 filter coefficients (denominator). */
79   double a[5];
80   /** one filter_state per channel. */
81   filter_state* v;
82   /** Linked list of block energies. */
83   struct ebur128_double_queue block_list;
84   unsigned long block_list_max;
85   unsigned long block_list_size;
86   /** Linked list of 3s-block energies, used to calculate LRA. */
87   struct ebur128_double_queue short_term_block_list;
88   unsigned long st_block_list_max;
89   unsigned long st_block_list_size;
90   int use_histogram;
91   unsigned long* block_energy_histogram;
92   unsigned long* short_term_block_energy_histogram;
93   /** Keeps track of when a new short term block is needed. */
94   size_t short_term_frame_counter;
95   /** Maximum sample peak, one per channel */
96   double* sample_peak;
97   double* prev_sample_peak;
98   /** Maximum true peak, one per channel */
99   double* true_peak;
100   double* prev_true_peak;
101   interpolator* interp;
102   float* resampler_buffer_input;
103   size_t resampler_buffer_input_frames;
104   float* resampler_buffer_output;
105   size_t resampler_buffer_output_frames;
106   /** The maximum window duration in ms. */
107   unsigned long window;
108   unsigned long history;
109 };
110 
111 static double relative_gate = -10.0;
112 
113 /* Those will be calculated when initializing the library */
114 static double relative_gate_factor;
115 static double minus_twenty_decibels;
116 static double histogram_energies[1000];
117 static double histogram_energy_boundaries[1001];
118 
119 static interpolator*
interp_create(unsigned int taps,unsigned int factor,unsigned int channels)120 interp_create(unsigned int taps, unsigned int factor, unsigned int channels) {
121   int errcode; /* unused */
122   interpolator* interp;
123   unsigned int j;
124 
125   interp = (interpolator*) calloc(1, sizeof(interpolator));
126   CHECK_ERROR(!interp, 0, exit);
127 
128   interp->taps = taps;
129   interp->factor = factor;
130   interp->channels = channels;
131   interp->delay = (interp->taps + interp->factor - 1) / interp->factor;
132 
133   /* Initialize the filter memory
134    * One subfilter per interpolation factor. */
135   interp->filter =
136       (interp_filter*) calloc(interp->factor, sizeof(*interp->filter));
137   CHECK_ERROR(!interp->filter, 0, free_interp);
138 
139   for (j = 0; j < interp->factor; j++) {
140     interp->filter[j].index =
141         (unsigned int*) calloc(interp->delay, sizeof(unsigned int));
142     interp->filter[j].coeff = (double*) calloc(interp->delay, sizeof(double));
143     CHECK_ERROR(!interp->filter[j].index || !interp->filter[j].coeff, 0,
144                 free_filter_index_coeff);
145   }
146 
147   /* One delay buffer per channel. */
148   interp->z = (float**) calloc(interp->channels, sizeof(float*));
149   CHECK_ERROR(!interp->z, 0, free_filter_index_coeff);
150   for (j = 0; j < interp->channels; j++) {
151     interp->z[j] = (float*) calloc(interp->delay, sizeof(float));
152     CHECK_ERROR(!interp->z[j], 0, free_filter_z);
153   }
154 
155   /* Calculate the filter coefficients */
156   for (j = 0; j < interp->taps; j++) {
157     /* Calculate sinc */
158     double m = (double) j - (double) (interp->taps - 1) / 2.0;
159     double c = 1.0;
160     if (fabs(m) > ALMOST_ZERO) {
161       c = sin(m * M_PI / interp->factor) / (m * M_PI / interp->factor);
162     }
163     /* Apply Hanning window */
164     c *= 0.5 * (1 - cos(2 * M_PI * j / (interp->taps - 1)));
165 
166     if (fabs(c) > ALMOST_ZERO) { /* Ignore any zero coeffs. */
167       /* Put the coefficient into the correct subfilter */
168       unsigned int f = j % interp->factor;
169       unsigned int t = interp->filter[f].count++;
170       interp->filter[f].coeff[t] = c;
171       interp->filter[f].index[t] = j / interp->factor;
172     }
173   }
174   return interp;
175 
176 free_filter_z:
177   for (j = 0; j < interp->channels; j++) {
178     free(interp->z[j]);
179   }
180   free(interp->z);
181 free_filter_index_coeff:
182   for (j = 0; j < interp->factor; j++) {
183     free(interp->filter[j].index);
184     free(interp->filter[j].coeff);
185   }
186   free(interp->filter);
187 free_interp:
188   free(interp);
189 exit:
190   return NULL;
191 }
192 
interp_destroy(interpolator * interp)193 static void interp_destroy(interpolator* interp) {
194   unsigned int j = 0;
195   if (!interp) {
196     return;
197   }
198   for (j = 0; j < interp->factor; j++) {
199     free(interp->filter[j].index);
200     free(interp->filter[j].coeff);
201   }
202   free(interp->filter);
203   for (j = 0; j < interp->channels; j++) {
204     free(interp->z[j]);
205   }
206   free(interp->z);
207   free(interp);
208 }
209 
210 static size_t
interp_process(interpolator * interp,size_t frames,float * in,float * out)211 interp_process(interpolator* interp, size_t frames, float* in, float* out) {
212   size_t frame = 0;
213   unsigned int chan = 0;
214   unsigned int f = 0;
215   unsigned int t = 0;
216   unsigned int out_stride = interp->channels * interp->factor;
217   float* outp = 0;
218   double acc = 0;
219   double c = 0;
220 
221   for (frame = 0; frame < frames; frame++) {
222     for (chan = 0; chan < interp->channels; chan++) {
223       /* Add sample to delay buffer */
224       interp->z[chan][interp->zi] = *in++;
225       /* Apply coefficients */
226       outp = out + chan;
227       for (f = 0; f < interp->factor; f++) {
228         acc = 0.0;
229         for (t = 0; t < interp->filter[f].count; t++) {
230           int i = (int) interp->zi - (int) interp->filter[f].index[t];
231           if (i < 0) {
232             i += (int) interp->delay;
233           }
234           c = interp->filter[f].coeff[t];
235           acc += (double) interp->z[chan][i] * c;
236         }
237         *outp = (float) acc;
238         outp += interp->channels;
239       }
240     }
241     out += out_stride;
242     interp->zi++;
243     if (interp->zi == interp->delay) {
244       interp->zi = 0;
245     }
246   }
247 
248   return frames * interp->factor;
249 }
250 
ebur128_init_filter(ebur128_state * st)251 static int ebur128_init_filter(ebur128_state* st) {
252   int errcode = EBUR128_SUCCESS;
253   int i, j;
254 
255   double f0 = 1681.974450955533;
256   double G = 3.999843853973347;
257   double Q = 0.7071752369554196;
258 
259   double K = tan(M_PI * f0 / (double) st->samplerate);
260   double Vh = pow(10.0, G / 20.0);
261   double Vb = pow(Vh, 0.4996667741545416);
262 
263   double pb[3] = { 0.0, 0.0, 0.0 };
264   double pa[3] = { 1.0, 0.0, 0.0 };
265   double rb[3] = { 1.0, -2.0, 1.0 };
266   double ra[3] = { 1.0, 0.0, 0.0 };
267 
268   double a0 = 1.0 + K / Q + K * K;
269   pb[0] = (Vh + Vb * K / Q + K * K) / a0;
270   pb[1] = 2.0 * (K * K - Vh) / a0;
271   pb[2] = (Vh - Vb * K / Q + K * K) / a0;
272   pa[1] = 2.0 * (K * K - 1.0) / a0;
273   pa[2] = (1.0 - K / Q + K * K) / a0;
274 
275   /* fprintf(stderr, "%.14f %.14f %.14f %.14f %.14f\n",
276                      b1[0], b1[1], b1[2], a1[1], a1[2]); */
277 
278   f0 = 38.13547087602444;
279   Q = 0.5003270373238773;
280   K = tan(M_PI * f0 / (double) st->samplerate);
281 
282   ra[1] = 2.0 * (K * K - 1.0) / (1.0 + K / Q + K * K);
283   ra[2] = (1.0 - K / Q + K * K) / (1.0 + K / Q + K * K);
284 
285   /* fprintf(stderr, "%.14f %.14f\n", a2[1], a2[2]); */
286 
287   st->d->b[0] = pb[0] * rb[0];
288   st->d->b[1] = pb[0] * rb[1] + pb[1] * rb[0];
289   st->d->b[2] = pb[0] * rb[2] + pb[1] * rb[1] + pb[2] * rb[0];
290   st->d->b[3] = pb[1] * rb[2] + pb[2] * rb[1];
291   st->d->b[4] = pb[2] * rb[2];
292 
293   st->d->a[0] = pa[0] * ra[0];
294   st->d->a[1] = pa[0] * ra[1] + pa[1] * ra[0];
295   st->d->a[2] = pa[0] * ra[2] + pa[1] * ra[1] + pa[2] * ra[0];
296   st->d->a[3] = pa[1] * ra[2] + pa[2] * ra[1];
297   st->d->a[4] = pa[2] * ra[2];
298 
299   st->d->v = (filter_state*) malloc(st->channels * sizeof(filter_state));
300   CHECK_ERROR(!st->d->v, EBUR128_ERROR_NOMEM, exit);
301   for (i = 0; i < (int) st->channels; ++i) {
302     for (j = 0; j < FILTER_STATE_SIZE; ++j) {
303       st->d->v[i][j] = 0.0;
304     }
305   }
306 
307 exit:
308   return errcode;
309 }
310 
ebur128_init_channel_map(ebur128_state * st)311 static int ebur128_init_channel_map(ebur128_state* st) {
312   size_t i;
313   st->d->channel_map = (int*) malloc(st->channels * sizeof(int));
314   if (!st->d->channel_map) {
315     return EBUR128_ERROR_NOMEM;
316   }
317   if (st->channels == 4) {
318     st->d->channel_map[0] = EBUR128_LEFT;
319     st->d->channel_map[1] = EBUR128_RIGHT;
320     st->d->channel_map[2] = EBUR128_LEFT_SURROUND;
321     st->d->channel_map[3] = EBUR128_RIGHT_SURROUND;
322   } else if (st->channels == 5) {
323     st->d->channel_map[0] = EBUR128_LEFT;
324     st->d->channel_map[1] = EBUR128_RIGHT;
325     st->d->channel_map[2] = EBUR128_CENTER;
326     st->d->channel_map[3] = EBUR128_LEFT_SURROUND;
327     st->d->channel_map[4] = EBUR128_RIGHT_SURROUND;
328   } else {
329     for (i = 0; i < st->channels; ++i) {
330       switch (i) {
331         case 0: st->d->channel_map[i] = EBUR128_LEFT; break;
332         case 1: st->d->channel_map[i] = EBUR128_RIGHT; break;
333         case 2: st->d->channel_map[i] = EBUR128_CENTER; break;
334         case 3: st->d->channel_map[i] = EBUR128_UNUSED; break;
335         case 4: st->d->channel_map[i] = EBUR128_LEFT_SURROUND; break;
336         case 5: st->d->channel_map[i] = EBUR128_RIGHT_SURROUND; break;
337         default: st->d->channel_map[i] = EBUR128_UNUSED; break;
338       }
339     }
340   }
341   return EBUR128_SUCCESS;
342 }
343 
ebur128_init_resampler(ebur128_state * st)344 static int ebur128_init_resampler(ebur128_state* st) {
345   int errcode = EBUR128_SUCCESS;
346 
347   if (st->samplerate < 96000) {
348     st->d->interp = interp_create(49, 4, st->channels);
349     CHECK_ERROR(!st->d->interp, EBUR128_ERROR_NOMEM, exit)
350   } else if (st->samplerate < 192000) {
351     st->d->interp = interp_create(49, 2, st->channels);
352     CHECK_ERROR(!st->d->interp, EBUR128_ERROR_NOMEM, exit)
353   } else {
354     st->d->resampler_buffer_input = NULL;
355     st->d->resampler_buffer_output = NULL;
356     st->d->interp = NULL;
357     goto exit;
358   }
359 
360   st->d->resampler_buffer_input_frames = st->d->samples_in_100ms * 4;
361   st->d->resampler_buffer_input = (float*) malloc(
362       st->d->resampler_buffer_input_frames * st->channels * sizeof(float));
363   CHECK_ERROR(!st->d->resampler_buffer_input, EBUR128_ERROR_NOMEM, free_interp)
364 
365   st->d->resampler_buffer_output_frames =
366       st->d->resampler_buffer_input_frames * st->d->interp->factor;
367   st->d->resampler_buffer_output = (float*) malloc(
368       st->d->resampler_buffer_output_frames * st->channels * sizeof(float));
369   CHECK_ERROR(!st->d->resampler_buffer_output, EBUR128_ERROR_NOMEM, free_input)
370 
371   return errcode;
372 
373 free_interp:
374   interp_destroy(st->d->interp);
375   st->d->interp = NULL;
376 free_input:
377   free(st->d->resampler_buffer_input);
378   st->d->resampler_buffer_input = NULL;
379 exit:
380   return errcode;
381 }
382 
ebur128_destroy_resampler(ebur128_state * st)383 static void ebur128_destroy_resampler(ebur128_state* st) {
384   free(st->d->resampler_buffer_input);
385   st->d->resampler_buffer_input = NULL;
386   free(st->d->resampler_buffer_output);
387   st->d->resampler_buffer_output = NULL;
388   interp_destroy(st->d->interp);
389   st->d->interp = NULL;
390 }
391 
ebur128_get_version(int * major,int * minor,int * patch)392 void ebur128_get_version(int* major, int* minor, int* patch) {
393   *major = EBUR128_VERSION_MAJOR;
394   *minor = EBUR128_VERSION_MINOR;
395   *patch = EBUR128_VERSION_PATCH;
396 }
397 
398 #define VALIDATE_MAX_CHANNELS (64)
399 #define VALIDATE_MAX_SAMPLERATE (2822400)
400 
401 #define VALIDATE_CHANNELS_AND_SAMPLERATE(err)                                  \
402   do {                                                                         \
403     if (channels == 0 || channels > VALIDATE_MAX_CHANNELS) {                   \
404       return (err);                                                            \
405     }                                                                          \
406                                                                                \
407     if (samplerate < 16 || samplerate > VALIDATE_MAX_SAMPLERATE) {             \
408       return (err);                                                            \
409     }                                                                          \
410   } while (0);
411 
412 ebur128_state*
ebur128_init(unsigned int channels,unsigned long samplerate,int mode)413 ebur128_init(unsigned int channels, unsigned long samplerate, int mode) {
414   int result;
415   int errcode;
416   ebur128_state* st;
417   unsigned int i;
418   size_t j;
419 
420   VALIDATE_CHANNELS_AND_SAMPLERATE(NULL);
421 
422   st = (ebur128_state*) malloc(sizeof(ebur128_state));
423   CHECK_ERROR(!st, 0, exit)
424   st->d = (struct ebur128_state_internal*) malloc(
425       sizeof(struct ebur128_state_internal));
426   CHECK_ERROR(!st->d, 0, free_state)
427   st->channels = channels;
428   errcode = ebur128_init_channel_map(st);
429   CHECK_ERROR(errcode, 0, free_internal)
430 
431   st->d->sample_peak = (double*) malloc(channels * sizeof(double));
432   CHECK_ERROR(!st->d->sample_peak, 0, free_channel_map)
433   st->d->prev_sample_peak = (double*) malloc(channels * sizeof(double));
434   CHECK_ERROR(!st->d->prev_sample_peak, 0, free_sample_peak)
435   st->d->true_peak = (double*) malloc(channels * sizeof(double));
436   CHECK_ERROR(!st->d->true_peak, 0, free_prev_sample_peak)
437   st->d->prev_true_peak = (double*) malloc(channels * sizeof(double));
438   CHECK_ERROR(!st->d->prev_true_peak, 0, free_true_peak)
439   for (i = 0; i < channels; ++i) {
440     st->d->sample_peak[i] = 0.0;
441     st->d->prev_sample_peak[i] = 0.0;
442     st->d->true_peak[i] = 0.0;
443     st->d->prev_true_peak[i] = 0.0;
444   }
445 
446   st->d->use_histogram = mode & EBUR128_MODE_HISTOGRAM ? 1 : 0;
447   st->d->history = ULONG_MAX;
448   st->samplerate = samplerate;
449   st->d->samples_in_100ms = (st->samplerate + 5) / 10;
450   st->mode = mode;
451   if ((mode & EBUR128_MODE_S) == EBUR128_MODE_S) {
452     st->d->window = 3000;
453   } else if ((mode & EBUR128_MODE_M) == EBUR128_MODE_M) {
454     st->d->window = 400;
455   } else {
456     goto free_prev_true_peak;
457   }
458   st->d->audio_data_frames = st->samplerate * st->d->window / 1000;
459   if (st->d->audio_data_frames % st->d->samples_in_100ms) {
460     /* round up to multiple of samples_in_100ms */
461     st->d->audio_data_frames =
462         (st->d->audio_data_frames + st->d->samples_in_100ms) -
463         (st->d->audio_data_frames % st->d->samples_in_100ms);
464   }
465   st->d->audio_data = (double*) malloc(st->d->audio_data_frames * st->channels *
466                                        sizeof(double));
467   CHECK_ERROR(!st->d->audio_data, 0, free_prev_true_peak)
468   for (j = 0; j < st->d->audio_data_frames * st->channels; ++j) {
469     st->d->audio_data[j] = 0.0;
470   }
471 
472   errcode = ebur128_init_filter(st);
473   CHECK_ERROR(errcode, 0, free_audio_data)
474 
475   if (st->d->use_histogram) {
476     st->d->block_energy_histogram =
477         (unsigned long*) malloc(1000 * sizeof(unsigned long));
478     CHECK_ERROR(!st->d->block_energy_histogram, 0, free_filter)
479     for (i = 0; i < 1000; ++i) {
480       st->d->block_energy_histogram[i] = 0;
481     }
482   } else {
483     st->d->block_energy_histogram = NULL;
484   }
485   if (st->d->use_histogram) {
486     st->d->short_term_block_energy_histogram =
487         (unsigned long*) malloc(1000 * sizeof(unsigned long));
488     CHECK_ERROR(!st->d->short_term_block_energy_histogram, 0,
489                 free_block_energy_histogram)
490     for (i = 0; i < 1000; ++i) {
491       st->d->short_term_block_energy_histogram[i] = 0;
492     }
493   } else {
494     st->d->short_term_block_energy_histogram = NULL;
495   }
496   STAILQ_INIT(&st->d->block_list);
497   st->d->block_list_size = 0;
498   st->d->block_list_max = st->d->history / 100;
499   STAILQ_INIT(&st->d->short_term_block_list);
500   st->d->st_block_list_size = 0;
501   st->d->st_block_list_max = st->d->history / 3000;
502   st->d->short_term_frame_counter = 0;
503 
504   result = ebur128_init_resampler(st);
505   CHECK_ERROR(result, 0, free_short_term_block_energy_histogram)
506 
507   /* the first block needs 400ms of audio data */
508   st->d->needed_frames = st->d->samples_in_100ms * 4;
509   /* start at the beginning of the buffer */
510   st->d->audio_data_index = 0;
511 
512   /* initialize static constants */
513   relative_gate_factor = pow(10.0, relative_gate / 10.0);
514   minus_twenty_decibels = pow(10.0, -20.0 / 10.0);
515   histogram_energy_boundaries[0] = pow(10.0, (-70.0 + 0.691) / 10.0);
516   if (st->d->use_histogram) {
517     for (i = 0; i < 1000; ++i) {
518       histogram_energies[i] =
519           pow(10.0, ((double) i / 10.0 - 69.95 + 0.691) / 10.0);
520     }
521     for (i = 1; i < 1001; ++i) {
522       histogram_energy_boundaries[i] =
523           pow(10.0, ((double) i / 10.0 - 70.0 + 0.691) / 10.0);
524     }
525   }
526 
527   return st;
528 
529 free_short_term_block_energy_histogram:
530   free(st->d->short_term_block_energy_histogram);
531 free_block_energy_histogram:
532   free(st->d->block_energy_histogram);
533 free_filter:
534   free(st->d->v);
535 free_audio_data:
536   free(st->d->audio_data);
537 free_prev_true_peak:
538   free(st->d->prev_true_peak);
539 free_true_peak:
540   free(st->d->true_peak);
541 free_prev_sample_peak:
542   free(st->d->prev_sample_peak);
543 free_sample_peak:
544   free(st->d->sample_peak);
545 free_channel_map:
546   free(st->d->channel_map);
547 free_internal:
548   free(st->d);
549 free_state:
550   free(st);
551 exit:
552   return NULL;
553 }
554 
ebur128_destroy(ebur128_state ** st)555 void ebur128_destroy(ebur128_state** st) {
556   struct ebur128_dq_entry* entry;
557   free((*st)->d->short_term_block_energy_histogram);
558   free((*st)->d->block_energy_histogram);
559   free((*st)->d->v);
560   free((*st)->d->audio_data);
561   free((*st)->d->channel_map);
562   free((*st)->d->sample_peak);
563   free((*st)->d->prev_sample_peak);
564   free((*st)->d->true_peak);
565   free((*st)->d->prev_true_peak);
566   while (!STAILQ_EMPTY(&(*st)->d->block_list)) {
567     entry = STAILQ_FIRST(&(*st)->d->block_list);
568     STAILQ_REMOVE_HEAD(&(*st)->d->block_list, entries);
569     free(entry);
570   }
571   while (!STAILQ_EMPTY(&(*st)->d->short_term_block_list)) {
572     entry = STAILQ_FIRST(&(*st)->d->short_term_block_list);
573     STAILQ_REMOVE_HEAD(&(*st)->d->short_term_block_list, entries);
574     free(entry);
575   }
576   ebur128_destroy_resampler(*st);
577   free((*st)->d);
578   free(*st);
579   *st = NULL;
580 }
581 
ebur128_check_true_peak(ebur128_state * st,size_t frames)582 static void ebur128_check_true_peak(ebur128_state* st, size_t frames) {
583   size_t c, i, frames_out;
584 
585   frames_out =
586       interp_process(st->d->interp, frames, st->d->resampler_buffer_input,
587                      st->d->resampler_buffer_output);
588 
589   for (i = 0; i < frames_out; ++i) {
590     for (c = 0; c < st->channels; ++c) {
591       double val =
592           (double) st->d->resampler_buffer_output[i * st->channels + c];
593 
594       if (EBUR128_MAX(val, -val) > st->d->prev_true_peak[c]) {
595         st->d->prev_true_peak[c] = EBUR128_MAX(val, -val);
596       }
597     }
598   }
599 }
600 
601 #if defined(__SSE2_MATH__) || defined(_M_X64) || _M_IX86_FP >= 2
602 #include <xmmintrin.h>
603 #define TURN_ON_FTZ                                                            \
604   unsigned int mxcsr = _mm_getcsr();                                           \
605   _mm_setcsr(mxcsr | _MM_FLUSH_ZERO_ON);
606 #define TURN_OFF_FTZ _mm_setcsr(mxcsr);
607 #define FLUSH_MANUALLY
608 #else
609 #warning "manual FTZ is being used, please enable SSE2 (-msse2 -mfpmath=sse)"
610 #define TURN_ON_FTZ
611 #define TURN_OFF_FTZ
612 #define FLUSH_MANUALLY                                                         \
613   st->d->v[c][4] = fabs(st->d->v[c][4]) < DBL_MIN ? 0.0 : st->d->v[c][4];      \
614   st->d->v[c][3] = fabs(st->d->v[c][3]) < DBL_MIN ? 0.0 : st->d->v[c][3];      \
615   st->d->v[c][2] = fabs(st->d->v[c][2]) < DBL_MIN ? 0.0 : st->d->v[c][2];      \
616   st->d->v[c][1] = fabs(st->d->v[c][1]) < DBL_MIN ? 0.0 : st->d->v[c][1];
617 #endif
618 
619 #define EBUR128_FILTER(type, min_scale, max_scale)                             \
620   static void ebur128_filter_##type(ebur128_state* st, const type* src,        \
621                                     size_t frames) {                           \
622     static double scaling_factor =                                             \
623         EBUR128_MAX(-((double) (min_scale)), (double) (max_scale));            \
624                                                                                \
625     double* audio_data = st->d->audio_data + st->d->audio_data_index;          \
626     size_t i, c;                                                               \
627                                                                                \
628     TURN_ON_FTZ                                                                \
629                                                                                \
630     if ((st->mode & EBUR128_MODE_SAMPLE_PEAK) == EBUR128_MODE_SAMPLE_PEAK) {   \
631       for (c = 0; c < st->channels; ++c) {                                     \
632         double max = 0.0;                                                      \
633         for (i = 0; i < frames; ++i) {                                         \
634           double cur = (double) src[i * st->channels + c];                     \
635           if (EBUR128_MAX(cur, -cur) > max) {                                  \
636             max = EBUR128_MAX(cur, -cur);                                      \
637           }                                                                    \
638         }                                                                      \
639         max /= scaling_factor;                                                 \
640         if (max > st->d->prev_sample_peak[c]) {                                \
641           st->d->prev_sample_peak[c] = max;                                    \
642         }                                                                      \
643       }                                                                        \
644     }                                                                          \
645     if ((st->mode & EBUR128_MODE_TRUE_PEAK) == EBUR128_MODE_TRUE_PEAK &&       \
646         st->d->interp) {                                                       \
647       for (i = 0; i < frames; ++i) {                                           \
648         for (c = 0; c < st->channels; ++c) {                                   \
649           st->d->resampler_buffer_input[i * st->channels + c] =                \
650               (float) ((double) src[i * st->channels + c] / scaling_factor);   \
651         }                                                                      \
652       }                                                                        \
653       ebur128_check_true_peak(st, frames);                                     \
654     }                                                                          \
655     for (c = 0; c < st->channels; ++c) {                                       \
656       if (st->d->channel_map[c] == EBUR128_UNUSED) {                           \
657         continue;                                                              \
658       }                                                                        \
659       for (i = 0; i < frames; ++i) {                                           \
660         st->d->v[c][0] =                                                       \
661             (double) ((double) src[i * st->channels + c] / scaling_factor) -   \
662             st->d->a[1] * st->d->v[c][1] - /**/                                \
663             st->d->a[2] * st->d->v[c][2] - /**/                                \
664             st->d->a[3] * st->d->v[c][3] - /**/                                \
665             st->d->a[4] * st->d->v[c][4];                                      \
666         audio_data[i * st->channels + c] = /**/                                \
667             st->d->b[0] * st->d->v[c][0] + /**/                                \
668             st->d->b[1] * st->d->v[c][1] + /**/                                \
669             st->d->b[2] * st->d->v[c][2] + /**/                                \
670             st->d->b[3] * st->d->v[c][3] + /**/                                \
671             st->d->b[4] * st->d->v[c][4];                                      \
672         st->d->v[c][4] = st->d->v[c][3];                                       \
673         st->d->v[c][3] = st->d->v[c][2];                                       \
674         st->d->v[c][2] = st->d->v[c][1];                                       \
675         st->d->v[c][1] = st->d->v[c][0];                                       \
676       }                                                                        \
677       FLUSH_MANUALLY                                                           \
678     }                                                                          \
679     TURN_OFF_FTZ                                                               \
680   }
681 
EBUR128_FILTER(short,SHRT_MIN,SHRT_MAX)682 EBUR128_FILTER(short, SHRT_MIN, SHRT_MAX)
683 EBUR128_FILTER(int, INT_MIN, INT_MAX)
684 EBUR128_FILTER(float, -1.0f, 1.0f)
685 EBUR128_FILTER(double, -1.0, 1.0)
686 
687 static double ebur128_energy_to_loudness(double energy) {
688   return 10 * (log(energy) / log(10.0)) - 0.691;
689 }
690 
find_histogram_index(double energy)691 static size_t find_histogram_index(double energy) {
692   size_t index_min = 0;
693   size_t index_max = 1000;
694   size_t index_mid;
695 
696   do {
697     index_mid = (index_min + index_max) / 2;
698     if (energy >= histogram_energy_boundaries[index_mid]) {
699       index_min = index_mid;
700     } else {
701       index_max = index_mid;
702     }
703   } while (index_max - index_min != 1);
704 
705   return index_min;
706 }
707 
ebur128_calc_gating_block(ebur128_state * st,size_t frames_per_block,double * optional_output)708 static int ebur128_calc_gating_block(ebur128_state* st,
709                                      size_t frames_per_block,
710                                      double* optional_output) {
711   size_t i, c;
712   double sum = 0.0;
713   double channel_sum;
714   for (c = 0; c < st->channels; ++c) {
715     if (st->d->channel_map[c] == EBUR128_UNUSED) {
716       continue;
717     }
718     channel_sum = 0.0;
719     if (st->d->audio_data_index < frames_per_block * st->channels) {
720       for (i = 0; i < st->d->audio_data_index / st->channels; ++i) {
721         channel_sum += st->d->audio_data[i * st->channels + c] *
722                        st->d->audio_data[i * st->channels + c];
723       }
724       for (i = st->d->audio_data_frames -
725                (frames_per_block - st->d->audio_data_index / st->channels);
726            i < st->d->audio_data_frames; ++i) {
727         channel_sum += st->d->audio_data[i * st->channels + c] *
728                        st->d->audio_data[i * st->channels + c];
729       }
730     } else {
731       for (i = st->d->audio_data_index / st->channels - frames_per_block;
732            i < st->d->audio_data_index / st->channels; ++i) {
733         channel_sum += st->d->audio_data[i * st->channels + c] *
734                        st->d->audio_data[i * st->channels + c];
735       }
736     }
737     if (st->d->channel_map[c] == EBUR128_Mp110 ||
738         st->d->channel_map[c] == EBUR128_Mm110 ||
739         st->d->channel_map[c] == EBUR128_Mp060 ||
740         st->d->channel_map[c] == EBUR128_Mm060 ||
741         st->d->channel_map[c] == EBUR128_Mp090 ||
742         st->d->channel_map[c] == EBUR128_Mm090) {
743       channel_sum *= 1.41;
744     } else if (st->d->channel_map[c] == EBUR128_DUAL_MONO) {
745       channel_sum *= 2.0;
746     }
747     sum += channel_sum;
748   }
749 
750   sum /= (double) frames_per_block;
751 
752   if (optional_output) {
753     *optional_output = sum;
754     return EBUR128_SUCCESS;
755   }
756 
757   if (sum >= histogram_energy_boundaries[0]) {
758     if (st->d->use_histogram) {
759       ++st->d->block_energy_histogram[find_histogram_index(sum)];
760     } else {
761       struct ebur128_dq_entry* block;
762       if (st->d->block_list_size == st->d->block_list_max) {
763         block = STAILQ_FIRST(&st->d->block_list);
764         STAILQ_REMOVE_HEAD(&st->d->block_list, entries);
765       } else {
766         block =
767             (struct ebur128_dq_entry*) malloc(sizeof(struct ebur128_dq_entry));
768         if (!block) {
769           return EBUR128_ERROR_NOMEM;
770         }
771         st->d->block_list_size++;
772       }
773       block->z = sum;
774       STAILQ_INSERT_TAIL(&st->d->block_list, block, entries);
775     }
776   }
777 
778   return EBUR128_SUCCESS;
779 }
780 
ebur128_set_channel(ebur128_state * st,unsigned int channel_number,int value)781 int ebur128_set_channel(ebur128_state* st,
782                         unsigned int channel_number,
783                         int value) {
784   if (channel_number >= st->channels) {
785     return EBUR128_ERROR_INVALID_CHANNEL_INDEX;
786   }
787   if (value == EBUR128_DUAL_MONO &&
788       (st->channels != 1 || channel_number != 0)) {
789     fprintf(stderr, "EBUR128_DUAL_MONO only works with mono files!\n");
790     return EBUR128_ERROR_INVALID_CHANNEL_INDEX;
791   }
792   st->d->channel_map[channel_number] = value;
793   return EBUR128_SUCCESS;
794 }
795 
ebur128_change_parameters(ebur128_state * st,unsigned int channels,unsigned long samplerate)796 int ebur128_change_parameters(ebur128_state* st,
797                               unsigned int channels,
798                               unsigned long samplerate) {
799   int errcode = EBUR128_SUCCESS;
800   size_t j;
801 
802   /* This is needed to suppress a clang-tidy warning. */
803 #ifndef __has_builtin
804 #define __has_builtin(x) 0
805 #endif
806 #if __has_builtin(__builtin_unreachable)
807   if (st->channels == 0) {
808     __builtin_unreachable();
809   }
810 #endif
811 
812   VALIDATE_CHANNELS_AND_SAMPLERATE(EBUR128_ERROR_NOMEM);
813 
814   if (channels == st->channels && samplerate == st->samplerate) {
815     return EBUR128_ERROR_NO_CHANGE;
816   }
817 
818   free(st->d->audio_data);
819   st->d->audio_data = NULL;
820 
821   if (channels != st->channels) {
822     unsigned int i;
823 
824     free(st->d->channel_map);
825     st->d->channel_map = NULL;
826     free(st->d->sample_peak);
827     st->d->sample_peak = NULL;
828     free(st->d->prev_sample_peak);
829     st->d->prev_sample_peak = NULL;
830     free(st->d->true_peak);
831     st->d->true_peak = NULL;
832     free(st->d->prev_true_peak);
833     st->d->prev_true_peak = NULL;
834     st->channels = channels;
835 
836     errcode = ebur128_init_channel_map(st);
837     CHECK_ERROR(errcode, EBUR128_ERROR_NOMEM, exit)
838 
839     st->d->sample_peak = (double*) malloc(channels * sizeof(double));
840     CHECK_ERROR(!st->d->sample_peak, EBUR128_ERROR_NOMEM, exit)
841     st->d->prev_sample_peak = (double*) malloc(channels * sizeof(double));
842     CHECK_ERROR(!st->d->prev_sample_peak, EBUR128_ERROR_NOMEM, exit)
843     st->d->true_peak = (double*) malloc(channels * sizeof(double));
844     CHECK_ERROR(!st->d->true_peak, EBUR128_ERROR_NOMEM, exit)
845     st->d->prev_true_peak = (double*) malloc(channels * sizeof(double));
846     CHECK_ERROR(!st->d->prev_true_peak, EBUR128_ERROR_NOMEM, exit)
847     for (i = 0; i < channels; ++i) {
848       st->d->sample_peak[i] = 0.0;
849       st->d->prev_sample_peak[i] = 0.0;
850       st->d->true_peak[i] = 0.0;
851       st->d->prev_true_peak[i] = 0.0;
852     }
853   }
854   if (samplerate != st->samplerate) {
855     st->samplerate = samplerate;
856     st->d->samples_in_100ms = (st->samplerate + 5) / 10;
857   }
858 
859   /* If we're here, either samplerate or channels
860    * have changed. Re-init filter. */
861   free(st->d->v);
862   st->d->v = NULL;
863   errcode = ebur128_init_filter(st);
864   CHECK_ERROR(errcode, EBUR128_ERROR_NOMEM, exit)
865 
866   st->d->audio_data_frames = st->samplerate * st->d->window / 1000;
867   if (st->d->audio_data_frames % st->d->samples_in_100ms) {
868     /* round up to multiple of samples_in_100ms */
869     st->d->audio_data_frames =
870         (st->d->audio_data_frames + st->d->samples_in_100ms) -
871         (st->d->audio_data_frames % st->d->samples_in_100ms);
872   }
873   st->d->audio_data = (double*) malloc(st->d->audio_data_frames * st->channels *
874                                        sizeof(double));
875   CHECK_ERROR(!st->d->audio_data, EBUR128_ERROR_NOMEM, exit)
876   for (j = 0; j < st->d->audio_data_frames * st->channels; ++j) {
877     st->d->audio_data[j] = 0.0;
878   }
879 
880   ebur128_destroy_resampler(st);
881   errcode = ebur128_init_resampler(st);
882   CHECK_ERROR(errcode, EBUR128_ERROR_NOMEM, exit)
883 
884   /* the first block needs 400ms of audio data */
885   st->d->needed_frames = st->d->samples_in_100ms * 4;
886   /* start at the beginning of the buffer */
887   st->d->audio_data_index = 0;
888   /* reset short term frame counter */
889   st->d->short_term_frame_counter = 0;
890 
891 exit:
892   return errcode;
893 }
894 
ebur128_set_max_window(ebur128_state * st,unsigned long window)895 int ebur128_set_max_window(ebur128_state* st, unsigned long window) {
896   int errcode = EBUR128_SUCCESS;
897   size_t j;
898 
899   if ((st->mode & EBUR128_MODE_S) == EBUR128_MODE_S && window < 3000) {
900     window = 3000;
901   } else if ((st->mode & EBUR128_MODE_M) == EBUR128_MODE_M && window < 400) {
902     window = 400;
903   }
904 
905   if (window == st->d->window) {
906     return EBUR128_ERROR_NO_CHANGE;
907   }
908 
909   size_t new_audio_data_frames;
910   if (safe_size_mul(st->samplerate, window, &new_audio_data_frames) != 0 ||
911       new_audio_data_frames > ((size_t) -1) - st->d->samples_in_100ms) {
912     return EBUR128_ERROR_NOMEM;
913   }
914   if (new_audio_data_frames % st->d->samples_in_100ms) {
915     /* round up to multiple of samples_in_100ms */
916     new_audio_data_frames = (new_audio_data_frames + st->d->samples_in_100ms) -
917                             (new_audio_data_frames % st->d->samples_in_100ms);
918   }
919 
920   size_t new_audio_data_size;
921   if (safe_size_mul(new_audio_data_frames, st->channels * sizeof(double),
922                     &new_audio_data_size) != 0) {
923     return EBUR128_ERROR_NOMEM;
924   }
925 
926   double* new_audio_data = (double*) malloc(new_audio_data_size);
927   CHECK_ERROR(!new_audio_data, EBUR128_ERROR_NOMEM, exit)
928 
929   st->d->window = window;
930   free(st->d->audio_data);
931   st->d->audio_data = new_audio_data;
932   st->d->audio_data_frames = new_audio_data_frames;
933   for (j = 0; j < st->d->audio_data_frames * st->channels; ++j) {
934     st->d->audio_data[j] = 0.0;
935   }
936 
937   /* the first block needs 400ms of audio data */
938   st->d->needed_frames = st->d->samples_in_100ms * 4;
939   /* start at the beginning of the buffer */
940   st->d->audio_data_index = 0;
941   /* reset short term frame counter */
942   st->d->short_term_frame_counter = 0;
943 
944 exit:
945   return errcode;
946 }
947 
ebur128_set_max_history(ebur128_state * st,unsigned long history)948 int ebur128_set_max_history(ebur128_state* st, unsigned long history) {
949   if ((st->mode & EBUR128_MODE_LRA) == EBUR128_MODE_LRA && history < 3000) {
950     history = 3000;
951   } else if ((st->mode & EBUR128_MODE_M) == EBUR128_MODE_M && history < 400) {
952     history = 400;
953   }
954   if (history == st->d->history) {
955     return EBUR128_ERROR_NO_CHANGE;
956   }
957   st->d->history = history;
958   st->d->block_list_max = st->d->history / 100;
959   st->d->st_block_list_max = st->d->history / 3000;
960   while (st->d->block_list_size > st->d->block_list_max) {
961     struct ebur128_dq_entry* block = STAILQ_FIRST(&st->d->block_list);
962     STAILQ_REMOVE_HEAD(&st->d->block_list, entries);
963     free(block);
964     st->d->block_list_size--;
965   }
966   while (st->d->st_block_list_size > st->d->st_block_list_max) {
967     struct ebur128_dq_entry* block =
968         STAILQ_FIRST(&st->d->short_term_block_list);
969     STAILQ_REMOVE_HEAD(&st->d->short_term_block_list, entries);
970     free(block);
971     st->d->st_block_list_size--;
972   }
973   return EBUR128_SUCCESS;
974 }
975 
976 static int ebur128_energy_shortterm(ebur128_state* st, double* out);
977 #define EBUR128_ADD_FRAMES(type)                                               \
978   int ebur128_add_frames_##type(ebur128_state* st, const type* src,            \
979                                 size_t frames) {                               \
980     size_t src_index = 0;                                                      \
981     unsigned int c = 0;                                                        \
982     for (c = 0; c < st->channels; c++) {                                       \
983       st->d->prev_sample_peak[c] = 0.0;                                        \
984       st->d->prev_true_peak[c] = 0.0;                                          \
985     }                                                                          \
986     while (frames > 0) {                                                       \
987       if (frames >= st->d->needed_frames) {                                    \
988         ebur128_filter_##type(st, src + src_index, st->d->needed_frames);      \
989         src_index += st->d->needed_frames * st->channels;                      \
990         frames -= st->d->needed_frames;                                        \
991         st->d->audio_data_index += st->d->needed_frames * st->channels;        \
992         /* calculate the new gating block */                                   \
993         if ((st->mode & EBUR128_MODE_I) == EBUR128_MODE_I) {                   \
994           if (ebur128_calc_gating_block(st, st->d->samples_in_100ms * 4,       \
995                                         NULL)) {                               \
996             return EBUR128_ERROR_NOMEM;                                        \
997           }                                                                    \
998         }                                                                      \
999         if ((st->mode & EBUR128_MODE_LRA) == EBUR128_MODE_LRA) {               \
1000           st->d->short_term_frame_counter += st->d->needed_frames;             \
1001           if (st->d->short_term_frame_counter ==                               \
1002               st->d->samples_in_100ms * 30) {                                  \
1003             struct ebur128_dq_entry* block;                                    \
1004             double st_energy;                                                  \
1005             if (ebur128_energy_shortterm(st, &st_energy) == EBUR128_SUCCESS && \
1006                 st_energy >= histogram_energy_boundaries[0]) {                 \
1007               if (st->d->use_histogram) {                                      \
1008                 ++st->d->short_term_block_energy_histogram                     \
1009                       [find_histogram_index(st_energy)];                       \
1010               } else {                                                         \
1011                 if (st->d->st_block_list_size == st->d->st_block_list_max) {   \
1012                   block = STAILQ_FIRST(&st->d->short_term_block_list);         \
1013                   STAILQ_REMOVE_HEAD(&st->d->short_term_block_list, entries);  \
1014                 } else {                                                       \
1015                   block = (struct ebur128_dq_entry*) malloc(                   \
1016                       sizeof(struct ebur128_dq_entry));                        \
1017                   if (!block) {                                                \
1018                     return EBUR128_ERROR_NOMEM;                                \
1019                   }                                                            \
1020                   st->d->st_block_list_size++;                                 \
1021                 }                                                              \
1022                 block->z = st_energy;                                          \
1023                 STAILQ_INSERT_TAIL(&st->d->short_term_block_list, block,       \
1024                                    entries);                                   \
1025               }                                                                \
1026             }                                                                  \
1027             st->d->short_term_frame_counter = st->d->samples_in_100ms * 20;    \
1028           }                                                                    \
1029         }                                                                      \
1030         /* 100ms are needed for all blocks besides the first one */            \
1031         st->d->needed_frames = st->d->samples_in_100ms;                        \
1032         /* reset audio_data_index when buffer full */                          \
1033         if (st->d->audio_data_index ==                                         \
1034             st->d->audio_data_frames * st->channels) {                         \
1035           st->d->audio_data_index = 0;                                         \
1036         }                                                                      \
1037       } else {                                                                 \
1038         ebur128_filter_##type(st, src + src_index, frames);                    \
1039         st->d->audio_data_index += frames * st->channels;                      \
1040         if ((st->mode & EBUR128_MODE_LRA) == EBUR128_MODE_LRA) {               \
1041           st->d->short_term_frame_counter += frames;                           \
1042         }                                                                      \
1043         st->d->needed_frames -= (unsigned long) frames;                        \
1044         frames = 0;                                                            \
1045       }                                                                        \
1046     }                                                                          \
1047     for (c = 0; c < st->channels; c++) {                                       \
1048       if (st->d->prev_sample_peak[c] > st->d->sample_peak[c]) {                \
1049         st->d->sample_peak[c] = st->d->prev_sample_peak[c];                    \
1050       }                                                                        \
1051       if (st->d->prev_true_peak[c] > st->d->true_peak[c]) {                    \
1052         st->d->true_peak[c] = st->d->prev_true_peak[c];                        \
1053       }                                                                        \
1054     }                                                                          \
1055     return EBUR128_SUCCESS;                                                    \
1056   }
1057 
1058 EBUR128_ADD_FRAMES(short)
EBUR128_ADD_FRAMES(int)1059 EBUR128_ADD_FRAMES(int)
1060 EBUR128_ADD_FRAMES(float)
1061 EBUR128_ADD_FRAMES(double)
1062 
1063 static int ebur128_calc_relative_threshold(ebur128_state* st,
1064                                            size_t* above_thresh_counter,
1065                                            double* relative_threshold) {
1066   struct ebur128_dq_entry* it;
1067   size_t i;
1068 
1069   if (st->d->use_histogram) {
1070     for (i = 0; i < 1000; ++i) {
1071       *relative_threshold +=
1072           st->d->block_energy_histogram[i] * histogram_energies[i];
1073       *above_thresh_counter += st->d->block_energy_histogram[i];
1074     }
1075   } else {
1076     STAILQ_FOREACH(it, &st->d->block_list, entries) {
1077       ++*above_thresh_counter;
1078       *relative_threshold += it->z;
1079     }
1080   }
1081 
1082   return EBUR128_SUCCESS;
1083 }
1084 
1085 static int
ebur128_gated_loudness(ebur128_state ** sts,size_t size,double * out)1086 ebur128_gated_loudness(ebur128_state** sts, size_t size, double* out) {
1087   struct ebur128_dq_entry* it;
1088   double gated_loudness = 0.0;
1089   double relative_threshold = 0.0;
1090   size_t above_thresh_counter = 0;
1091   size_t i, j, start_index;
1092 
1093   for (i = 0; i < size; i++) {
1094     if (sts[i] && (sts[i]->mode & EBUR128_MODE_I) != EBUR128_MODE_I) {
1095       return EBUR128_ERROR_INVALID_MODE;
1096     }
1097   }
1098 
1099   for (i = 0; i < size; i++) {
1100     if (!sts[i]) {
1101       continue;
1102     }
1103     ebur128_calc_relative_threshold(sts[i], &above_thresh_counter,
1104                                     &relative_threshold);
1105   }
1106   if (!above_thresh_counter) {
1107     *out = -HUGE_VAL;
1108     return EBUR128_SUCCESS;
1109   }
1110 
1111   relative_threshold /= (double) above_thresh_counter;
1112   relative_threshold *= relative_gate_factor;
1113 
1114   above_thresh_counter = 0;
1115   if (relative_threshold < histogram_energy_boundaries[0]) {
1116     start_index = 0;
1117   } else {
1118     start_index = find_histogram_index(relative_threshold);
1119     if (relative_threshold > histogram_energies[start_index]) {
1120       ++start_index;
1121     }
1122   }
1123   for (i = 0; i < size; i++) {
1124     if (!sts[i]) {
1125       continue;
1126     }
1127     if (sts[i]->d->use_histogram) {
1128       for (j = start_index; j < 1000; ++j) {
1129         gated_loudness +=
1130             sts[i]->d->block_energy_histogram[j] * histogram_energies[j];
1131         above_thresh_counter += sts[i]->d->block_energy_histogram[j];
1132       }
1133     } else {
1134       STAILQ_FOREACH(it, &sts[i]->d->block_list, entries) {
1135         if (it->z >= relative_threshold) {
1136           ++above_thresh_counter;
1137           gated_loudness += it->z;
1138         }
1139       }
1140     }
1141   }
1142   if (!above_thresh_counter) {
1143     *out = -HUGE_VAL;
1144     return EBUR128_SUCCESS;
1145   }
1146   gated_loudness /= (double) above_thresh_counter;
1147   *out = ebur128_energy_to_loudness(gated_loudness);
1148   return EBUR128_SUCCESS;
1149 }
1150 
ebur128_relative_threshold(ebur128_state * st,double * out)1151 int ebur128_relative_threshold(ebur128_state* st, double* out) {
1152   double relative_threshold = 0.0;
1153   size_t above_thresh_counter = 0;
1154 
1155   if ((st->mode & EBUR128_MODE_I) != EBUR128_MODE_I) {
1156     return EBUR128_ERROR_INVALID_MODE;
1157   }
1158 
1159   ebur128_calc_relative_threshold(st, &above_thresh_counter,
1160                                   &relative_threshold);
1161 
1162   if (!above_thresh_counter) {
1163     *out = -70.0;
1164     return EBUR128_SUCCESS;
1165   }
1166 
1167   relative_threshold /= (double) above_thresh_counter;
1168   relative_threshold *= relative_gate_factor;
1169 
1170   *out = ebur128_energy_to_loudness(relative_threshold);
1171   return EBUR128_SUCCESS;
1172 }
1173 
ebur128_loudness_global(ebur128_state * st,double * out)1174 int ebur128_loudness_global(ebur128_state* st, double* out) {
1175   return ebur128_gated_loudness(&st, 1, out);
1176 }
1177 
ebur128_loudness_global_multiple(ebur128_state ** sts,size_t size,double * out)1178 int ebur128_loudness_global_multiple(ebur128_state** sts,
1179                                      size_t size,
1180                                      double* out) {
1181   return ebur128_gated_loudness(sts, size, out);
1182 }
1183 
ebur128_energy_in_interval(ebur128_state * st,size_t interval_frames,double * out)1184 static int ebur128_energy_in_interval(ebur128_state* st,
1185                                       size_t interval_frames,
1186                                       double* out) {
1187   if (interval_frames > st->d->audio_data_frames) {
1188     return EBUR128_ERROR_INVALID_MODE;
1189   }
1190   ebur128_calc_gating_block(st, interval_frames, out);
1191   return EBUR128_SUCCESS;
1192 }
1193 
ebur128_energy_shortterm(ebur128_state * st,double * out)1194 static int ebur128_energy_shortterm(ebur128_state* st, double* out) {
1195   return ebur128_energy_in_interval(st, st->d->samples_in_100ms * 30, out);
1196 }
1197 
ebur128_loudness_momentary(ebur128_state * st,double * out)1198 int ebur128_loudness_momentary(ebur128_state* st, double* out) {
1199   double energy;
1200   int error;
1201 
1202   error = ebur128_energy_in_interval(st, st->d->samples_in_100ms * 4, &energy);
1203   if (error) {
1204     return error;
1205   }
1206 
1207   if (energy <= 0.0) {
1208     *out = -HUGE_VAL;
1209     return EBUR128_SUCCESS;
1210   }
1211 
1212   *out = ebur128_energy_to_loudness(energy);
1213   return EBUR128_SUCCESS;
1214 }
1215 
ebur128_loudness_shortterm(ebur128_state * st,double * out)1216 int ebur128_loudness_shortterm(ebur128_state* st, double* out) {
1217   double energy;
1218   int error;
1219 
1220   error = ebur128_energy_shortterm(st, &energy);
1221   if (error) {
1222     return error;
1223   }
1224 
1225   if (energy <= 0.0) {
1226     *out = -HUGE_VAL;
1227     return EBUR128_SUCCESS;
1228   }
1229 
1230   *out = ebur128_energy_to_loudness(energy);
1231   return EBUR128_SUCCESS;
1232 }
1233 
ebur128_loudness_window(ebur128_state * st,unsigned long window,double * out)1234 int ebur128_loudness_window(ebur128_state* st,
1235                             unsigned long window,
1236                             double* out) {
1237   double energy;
1238   size_t interval_frames;
1239   int error;
1240 
1241   if (window > st->d->window) {
1242     return EBUR128_ERROR_INVALID_MODE;
1243   }
1244 
1245   interval_frames = st->samplerate * window / 1000;
1246   error = ebur128_energy_in_interval(st, interval_frames, &energy);
1247   if (error) {
1248     return error;
1249   }
1250 
1251   if (energy <= 0.0) {
1252     *out = -HUGE_VAL;
1253     return EBUR128_SUCCESS;
1254   }
1255 
1256   *out = ebur128_energy_to_loudness(energy);
1257   return EBUR128_SUCCESS;
1258 }
1259 
ebur128_double_cmp(const void * p1,const void * p2)1260 static int ebur128_double_cmp(const void* p1, const void* p2) {
1261   const double* d1 = (const double*) p1;
1262   const double* d2 = (const double*) p2;
1263   return (*d1 > *d2) - (*d1 < *d2);
1264 }
1265 
1266 /* EBU - TECH 3342 */
ebur128_loudness_range_multiple(ebur128_state ** sts,size_t size,double * out)1267 int ebur128_loudness_range_multiple(ebur128_state** sts,
1268                                     size_t size,
1269                                     double* out) {
1270   size_t i, j;
1271   struct ebur128_dq_entry* it;
1272   double* stl_vector;
1273   size_t stl_size;
1274   double* stl_relgated;
1275   size_t stl_relgated_size;
1276   double stl_power, stl_integrated;
1277   /* High and low percentile energy */
1278   double h_en, l_en;
1279   int use_histogram = 0;
1280 
1281   for (i = 0; i < size; ++i) {
1282     if (sts[i]) {
1283       if ((sts[i]->mode & EBUR128_MODE_LRA) != EBUR128_MODE_LRA) {
1284         return EBUR128_ERROR_INVALID_MODE;
1285       }
1286       if (i == 0 && sts[i]->mode & EBUR128_MODE_HISTOGRAM) {
1287         use_histogram = 1;
1288       } else if (use_histogram != !!(sts[i]->mode & EBUR128_MODE_HISTOGRAM)) {
1289         return EBUR128_ERROR_INVALID_MODE;
1290       }
1291     }
1292   }
1293 
1294   if (use_histogram) {
1295     unsigned long hist[1000] = { 0 };
1296     size_t percentile_low, percentile_high;
1297     size_t index;
1298 
1299     stl_size = 0;
1300     stl_power = 0.0;
1301     for (i = 0; i < size; ++i) {
1302       if (!sts[i]) {
1303         continue;
1304       }
1305       for (j = 0; j < 1000; ++j) {
1306         hist[j] += sts[i]->d->short_term_block_energy_histogram[j];
1307         stl_size += sts[i]->d->short_term_block_energy_histogram[j];
1308         stl_power += sts[i]->d->short_term_block_energy_histogram[j] *
1309                      histogram_energies[j];
1310       }
1311     }
1312     if (!stl_size) {
1313       *out = 0.0;
1314       return EBUR128_SUCCESS;
1315     }
1316 
1317     stl_power /= stl_size;
1318     stl_integrated = minus_twenty_decibels * stl_power;
1319 
1320     if (stl_integrated < histogram_energy_boundaries[0]) {
1321       index = 0;
1322     } else {
1323       index = find_histogram_index(stl_integrated);
1324       if (stl_integrated > histogram_energies[index]) {
1325         ++index;
1326       }
1327     }
1328     stl_size = 0;
1329     for (j = index; j < 1000; ++j) {
1330       stl_size += hist[j];
1331     }
1332     if (!stl_size) {
1333       *out = 0.0;
1334       return EBUR128_SUCCESS;
1335     }
1336 
1337     percentile_low = (size_t) ((stl_size - 1) * 0.1 + 0.5);
1338     percentile_high = (size_t) ((stl_size - 1) * 0.95 + 0.5);
1339 
1340     stl_size = 0;
1341     j = index;
1342     while (stl_size <= percentile_low) {
1343       stl_size += hist[j++];
1344     }
1345     l_en = histogram_energies[j - 1];
1346     while (stl_size <= percentile_high) {
1347       stl_size += hist[j++];
1348     }
1349     h_en = histogram_energies[j - 1];
1350 
1351     *out = ebur128_energy_to_loudness(h_en) - ebur128_energy_to_loudness(l_en);
1352     return EBUR128_SUCCESS;
1353   }
1354 
1355   stl_size = 0;
1356   for (i = 0; i < size; ++i) {
1357     if (!sts[i]) {
1358       continue;
1359     }
1360     STAILQ_FOREACH(it, &sts[i]->d->short_term_block_list, entries) {
1361       ++stl_size;
1362     }
1363   }
1364   if (!stl_size) {
1365     *out = 0.0;
1366     return EBUR128_SUCCESS;
1367   }
1368   stl_vector = (double*) malloc(stl_size * sizeof(double));
1369   if (!stl_vector) {
1370     return EBUR128_ERROR_NOMEM;
1371   }
1372 
1373   j = 0;
1374   for (i = 0; i < size; ++i) {
1375     if (!sts[i]) {
1376       continue;
1377     }
1378     STAILQ_FOREACH(it, &sts[i]->d->short_term_block_list, entries) {
1379       stl_vector[j] = it->z;
1380       ++j;
1381     }
1382   }
1383   qsort(stl_vector, stl_size, sizeof(double), ebur128_double_cmp);
1384   stl_power = 0.0;
1385   for (i = 0; i < stl_size; ++i) {
1386     stl_power += stl_vector[i];
1387   }
1388   stl_power /= (double) stl_size;
1389   stl_integrated = minus_twenty_decibels * stl_power;
1390 
1391   stl_relgated = stl_vector;
1392   stl_relgated_size = stl_size;
1393   while (stl_relgated_size > 0 && *stl_relgated < stl_integrated) {
1394     ++stl_relgated;
1395     --stl_relgated_size;
1396   }
1397 
1398   if (stl_relgated_size) {
1399     h_en = stl_relgated[(size_t) ((stl_relgated_size - 1) * 0.95 + 0.5)];
1400     l_en = stl_relgated[(size_t) ((stl_relgated_size - 1) * 0.1 + 0.5)];
1401     free(stl_vector);
1402     *out = ebur128_energy_to_loudness(h_en) - ebur128_energy_to_loudness(l_en);
1403   } else {
1404     free(stl_vector);
1405     *out = 0.0;
1406   }
1407 
1408   return EBUR128_SUCCESS;
1409 }
1410 
ebur128_loudness_range(ebur128_state * st,double * out)1411 int ebur128_loudness_range(ebur128_state* st, double* out) {
1412   return ebur128_loudness_range_multiple(&st, 1, out);
1413 }
1414 
ebur128_sample_peak(ebur128_state * st,unsigned int channel_number,double * out)1415 int ebur128_sample_peak(ebur128_state* st,
1416                         unsigned int channel_number,
1417                         double* out) {
1418   if ((st->mode & EBUR128_MODE_SAMPLE_PEAK) != EBUR128_MODE_SAMPLE_PEAK) {
1419     return EBUR128_ERROR_INVALID_MODE;
1420   }
1421 
1422   if (channel_number >= st->channels) {
1423     return EBUR128_ERROR_INVALID_CHANNEL_INDEX;
1424   }
1425 
1426   *out = st->d->sample_peak[channel_number];
1427   return EBUR128_SUCCESS;
1428 }
1429 
ebur128_prev_sample_peak(ebur128_state * st,unsigned int channel_number,double * out)1430 int ebur128_prev_sample_peak(ebur128_state* st,
1431                              unsigned int channel_number,
1432                              double* out) {
1433   if ((st->mode & EBUR128_MODE_SAMPLE_PEAK) != EBUR128_MODE_SAMPLE_PEAK) {
1434     return EBUR128_ERROR_INVALID_MODE;
1435   }
1436 
1437   if (channel_number >= st->channels) {
1438     return EBUR128_ERROR_INVALID_CHANNEL_INDEX;
1439   }
1440 
1441   *out = st->d->prev_sample_peak[channel_number];
1442   return EBUR128_SUCCESS;
1443 }
1444 
ebur128_true_peak(ebur128_state * st,unsigned int channel_number,double * out)1445 int ebur128_true_peak(ebur128_state* st,
1446                       unsigned int channel_number,
1447                       double* out) {
1448   if ((st->mode & EBUR128_MODE_TRUE_PEAK) != EBUR128_MODE_TRUE_PEAK) {
1449     return EBUR128_ERROR_INVALID_MODE;
1450   }
1451 
1452   if (channel_number >= st->channels) {
1453     return EBUR128_ERROR_INVALID_CHANNEL_INDEX;
1454   }
1455 
1456   *out = EBUR128_MAX(st->d->true_peak[channel_number],
1457                      st->d->sample_peak[channel_number]);
1458   return EBUR128_SUCCESS;
1459 }
1460 
ebur128_prev_true_peak(ebur128_state * st,unsigned int channel_number,double * out)1461 int ebur128_prev_true_peak(ebur128_state* st,
1462                            unsigned int channel_number,
1463                            double* out) {
1464   if ((st->mode & EBUR128_MODE_TRUE_PEAK) != EBUR128_MODE_TRUE_PEAK) {
1465     return EBUR128_ERROR_INVALID_MODE;
1466   }
1467 
1468   if (channel_number >= st->channels) {
1469     return EBUR128_ERROR_INVALID_CHANNEL_INDEX;
1470   }
1471 
1472   *out = EBUR128_MAX(st->d->prev_true_peak[channel_number],
1473                      st->d->prev_sample_peak[channel_number]);
1474   return EBUR128_SUCCESS;
1475 }
1476