libs/vamp-plugins/ebu_r128_proc.cc

   1 // ------------------------------------------------------------------------
   2 //
   3 //  Copyright (C) 2010-2011 Fons Adriaensen <fons@linuxaudio.org>
   4 //  Copyright (C) 2015 Robin Gareus <robin@gareus.org>
   5 //
   6 //  This program is free software; you can redistribute it and/or modify
   7 //  it under the terms of the GNU General Public License as published by
   8 //  the Free Software Foundation; either version 2 of the License, or
   9 //  (at your option) any later version.
  10 //
  11 //  This program is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 //  GNU General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU General Public License
  17 //  along with this program; if not, write to the Free Software
  18 //  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  19 //
  20 // ------------------------------------------------------------------------
  21
  22
  23 #include <string.h>
  24 #include <math.h>
  25 #include "ebu_r128_proc.h"
  26
  27 #ifdef COMPILER_MSVC
  28 #include <float.h>
  29 // C99 'isfinite()' is not available in MSVC.
  30 #define isfinite_local(val) (bool)_finite((double)val)
  31 #else
  32 #define isfinite_local isfinite
  33 #endif
  34
  35 namespace Fons {
  36
  37 float Ebu_r128_hist::_bin_power [100] = { 0.0f };
  38 float Ebu_r128_proc::_chan_gain [5] = { 1.0f, 1.0f, 1.0f, 1.41f, 1.41f };
  39
  40
  41 Ebu_r128_hist::Ebu_r128_hist (void)
  42 {
  43     _histc = new int [751];
  44     initstat ();
  45     reset ();
  46 }
  47
  48
  49 Ebu_r128_hist::~Ebu_r128_hist (void)
  50 {
  51     delete[] _histc;
  52 }
  53
  54
  55 void Ebu_r128_hist::reset (void)
  56 {
  57     memset (_histc, 0, 751 * sizeof (float));
  58     _count = 0;
  59     _error = 0;
  60 }
  61
  62
  63 void Ebu_r128_hist::initstat (void)
  64 {
  65     int i;
  66
  67     if (_bin_power [0]) return;
  68     for (i = 0; i < 100; i++)
  69     {
  70         _bin_power [i] = powf (10.0f, i / 100.0f);
  71     }
  72 }
  73
  74
  75 void Ebu_r128_hist::addpoint (float v)
  76 {
  77     int k;
  78
  79     k = (int) floorf (10 * v + 700.5f);
  80     if (k < 0) return;
  81     if (k > 750)
  82     {
  83         k = 750;
  84         _error++;
  85     }
  86     _histc [k]++;
  87     _count++;
  88 }
  89
  90
  91 float Ebu_r128_hist::integrate (int i)
  92 {
  93     int   j, k, n;
  94     float s;
  95
  96     j = i % 100;
  97     n = 0;
  98     s = 0;
  99     while (i <= 750)
 100     {
 101         k = _histc [i++];
 102         n += k;
 103         s += k * _bin_power [j++];
 104         if (j == 100)
 105         {
 106             j = 0;
 107             s /= 10.0f;
 108         }
 109     }
 110     return s / n;
 111 }
 112
 113
 114 void Ebu_r128_hist::calc_integ (float *vi, float *th)
 115 {
 116     int   k;
 117     float s;
 118
 119     if (_count < 50)
 120     {
 121         *vi = -200.0f;
 122         return;
 123     }
 124     s = integrate (0);
 125 //  Original threshold was -8 dB below result of first integration
 126 //    if (th) *th = 10 * log10f (s) - 8.0f;
 127 //    k = (int)(floorf (100 * log10f (s) + 0.5f)) + 620;
 128 //  Threshold redefined to -10 dB below result of first integration
 129     if (th) *th = 10 * log10f (s) - 10.0f;
 130     k = (int)(floorf (100 * log10f (s) + 0.5f)) + 600;
 131     if (k < 0) k = 0;
 132     s = integrate (k);
 133     *vi = 10 * log10f (s);
 134 }
 135
 136
 137 void Ebu_r128_hist::calc_range (float *v0, float *v1, float *th)
 138 {
 139     int   i, j, k, n;
 140     float a, b, s;
 141
 142     if (_count < 20)
 143     {
 144         *v0 = -200.0f;
 145         *v1 = -200.0f;
 146         return;
 147     }
 148     s = integrate (0);
 149     if (th) *th = 10 * log10f (s) - 20.0f;
 150     k = (int)(floorf (100 * log10f (s) + 0.5)) + 500;
 151     if (k < 0) k = 0;
 152     for (i = k, n = 0; i <= 750; i++) n += _histc [i];
 153     a = 0.10f * n;
 154     b = 0.95f * n;
 155     for (i =   k, s = 0; s < a; i++) s += _histc [i];
 156     for (j = 750, s = n; s > b; j--) s -= _histc [j];
 157     *v0 = (i - 701) / 10.0f;
 158     *v1 = (j - 699) / 10.0f;
 159 }
 160
 161
 162
 163
 164 Ebu_r128_proc::Ebu_r128_proc (void)
 165 {
 166     reset ();
 167 }
 168
 169
 170 Ebu_r128_proc::~Ebu_r128_proc (void)
 171 {
 172 }
 173
 174
 175 void Ebu_r128_proc::init (int nchan, float fsamp)
 176 {
 177     _nchan = nchan;
 178     _fsamp = fsamp;
 179     _fragm = (int) fsamp / 20;
 180     detect_init (_fsamp);
 181     reset ();
 182 }
 183
 184
 185 void Ebu_r128_proc::reset (void)
 186 {
 187     _integr = false;
 188     _frcnt = _fragm;
 189     _frpwr = 1e-30f;
 190     _wrind  = 0;
 191     _div1 = 0;
 192     _div2 = 0;
 193     _loudness_M = -200.0f;
 194     _loudness_S = -200.0f;
 195     memset (_power, 0, 64 * sizeof (float));
 196     integr_reset ();
 197     detect_reset ();
 198 }
 199
 200
 201 void Ebu_r128_proc::integr_reset (void)
 202 {
 203     _hist_M.reset ();
 204     _hist_S.reset ();
 205     _maxloudn_M = -200.0f;
 206     _maxloudn_S = -200.0f;
 207     _integrated = -200.0f;
 208     _integ_thr  = -200.0f;
 209     _range_min  = -200.0f;
 210     _range_max  = -200.0f;
 211     _range_thr  = -200.0f;
 212     _div1 = _div2 = 0;
 213 }
 214
 215
 216 void Ebu_r128_proc::process (int nfram, const float *const *input)
 217 {
 218     int  i, k;
 219
 220     for (i = 0; i < _nchan; i++) _ipp [i] = input [i];
 221     while (nfram)
 222     {
 223         k = (_frcnt < nfram) ? _frcnt : nfram;
 224         _frpwr += detect_process (k);
 225         _frcnt -= k;
 226         if (_frcnt == 0)
 227         {
 228             _power [_wrind++] = _frpwr / _fragm;
 229             _frcnt = _fragm;
 230             _frpwr = 1e-30f;
 231             _wrind &= 63;
 232             _loudness_M = addfrags (8);
 233             _loudness_S = addfrags (60);
 234             if (!isfinite_local(_loudness_M) || _loudness_M < -200.f) _loudness_M = -200.0f;
 235             if (!isfinite_local(_loudness_S) || _loudness_S < -200.f) _loudness_S = -200.0f;
 236             if (_loudness_M > _maxloudn_M) _maxloudn_M = _loudness_M;
 237             if (_loudness_S > _maxloudn_S) _maxloudn_S = _loudness_S;
 238             if (_integr)
 239             {
 240                 if (++_div1 == 2)
 241                 {
 242                     _hist_M.addpoint (_loudness_M);
 243                     _div1 = 0;
 244                 }
 245                 if (++_div2 == 10)
 246                 {
 247                     _hist_S.addpoint (_loudness_S);
 248                     _div2 = 0;
 249                     _hist_M.calc_integ (&_integrated, &_integ_thr);
 250                     _hist_S.calc_range (&_range_min, &_range_max, &_range_thr);
 251                 }
 252             }
 253         }
 254         for (i = 0; i < _nchan; i++) _ipp [i] += k;
 255         nfram -= k;
 256     }
 257 }
 258
 259
 260 float Ebu_r128_proc::addfrags (int nfrag)
 261 {
 262     int    i, k;
 263     float  s;
 264
 265     s = 0;
 266     k = (_wrind - nfrag) & 63;
 267     for (i = 0; i < nfrag; i++) s += _power [(i + k) & 63];
 268     return -0.6976f + 10 * log10f (s / nfrag);
 269 }
 270
 271
 272 void Ebu_r128_proc::detect_init (float fsamp)
 273 {
 274     float a, b, c, d, r, u1, u2, w1, w2;
 275
 276     r = 1 / tan (4712.3890f / fsamp);
 277     w1 = r / 1.12201f;
 278     w2 = r * 1.12201f;
 279     u1 = u2 = 1.4085f + 210.0f / fsamp;
 280     a = u1 * w1;
 281     b = w1 * w1;
 282     c = u2 * w2;
 283     d = w2 * w2;
 284     r = 1 + a + b;
 285     _a0 = (1 + c + d) / r;
 286     _a1 = (2 - 2 * d) / r;
 287     _a2 = (1 - c + d) / r;
 288     _b1 = (2 - 2 * b) / r;
 289     _b2 = (1 - a + b) / r;
 290     r = 48.0f / fsamp;
 291     a = 4.9886075f * r;
 292     b = 6.2298014f * r * r;
 293     r = 1 + a + b;
 294     a *= 2 / r;
 295     b *= 4 / r;
 296     _c3 = a + b;
 297     _c4 = b;
 298     r = 1.004995f / r;
 299     _a0 *= r;
 300     _a1 *= r;
 301     _a2 *= r;
 302 }
 303
 304
 305 void Ebu_r128_proc::detect_reset (void)
 306 {
 307     for (int i = 0; i < MAXCH; i++) _fst [i].reset ();
 308 }
 309
 310
 311 float Ebu_r128_proc::detect_process (int nfram)
 312 {
 313     int   i, j;
 314     float si, sj;
 315     float x, y, z1, z2, z3, z4;
 316     float const *p;
 317     Ebu_r128_fst *S;
 318
 319     si = 0;
 320     for (i = 0, S = _fst; i < _nchan; i++, S++)
 321     {
 322         z1 = S->_z1;
 323         z2 = S->_z2;
 324         z3 = S->_z3;
 325         z4 = S->_z4;
 326         p = _ipp [i];
 327         sj = 0;
 328         for (j = 0; j < nfram; j++)
 329         {
 330             x = p [j] - _b1 * z1 - _b2 * z2 + 1e-15f;
 331             y = _a0 * x + _a1 * z1 + _a2 * z2 - _c3 * z3 - _c4 * z4;
 332             z2 = z1;
 333             z1 = x;
 334             z4 += z3;
 335             z3 += y;
 336             sj += y * y;
 337         }
 338         if (_nchan == 1) si = 2 * sj;
 339         else si += _chan_gain [i] * sj;
 340         S->_z1 = !isfinite_local(z1) ? 0 : z1;
 341         S->_z2 = !isfinite_local(z2) ? 0 : z2;
 342         S->_z3 = !isfinite_local(z3) ? 0 : z3;
 343         S->_z4 = !isfinite_local(z4) ? 0 : z4;
 344     }
 345     return si;
 346 }
 347
 348 };