1 // ----------------------------------------------------------------------------
2 //
3 //  Copyright (C) 2010-2018 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 3 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, see <http://www.gnu.org/licenses/>.
17 //
18 // ----------------------------------------------------------------------------
19 
20 
21 #include <assert.h>
22 #include <string.h>
23 #include <math.h>
24 #include "peaklim.h"
25 
26 
init(int hlen)27 void Histmin::init (int hlen)
28 {
29     int i;
30 
31     assert (hlen <= SIZE);
32     _hlen = hlen;
33     _hold = hlen;
34     _wind = 0;
35     _vmin = 1;
36     for (i = 0; i < SIZE; i++) _hist [i] = _vmin;
37 }
38 
39 
write(float v)40 float Histmin::write (float v)
41 {
42     int i, j;
43 
44     i = _wind;
45     _hist [i] = v;
46     if (v <= _vmin)
47     {
48 	_vmin = v;
49 	_hold = _hlen;
50     }
51     else if (--_hold == 0)
52     {
53 	_vmin = v;
54 	_hold = _hlen;
55 	for (j = 1 - _hlen; j < 0; j++)
56 	{
57 	    v = _hist [(i + j) & MASK];
58 	    if (v < _vmin)
59 	    {
60 		_vmin = v;
61 		_hold = _hlen + j;
62 	    }
63 	}
64     }
65     _wind = ++i & MASK;
66     return _vmin;
67 }
68 
69 
70 
Peaklim(void)71 Peaklim::Peaklim (void) :
72     _fsamp (0),
73     _nchan (0),
74     _rstat (false),
75     _peak (0),
76     _gmax (1),
77     _gmin (1)
78 {
79     for (int i = 0; i < MAXCHAN; i++) _dbuff [i] = 0;
80 }
81 
82 
~Peaklim(void)83 Peaklim::~Peaklim (void)
84 {
85     fini ();
86 }
87 
88 
set_inpgain(float v)89 void Peaklim::set_inpgain (float v)
90 {
91     if (v >  30.0f) v =  30.0f;
92     if (v < -30.0f) v = -30.0f;
93     _g1 = powf (10.0f, 0.05f * v);
94 }
95 
96 
set_threshold(float v)97 void Peaklim::set_threshold (float v)
98 {
99     if (v >   0.0f) v =   0.0f;
100     if (v < -30.0f) v = -30.0f;
101     _gt = powf (10.0f, -0.05f * v);
102 }
103 
104 
set_release(float v)105 void Peaklim::set_release (float v)
106 {
107     if (v > 1.0f)  v = 1.0f;
108     if (v < 1e-3f) v = 1e-3f;
109     _w3 = 1.0f / (v * _fsamp);
110 }
111 
112 
init(float fsamp,int nchan)113 void Peaklim::init (float fsamp, int nchan)
114 {
115     int i, k1, k2;
116 
117     fini ();
118     if (nchan > MAXCHAN) nchan = MAXCHAN;
119     _fsamp = fsamp;
120     _nchan = nchan;
121     if      (fsamp > 130000) _div1 = 32;
122     else if (fsamp >  65000) _div1 = 16;
123     else 	             _div1 = 8;
124     _div2 = 8;
125     k1 = (int)(ceilf (1.2e-3f * _fsamp / _div1));
126     k2 = 12;
127     _delay = k1 * _div1;
128     for (_dsize = 64; _dsize < _delay + _div1; _dsize *= 2);
129     _dmask = _dsize - 1;
130     _delri = 0;
131     for (i = 0; i < _nchan; i++)
132     {
133          _dbuff [i] = new float [_dsize];
134 	 memset (_dbuff [i], 0, _dsize * sizeof (float));
135     }
136     _hist1.init (k1 + 1);
137     _hist2.init (k2);
138     _c1 = _div1;
139     _c2 = _div2;
140     _m1 = 0.0f;
141     _m2 = 0.0f;
142     _wlf = 6.28f * 500.0f / fsamp;
143     _w1 = 10.0f / _delay;
144     _w2 = _w1 / _div2;
145     _w3 = 1.0f / (0.01f * fsamp);
146     for (i = 0; i < _nchan; i++) _zlf [i] = 0.0f;
147     _z1 = 1.0f;
148     _z2 = 1.0f;
149     _z3 = 1.0f;
150     _gt = 1.0f;
151     _g0 = 1.0f;
152     _g1 = 1.0f;
153     _dg = 0.0f;
154     _gmax = 1.0f;
155     _gmin = 1.0f;
156 }
157 
158 
fini(void)159 void Peaklim::fini (void)
160 {
161     int i;
162 
163     for (i = 0; i < MAXCHAN; i++)
164     {
165 	delete[] _dbuff [i];
166 	_dbuff [i] = 0;
167     }
168     _nchan = 0;
169 }
170 
process(int nframes,float * inp[],float * out[])171 void Peaklim::process (int nframes, float *inp[], float *out[])
172 {
173     int     i, j, k, n, ri, wi;
174     float   g, d, h1, h2, m1, m2, x, z, z1, z2, z3, pk, t0, t1, *p;
175 
176     ri = _delri;
177     wi = (ri + _delay) & _dmask;
178     h1 = _hist1.vmin ();
179     h2 = _hist2.vmin ();
180     m1 = _m1;
181     m2 = _m2;
182     z1 = _z1;
183     z2 = _z2;
184     z3 = _z3;
185 
186     if (_rstat)
187     {
188 	_rstat = false;
189 	pk = 0;
190 	t0 = _gmax;
191 	t1 = _gmin;
192     }
193     else
194     {
195 	pk = _peak;
196 	t0 = _gmin;
197 	t1 = _gmax;
198     }
199 
200 
201     k = 0;
202     while (nframes)
203     {
204 	n = (_c1 < nframes) ? _c1 : nframes;
205 
206         g = _g0;
207 	for (j = 0; j < _nchan; j++)
208 	{
209 	    p = inp [j] + k;
210 	    z = _zlf [j];
211 	    g = _g0;
212 	    d = _dg;
213 	    for (i = 0; i < n; i++)
214 	    {
215 	        x = g * *p++;
216 		g += d;
217 	        _dbuff [j][wi + i] = x;
218 		z += _wlf * (x - z) + 1e-20f;
219                 x = fabsf (x);
220 	        if (x > m1) m1 = x;
221                 x = fabsf (z);
222 	        if (x > m2) m2 = x;
223 	    }
224 	    _zlf [j] = z;
225 	}
226 	_g0 = g;
227 
228 	_c1 -= n;
229 	if (_c1 == 0)
230 	{
231 	    m1 *= _gt;
232 	    if (m1 > pk) pk = m1;
233   	    h1 = (m1 > 1.0f) ? 1.0f / m1 : 1.0f;
234 	    h1 = _hist1.write (h1);
235 	    m1 = 0;
236 	    _c1 = _div1;
237 	    if (--_c2 == 0)
238 	    {
239 		m2 *= _gt;
240 	        h2 = (m2 > 1.0f) ? 1.0f / m2 : 1.0f;
241 	        h2 = _hist2.write (h2);
242 	        m2 = 0;
243 	        _c2 = _div2;
244 		_dg = _g1 - _g0;
245 		if (fabsf (_dg) < 1e-9f)
246 		{
247 		    _g0 = _g1;
248 		    _dg = 0;
249 		}
250 		else
251 		{
252 		    _dg /= _div1 * _div2;
253 		}
254 	    }
255 	}
256 
257 	for (i = 0; i < n; i++)
258 	{
259 	    z1 += _w1 * (h1 - z1);
260 	    z2 += _w2 * (h2 - z2);
261 	    z = (z2 < z1) ? z2 : z1;
262 	    if (z < z3)  z3 += _w1 * (z - z3);
263 	    else         z3 += _w3 * (z - z3);
264 	    if (z3 > t1) t1 = z3;
265 	    if (z3 < t0) t0 = z3;
266 	    for (j = 0; j < _nchan; j++)
267 	    {
268                 out [j][k + i] = z3 * _dbuff [j][ri + i];
269 	    }
270 	}
271 
272 	wi = (wi + n) & _dmask;
273 	ri = (ri + n) & _dmask;
274 	k += n;
275 	nframes -= n;
276     }
277 
278     _delri = ri;
279     _m1 = m1;
280     _m2 = m2;
281     _z1 = z1;
282     _z2 = z2;
283     _z3 = z3;
284     _peak = pk;
285     _gmin = t0;
286     _gmax = t1;
287 }
288