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