git.sesse.net Git - nageru/blob - ebu_r128_proc.cc

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