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 #include "ebu_r128_proc.h"
23 #include <math.h>
24 #include <string.h>
25
26 #ifdef COMPILER_MSVC
27 #include <float.h>
28 // C99 'isfinite()' is not available in MSVC.
29 #define isfinite_local(val) (bool)_finite((double)val)
30 #else
31 #define isfinite_local isfinite
32 #endif
33
34 using namespace FonsEBU;
35
36 float Ebu_r128_proc::Ebu_r128_hist::_bin_power[100] = { 0.0f };
37 float Ebu_r128_proc::_chan_gain[5] = { 1.0f, 1.0f, 1.0f, 1.41f, 1.41f };
38
Ebu_r128_hist()39 Ebu_r128_proc::Ebu_r128_hist::Ebu_r128_hist ()
40 {
41 _histc = new int[751];
42 initstat ();
43 reset ();
44 }
45
~Ebu_r128_hist()46 Ebu_r128_proc::Ebu_r128_hist::~Ebu_r128_hist ()
47 {
48 delete[] _histc;
49 }
50
51 void
reset()52 Ebu_r128_proc::Ebu_r128_hist::reset ()
53 {
54 memset (_histc, 0, 751 * sizeof (float));
55 _count = 0;
56 _error = 0;
57 }
58
59 void
initstat()60 Ebu_r128_proc::Ebu_r128_hist::initstat ()
61 {
62 if (_bin_power[0])
63 return;
64 for (int i = 0; i < 100; i++) {
65 _bin_power[i] = powf (10.0f, i / 100.0f);
66 }
67 }
68
69 void
addpoint(float v)70 Ebu_r128_proc::Ebu_r128_hist::addpoint (float v)
71 {
72 int k = (int)floorf (10 * v + 700.5f);
73 if (k < 0)
74 return;
75 if (k > 750) {
76 k = 750;
77 _error++;
78 }
79 _histc[k]++;
80 _count++;
81 }
82
83 float
integrate(int i)84 Ebu_r128_proc::Ebu_r128_hist::integrate (int i)
85 {
86 int j, k, n;
87 float s;
88
89 j = i % 100;
90 n = 0;
91 s = 0;
92 while (i <= 750) {
93 k = _histc[i++];
94 n += k;
95 s += k * _bin_power[j++];
96 if (j == 100) {
97 j = 0;
98 s /= 10.0f;
99 }
100 }
101 return s / n;
102 }
103
104 void
calc_integ(float * vi,float * th)105 Ebu_r128_proc::Ebu_r128_hist::calc_integ (float* vi, float* th)
106 {
107 int k;
108 float s;
109
110 if (_count < 50) {
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) {
120 *th = 10 * log10f (s) - 10.0f;
121 }
122 k = (int)(floorf (100 * log10f (s) + 0.5f)) + 600;
123 if (k < 0) {
124 k = 0;
125 }
126 s = integrate (k);
127 *vi = 10 * log10f (s);
128 }
129
130 void
calc_range(float * v0,float * v1,float * th)131 Ebu_r128_proc::Ebu_r128_hist::calc_range (float* v0, float* v1, float* th)
132 {
133 int i, j, k, n;
134 float a, b, s;
135
136 if (_count < 20) {
137 *v0 = -200.0f;
138 *v1 = -200.0f;
139 return;
140 }
141 s = integrate (0);
142 if (th) {
143 *th = 10 * log10f (s) - 20.0f;
144 }
145 k = (int)(floorf (100 * log10f (s) + 0.5)) + 500;
146 if (k < 0) {
147 k = 0;
148 }
149 for (i = k, n = 0; i <= 750; i++) {
150 n += _histc[i];
151 }
152 a = 0.10f * n;
153 b = 0.95f * n;
154 for (i = k, s = 0; s < a; i++) {
155 s += _histc[i];
156 }
157 for (j = 750, s = n; s > b; j--) {
158 s -= _histc[j];
159 }
160 *v0 = (i - 701) / 10.0f;
161 *v1 = (j - 699) / 10.0f;
162 }
163
Ebu_r128_proc()164 Ebu_r128_proc::Ebu_r128_proc ()
165 {
166 reset ();
167 }
168
~Ebu_r128_proc()169 Ebu_r128_proc::~Ebu_r128_proc ()
170 {
171 }
172
173 void
init(int nchan,float fsamp)174 Ebu_r128_proc::init (int nchan, float fsamp)
175 {
176 _nchan = nchan;
177 _fsamp = fsamp;
178 _fragm = (int)fsamp / 20;
179 detect_init (_fsamp);
180 reset ();
181 }
182
183 void
reset()184 Ebu_r128_proc::reset ()
185 {
186 _integr = false;
187 _frcnt = _fragm;
188 _frpwr = 1e-30f;
189 _wrind = 0;
190 _div1 = 0;
191 _div2 = 0;
192 _loudness_M = -200.0f;
193 _loudness_S = -200.0f;
194 memset (_power, 0, 64 * sizeof (float));
195 integr_reset ();
196 detect_reset ();
197 }
198
199 void
integr_reset()200 Ebu_r128_proc::integr_reset ()
201 {
202 _hist_M.reset ();
203 _hist_S.reset ();
204 _maxloudn_M = -200.0f;
205 _maxloudn_S = -200.0f;
206 _integrated = -200.0f;
207 _integ_thr = -200.0f;
208 _range_min = -200.0f;
209 _range_max = -200.0f;
210 _range_thr = -200.0f;
211 _div1 = _div2 = 0;
212 }
213
214 void
process(int nfram,const float * const * input)215 Ebu_r128_proc::process (int nfram, const float* const* input)
216 {
217 int i, k;
218
219 for (i = 0; i < _nchan; i++) {
220 _ipp[i] = input[i];
221 }
222 while (nfram) {
223 k = (_frcnt < nfram) ? _frcnt : nfram;
224 _frpwr += detect_process (k);
225 _frcnt -= k;
226 if (_frcnt == 0) {
227 _power[_wrind++] = _frpwr / _fragm;
228 _frcnt = _fragm;
229 _frpwr = 1e-30f;
230 _wrind &= 63;
231 _loudness_M = addfrags (8);
232 _loudness_S = addfrags (60);
233
234 if (!isfinite_local (_loudness_M) || _loudness_M < -200.f) {
235 _loudness_M = -200.0f;
236 }
237 if (!isfinite_local (_loudness_S) || _loudness_S < -200.f) {
238 _loudness_S = -200.0f;
239 }
240 if (_loudness_M > _maxloudn_M) {
241 _maxloudn_M = _loudness_M;
242 }
243 if (_loudness_S > _maxloudn_S) {
244 _maxloudn_S = _loudness_S;
245 }
246
247 if (_integr) {
248 if (++_div1 == 2) {
249 _hist_M.addpoint (_loudness_M);
250 _div1 = 0;
251 }
252 if (++_div2 == 10) {
253 _hist_S.addpoint (_loudness_S);
254 _div2 = 0;
255 _hist_M.calc_integ (&_integrated, &_integ_thr);
256 _hist_S.calc_range (&_range_min, &_range_max, &_range_thr);
257 }
258 }
259 }
260
261 for (i = 0; i < _nchan; i++) {
262 _ipp[i] += k;
263 }
264 nfram -= k;
265 }
266 }
267
268 float
addfrags(int nfrag)269 Ebu_r128_proc::addfrags (int nfrag)
270 {
271 int i, k;
272 float s;
273
274 s = 0;
275 k = (_wrind - nfrag) & 63;
276 for (i = 0; i < nfrag; i++) {
277 s += _power[(i + k) & 63];
278 }
279 return -0.6976f + 10 * log10f (s / nfrag);
280 }
281
282 void
detect_init(float fsamp)283 Ebu_r128_proc::detect_init (float fsamp)
284 {
285 float a, b, c, d, r, u1, u2, w1, w2;
286
287 r = 1 / tan (4712.3890f / fsamp);
288 w1 = r / 1.12201f;
289 w2 = r * 1.12201f;
290 u1 = u2 = 1.4085f + 210.0f / fsamp;
291
292 a = u1 * w1;
293 b = w1 * w1;
294 c = u2 * w2;
295 d = w2 * w2;
296
297 r = 1 + a + b;
298 _a0 = (1 + c + d) / r;
299 _a1 = (2 - 2 * d) / r;
300 _a2 = (1 - c + d) / r;
301 _b1 = (2 - 2 * b) / r;
302 _b2 = (1 - a + b) / r;
303
304 r = 48.0f / fsamp;
305 a = 4.9886075f * r;
306 b = 6.2298014f * r * r;
307 r = 1 + a + b;
308 a *= 2 / r;
309 b *= 4 / r;
310 _c3 = a + b;
311 _c4 = b;
312
313 r = 1.004995f / r;
314 _a0 *= r;
315 _a1 *= r;
316 _a2 *= r;
317 }
318
319 void
detect_reset()320 Ebu_r128_proc::detect_reset ()
321 {
322 for (int i = 0; i < MAXCH; i++)
323 _fst[i].reset ();
324 }
325
326 float
detect_process(int nfram)327 Ebu_r128_proc::detect_process (int nfram)
328 {
329 int i, j;
330 float si, sj;
331 float x, y, z1, z2, z3, z4;
332 float const* p;
333 Ebu_r128_fst* S;
334
335 si = 0;
336 for (i = 0, S = _fst; i < _nchan; i++, S++) {
337 z1 = S->_z1;
338 z2 = S->_z2;
339 z3 = S->_z3;
340 z4 = S->_z4;
341 p = _ipp[i];
342 sj = 0;
343 for (j = 0; j < nfram; j++) {
344 x = p[j] - _b1 * z1 - _b2 * z2 + 1e-15f;
345 y = _a0 * x + _a1 * z1 + _a2 * z2 - _c3 * z3 - _c4 * z4;
346 z2 = z1;
347 z1 = x;
348 z4 += z3;
349 z3 += y;
350 sj += y * y;
351 }
352 if (_nchan == 1) {
353 si = 2 * sj;
354 } else {
355 si += _chan_gain[i] * sj;
356 }
357 S->_z1 = !isfinite_local (z1) ? 0 : z1;
358 S->_z2 = !isfinite_local (z2) ? 0 : z2;
359 S->_z3 = !isfinite_local (z3) ? 0 : z3;
360 S->_z4 = !isfinite_local (z4) ? 0 : z4;
361 }
362 return si;
363 }
364