1 // ----------------------------------------------------------------------------
2 //
3 // Copyright (C) 2008-2017 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 <string.h>
22 #include <math.h>
23 #include <assert.h>
24 #include "dplimit1.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
Dplimit1(void)71 Dplimit1::Dplimit1 (void) :
72 _state (BYPASS),
73 _fsamp (0),
74 _nchan (0),
75 _rstat (false),
76 _peak (0),
77 _gmax (1),
78 _gmin (1)
79 {
80 for (int i = 0; i < MAXCHAN; i++) _dbuff [i] = 0;
81 }
82
83
~Dplimit1(void)84 Dplimit1::~Dplimit1 (void)
85 {
86 fini ();
87 }
88
89
set_threshd(float v)90 void Dplimit1::set_threshd (float v)
91 {
92 _gt = powf (10.0f, -0.05f * v);
93 }
94
95
set_reltime(float v)96 void Dplimit1::set_reltime (float v)
97 {
98 _w3 = 1.0f / (v * _fsamp);
99 }
100
101
init(float fsamp,int nchan,float threshd,float reltime)102 void Dplimit1::init (float fsamp, int nchan, float threshd, float reltime)
103 {
104 int i;
105
106 fini ();
107 if (nchan > MAXCHAN) nchan = MAXCHAN;
108 _fsamp = fsamp;
109 _nchan = nchan;
110 if (fsamp > 130000) _div1 = 32;
111 else if (fsamp > 65000) _div1 = 16;
112 else _div1 = 8;
113 _div2 = 8;
114 _len1 = (int)(ceilf (1.2e-3f * _fsamp / _div1));
115 _len2 = 12;
116 _delay = _len1 * _div1;
117 for (_dsize = 64; _dsize < _delay + _div1; _dsize *= 2);
118 _dmask = _dsize - 1;
119 _delri = 0;
120 for (i = 0; i < _nchan; i++)
121 {
122 _dbuff [i] = new float [_dsize];
123 memset (_dbuff [i], 0, _dsize * sizeof (float));
124 }
125 set_threshd (threshd);
126 set_reltime (reltime);
127 _hist1.init (_len1 + 1);
128 _hist2.init (_len2);
129 _c1 = _div1;
130 _c2 = _div2;
131 _m1 = 0.0f;
132 _m2 = 0.0f;
133 _wlf = 6.28f * 500.0f / fsamp;
134 _w1 = 10.0f / _delay;
135 _w2 = _w1 / _div2;
136 for (i = 0; i < _nchan; i++) _zlf [i] = 0.0f;
137 _z1 = 1.0f;
138 _z2 = 1.0f;
139 _z3 = 1.0f;
140 _gmax = 1.0f;
141 _gmin = 1.0f;
142 _state = ACTIVE;
143 }
144
145
fini(void)146 void Dplimit1::fini (void)
147 {
148 int i;
149
150 for (i = 0; i < MAXCHAN; i++)
151 {
152 delete[] _dbuff [i];
153 _dbuff [i] = 0;
154 }
155 _nchan = 0;
156 _state = BYPASS;
157 }
158
159
prepare(int nsamp)160 void Dplimit1::prepare (int nsamp)
161 {
162 }
163
164
process1(int nsamp,float * data[])165 void Dplimit1::process1 (int nsamp, float *data[])
166 {
167 int i, j, k, n, c1, c2, ri, wi;
168 float g1, g2, m1, m2, x, z, z1, z2, z3, pk, t0, t1, *p;
169 float *dp [MAXCHAN];
170
171 ri = _delri;
172 wi = (ri + _delay) & _dmask;
173 g1 = _hist1.vmin ();
174 c1 = _c1;
175 m1 = _m1;
176 g2 = _hist2.vmin ();
177 c2 = _c2;
178 m2 = _m2;
179 z1 = _z1;
180 z2 = _z2;
181 z3 = _z3;
182 n = _nchan;
183 for (j = 0; j < n; j++) dp [j] = data [j];
184
185 if (_rstat)
186 {
187 _rstat = false;
188 pk = 0;
189 t0 = _gmax;
190 t1 = _gmin;
191 }
192 else
193 {
194 pk = _peak;
195 t0 = _gmin;
196 t1 = _gmax;
197 }
198
199 while (nsamp)
200 {
201 k = (c1 < nsamp) ? c1 : nsamp;
202
203 for (j = 0; j < n; j++)
204 {
205 p = dp [j];
206 z = _zlf [j];
207 for (i = 0; i < k; i++)
208 {
209 x = *p++;
210 z += _wlf * (x - z) + 1e-20f;
211 _dbuff [j][wi + i] = x;
212 x = fabsf (x);
213 if (x > m1) m1 = x;
214 x = fabsf (z);
215 if (x > m2) m2 = x;
216 }
217 _zlf [j] = z;
218 }
219
220 c1 -= k;
221 if (c1 == 0)
222 {
223 m1 *= _gt;
224 if (m1 > pk) pk = m1;
225 g1 = (m1 > 1.0f) ? 1.0f / m1 : 1.0f;
226 g1 = _hist1.write (g1);
227 c1 = _div1;
228 m1 = 0;
229 if (--c2 == 0)
230 {
231 m2 *= _gt;
232 g2 = (m2 > 1.0f) ? 1.0f / m2 : 1.0f;
233 g2 = _hist2.write (g2);
234 c2 = _div2;
235 m2 = 0;
236 }
237 }
238
239 for (i = 0; i < k; i++)
240 {
241 z1 += _w1 * (g1 - z1);
242 z2 += _w2 * (g2 - z2);
243 z = (z2 < z1) ? z2 : z1;
244 if (z < z3) z3 += _w1 * (z - z3);
245 else z3 += _w3 * (z - z3);
246 if (z3 > t1) t1 = z3;
247 if (z3 < t0) t0 = z3;
248 for (j = 0; j < n; j++)
249 {
250 dp [j][i] = z3 * _dbuff [j][ri + i];
251 }
252 }
253
254 for (j = 0; j < n; j++) dp [j] += k;
255 wi = (wi + k) & _dmask;
256 ri = (ri + k) & _dmask;
257 nsamp -= k;
258 }
259
260 _delri = ri;
261 _c1 = c1;
262 _m1 = m1;
263 _c2 = c2;
264 _m2 = m2;
265 _z1 = z1;
266 _z2 = z2;
267 _z3 = z3;
268 _peak = pk;
269 _gmin = t0;
270 _gmax = t1;
271 }
272