3e2ba74004cc84a24d7be355fdb6da9d7fadb9f8
[ardour.git] / libs / vamp-plugins / TruePeak.cpp
1 /*
2  *  Copyright (C) 2013-2016 Robin Gareus <robin@gareus.org>
3  *  Copyright (C) 2006-2012 Fons Adriaensen <fons@linuxaudio.org>
4  *  Copyright (C) 2006 Chris Cannam
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 3 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, see <http://www.gnu.org/licenses/>.
18  */
19
20 #include <stdlib.h>
21 #include <stdio.h>
22 #include <string.h>
23 #include <math.h>
24
25 #include "TruePeak.h"
26
27 namespace TruePeakMeter {
28
29 static double sinc (double x)
30 {
31         x = fabs (x);
32         if (x < 1e-6) return 1.0;
33         x *= M_PI;
34         return sin (x) / x;
35 }
36
37 static double wind (double x)
38 {
39         x = fabs (x);
40         if (x >= 1.0) return 0.0f;
41         x *= M_PI;
42         return 0.384 + 0.500 * cos (x) + 0.116 * cos (2 * x);
43 }
44
45 Resampler_table  *Resampler_table::_list = 0;
46 Resampler_mutex   Resampler_table::_mutex;
47
48 Resampler_table::Resampler_table (double fr, unsigned int hl, unsigned int np)
49         : _next (0)
50         , _refc (0)
51         , _fr (fr)
52         , _hl (hl)
53         , _np (np)
54 {
55         unsigned int  i, j;
56         double        t;
57         float        *p;
58
59         _ctab = new float [hl * (np + 1)];
60         p = _ctab;
61         for (j = 0; j <= np; j++)
62         {
63                 t = (double) j / (double) np;
64                 for (i = 0; i < hl; i++)
65                 {
66                         p [hl - i - 1] = (float)(fr * sinc (t * fr) * wind (t / hl));
67                         t += 1;
68                 }
69                 p += hl;
70         }
71 }
72
73 Resampler_table::~Resampler_table (void)
74 {
75         delete[] _ctab;
76 }
77
78 Resampler_table *
79 Resampler_table::create (double fr, unsigned int hl, unsigned int np)
80 {
81         Resampler_table *P;
82
83         _mutex.lock ();
84         P = _list;
85         while (P)
86         {
87                 if ((fr >= P->_fr * 0.999) && (fr <= P->_fr * 1.001) && (hl == P->_hl) && (np == P->_np))
88                 {
89                         P->_refc++;
90                         _mutex.unlock ();
91                         return P;
92                 }
93                 P = P->_next;
94         }
95         P = new Resampler_table (fr, hl, np);
96         P->_refc = 1;
97         P->_next = _list;
98         _list = P;
99         _mutex.unlock ();
100         return P;
101 }
102
103 void
104 Resampler_table::destroy (Resampler_table *T)
105 {
106         Resampler_table *P, *Q;
107
108         _mutex.lock ();
109         if (T)
110         {
111                 T->_refc--;
112                 if (T->_refc == 0)
113                 {
114                         P = _list;
115                         Q = 0;
116                         while (P)
117                         {
118                                 if (P == T)
119                                 {
120                                         if (Q) Q->_next = T->_next;
121                                         else      _list = T->_next;
122                                         break;
123                                 }
124                                 Q = P;
125                                 P = P->_next;
126                         }
127                         delete T;
128                 }
129         }
130         _mutex.unlock ();
131 }
132
133 static unsigned int
134 gcd (unsigned int a, unsigned int b)
135 {
136         if (a == 0) return b;
137         if (b == 0) return a;
138         while (1)
139         {
140                 if (a > b)
141                 {
142                         a = a % b;
143                         if (a == 0) return b;
144                         if (a == 1) return 1;
145                 }
146                 else
147                 {
148                         b = b % a;
149                         if (b == 0) return a;
150                         if (b == 1) return 1;
151                 }
152         }
153         return 1;
154 }
155
156 Resampler::Resampler (void)
157         : _table (0)
158         , _nchan (0)
159         , _buff  (0)
160 {
161         reset ();
162 }
163
164 Resampler::~Resampler (void)
165 {
166         clear ();
167 }
168
169 int
170 Resampler::setup (unsigned int fs_inp,
171                   unsigned int fs_out,
172                   unsigned int nchan,
173                   unsigned int hlen)
174 {
175         if ((hlen < 8) || (hlen > 96)) return 1;
176         return setup (fs_inp, fs_out, nchan, hlen, 1.0 - 2.6 / hlen);
177 }
178
179 int
180 Resampler::setup (unsigned int fs_inp,
181                   unsigned int fs_out,
182                   unsigned int nchan,
183                   unsigned int hlen,
184                   double       frel)
185 {
186         unsigned int       g, h, k, n, s;
187         double             r;
188         float              *B = 0;
189         Resampler_table    *T = 0;
190
191         k = s = 0;
192         if (fs_inp && fs_out && nchan)
193         {
194                 r = (double) fs_out / (double) fs_inp;
195                 g = gcd (fs_out, fs_inp);
196                 n = fs_out / g;
197                 s = fs_inp / g;
198                 if ((16 * r >= 1) && (n <= 1000))
199                 {
200                         h = hlen;
201                         k = 250;
202                         if (r < 1)
203                         {
204                                 frel *= r;
205                                 h = (unsigned int)(ceil (h / r));
206                                 k = (unsigned int)(ceil (k / r));
207                         }
208                         T = Resampler_table::create (frel, h, n);
209                         B = new float [nchan * (2 * h - 1 + k)];
210                 }
211         }
212         clear ();
213         if (T)
214         {
215                 _table = T;
216                 _buff  = B;
217                 _nchan = nchan;
218                 _inmax = k;
219                 _pstep = s;
220                 return reset ();
221         }
222         else return 1;
223 }
224
225 void
226 Resampler::clear (void)
227 {
228         Resampler_table::destroy (_table);
229         delete[] _buff;
230         _buff  = 0;
231         _table = 0;
232         _nchan = 0;
233         _inmax = 0;
234         _pstep = 0;
235         reset ();
236 }
237
238 double
239 Resampler::inpdist (void) const
240 {
241         if (!_table) return 0;
242         return (int)(_table->_hl + 1 - _nread) - (double)_phase / _table->_np;
243 }
244
245 int
246 Resampler::inpsize (void) const
247 {
248         if (!_table) return 0;
249         return 2 * _table->_hl;
250 }
251
252 int
253 Resampler::reset (void)
254 {
255         if (!_table) return 1;
256
257         inp_count = 0;
258         out_count = 0;
259         inp_data = 0;
260         out_data = 0;
261         _index = 0;
262         _nread = 0;
263         _nzero = 0;
264         _phase = 0;
265         if (_table)
266         {
267                 _nread = 2 * _table->_hl;
268                 return 0;
269         }
270         return 1;
271 }
272
273 int
274 Resampler::process (void)
275 {
276         unsigned int   hl, ph, np, dp, in, nr, nz, i, n, c;
277         float          *p1, *p2;
278
279         if (!_table) return 1;
280
281         hl = _table->_hl;
282         np = _table->_np;
283         dp = _pstep;
284         in = _index;
285         nr = _nread;
286         ph = _phase;
287         nz = _nzero;
288         n = (2 * hl - nr) * _nchan;
289         p1 = _buff + in * _nchan;
290         p2 = p1 + n;
291
292         while (out_count)
293         {
294                 if (nr)
295                 {
296                         if (inp_count == 0) break;
297                         if (inp_data)
298                         {
299                                 for (c = 0; c < _nchan; c++) p2 [c] = inp_data [c];
300                                 inp_data += _nchan;
301                                 nz = 0;
302                         }
303                         else
304                         {
305                                 for (c = 0; c < _nchan; c++) p2 [c] = 0;
306                                 if (nz < 2 * hl) nz++;
307                         }
308                         nr--;
309                         p2 += _nchan;
310                         inp_count--;
311                 }
312                 else
313                 {
314                         if (out_data)
315                         {
316                                 if (nz < 2 * hl)
317                                 {
318                                         float *c1 = _table->_ctab + hl * ph;
319                                         float *c2 = _table->_ctab + hl * (np - ph);
320                                         for (c = 0; c < _nchan; c++)
321                                         {
322                                                 float *q1 = p1 + c;
323                                                 float *q2 = p2 + c;
324                                                 float s = 1e-20f;
325                                                 for (i = 0; i < hl; i++)
326                                                 {
327                                                         q2 -= _nchan;
328                                                         s += *q1 * c1 [i] + *q2 * c2 [i];
329                                                         q1 += _nchan;
330                                                 }
331                                                 *out_data++ = s - 1e-20f;
332                                         }
333                                 }
334                                 else
335                                 {
336                                         for (c = 0; c < _nchan; c++) *out_data++ = 0;
337                                 }
338                         }
339                         out_count--;
340
341                         ph += dp;
342                         if (ph >= np)
343                         {
344                                 nr = ph / np;
345                                 ph -= nr * np;
346                                 in += nr;
347                                 p1 += nr * _nchan;;
348                                 if (in >= _inmax)
349                                 {
350                                         n = (2 * hl - nr) * _nchan;
351                                         memcpy (_buff, p1, n * sizeof (float));
352                                         in = 0;
353                                         p1 = _buff;
354                                         p2 = p1 + n;
355                                 }
356                         }
357                 }
358         }
359         _index = in;
360         _nread = nr;
361         _phase = ph;
362         _nzero = nz;
363
364         return 0;
365 }
366
367 TruePeakdsp::TruePeakdsp (void)
368         : _m (0)
369         , _p (0)
370         , _res (true)
371         , _buf (NULL)
372 {
373 }
374
375 TruePeakdsp::~TruePeakdsp (void)
376 {
377         free(_buf);
378 }
379
380 void
381 TruePeakdsp::process (float const *d, int n)
382 {
383         _src.inp_count = n;
384         _src.inp_data = d;
385         _src.out_count = n * 4;
386         _src.out_data = _buf;
387         _src.process ();
388
389         float m = _res ? 0 : _m;
390         float p = _res ? 0 : _p;
391         float v;
392         float *b = _buf;
393         while (n--) {
394                 v = fabsf(*b++);
395                 if (v > m) m = v;
396                 if (v > p) p = v;
397                 v = fabsf(*b++);
398                 if (v > m) m = v;
399                 if (v > p) p = v;
400                 v = fabsf(*b++);
401                 if (v > m) m = v;
402                 if (v > p) p = v;
403                 v = fabsf(*b++);
404                 if (v > m) m = v;
405                 if (v > p) p = v;
406         }
407
408         if (_res) {
409                 _m = m;
410                 _p = p;
411                 _res = false;
412         } else {
413                 if (m > _m) { _m = m; }
414                 if (p > _p) { _p = p; }
415         }
416 }
417
418 float
419 TruePeakdsp::read (void)
420 {
421         _res = true;
422         return _m;
423 }
424
425 void
426 TruePeakdsp::read (float &m, float &p)
427 {
428         _res = true;
429         m = _m;
430         p = _p;
431 }
432
433 void
434 TruePeakdsp::reset ()
435 {
436         _res = true;
437         _m = 0;
438         _p = 0;
439 }
440
441 void
442 TruePeakdsp::init (float fsamp)
443 {
444         _src.setup(fsamp, fsamp * 4.0, 1, 24, 1.0);
445         _buf = (float*) malloc(32768 * sizeof(float));
446
447         /* q/d initialize */
448         float zero[8192];
449         for (int i = 0; i < 8192; ++i) {
450                 zero[i]= 0.0;
451         }
452         _src.inp_count = 8192;
453         _src.inp_data = zero;
454         _src.out_count = 32768;
455         _src.out_data = _buf;
456         _src.process ();
457 }
458
459 }
460
461 ///////////////////////////////////////////////////////////////////////////////
462
463 using std::string;
464 using std::vector;
465 using std::cerr;
466 using std::endl;
467 using namespace TruePeakMeter;
468
469 VampTruePeak::VampTruePeak(float inputSampleRate)
470     : Plugin(inputSampleRate)
471     , m_blockSize(0)
472 {
473 }
474
475 VampTruePeak::~VampTruePeak()
476 {
477 }
478
479 string
480 VampTruePeak::getIdentifier() const
481 {
482         return "dBTP";
483 }
484
485 string
486 VampTruePeak::getName() const
487 {
488         return "dBTP Meter";
489 }
490
491 string
492 VampTruePeak::getDescription() const
493 {
494         return "True Peak Meter (4x Oversampling)";
495 }
496
497 string
498 VampTruePeak::getMaker() const
499 {
500         return "Robin Gareus, Fons Adrianesen";
501 }
502
503 int
504 VampTruePeak::getPluginVersion() const
505 {
506         return 2;
507 }
508
509 string
510 VampTruePeak::getCopyright() const
511 {
512         return "GPL version 3 or later";
513 }
514
515 bool
516 VampTruePeak::initialise(size_t channels, size_t stepSize, size_t blockSize)
517 {
518         if (channels < getMinChannelCount() ||
519                         channels > getMaxChannelCount()) {
520                 return false;
521         }
522
523         if (blockSize > 8192) {
524                 return false;
525         }
526
527         m_blockSize = blockSize;
528
529         _meter.init (m_inputSampleRate);
530
531         return true;
532 }
533
534 void
535 VampTruePeak::reset()
536 {
537         _meter.reset ();
538 }
539
540 VampTruePeak::OutputList
541 VampTruePeak::getOutputDescriptors() const
542 {
543         OutputList list;
544
545         OutputDescriptor zc;
546         zc.identifier = "level";
547         zc.name = "TruePeak";
548         zc.description = "TruePeak (4x Oversampling)";
549         zc.unit = "dbTP";
550         zc.hasFixedBinCount = true;
551         zc.binCount = 0;
552         zc.hasKnownExtents = false;
553         zc.isQuantized = false;
554         zc.sampleType = OutputDescriptor::OneSamplePerStep;
555         list.push_back(zc);
556
557         return list;
558 }
559
560 VampTruePeak::FeatureSet
561 VampTruePeak::process(const float *const *inputBuffers,
562                       Vamp::RealTime timestamp)
563 {
564         if (m_blockSize == 0) {
565                 cerr << "ERROR: VampTruePeak::process: "
566                         << "VampTruePeak has not been initialised"
567                         << endl;
568                 return FeatureSet();
569         }
570
571         _meter.process (inputBuffers[0], m_blockSize);
572
573         // TODO return momentary
574         return FeatureSet();
575 }
576
577 VampTruePeak::FeatureSet
578 VampTruePeak::getRemainingFeatures()
579 {
580         FeatureSet returnFeatures;
581
582         float m, p;
583         _meter.read(m, p);
584
585         Feature dbtp;
586         dbtp.hasTimestamp = false;
587         dbtp.values.push_back(p);
588         returnFeatures[0].push_back(dbtp);
589
590         return returnFeatures;
591 }