]> git.sesse.net Git - nageru/blob - ebu_r128_proc.cc
Support searching for the theme in multiple directories.
[nageru] / 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