1 // -----------------------------------------------------------------------
2 //  Copyright (C) 2003-2010 Fons Adriaensen <fons@linuxaudio.org>
3 //
4 //  This program is free software; you can redistribute it and/or modify
5 //  it under the terms of the GNU General Public License as published by
6 //  the Free Software Foundation; either version 2 of the License, or
7 //  (at your option) any later version.
8 //
9 //  This program is distributed in the hope that it will be useful,
10 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
11 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 //  GNU General Public License for more details.
13 //
14 //  You should have received a copy of the GNU General Public License
15 //  along with this program; if not, write to the Free Software
16 //  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
17 // -----------------------------------------------------------------------
18 
19 #include <math.h>
20 #include "zita.h"
21 
22 namespace Ms {
23 
24 enum {
25       R_DELAY, R_XOVER, R_RTLOW, R_RTMID, R_FDAMP,
26       R_EQ1FR, R_EQ1GN,
27       R_EQ2FR, R_EQ2GN,
28       R_OPMIX
29       };
30 
31 static const std::vector<ParDescr> pd = {
32       { R_DELAY, "delay", false,  0.02f,       0.100f,       0.04f },
33       { R_XOVER, "xover", true,   logf(50.0),  logf(1000.0), 200.0 },
34       { R_RTLOW, "rtlow", true,   logf(1.0),   logf(8.0),    3.0   },
35       { R_RTMID, "rtmid", true,   logf(1.0),   logf(8.0),    2.0   },
36       { R_FDAMP, "fdamp", true,   logf(1.5e3), logf(24.0e3), 6.0e3 },
37 
38       { R_EQ1FR, "eq1fr", true,   logf(40.0),  logf(2.5e3), 160.0 },
39       { R_EQ1GN, "eq1gn", false, -15.0,        15.0,        0.0   },
40 
41       { R_EQ2FR, "eq2fr", true,   logf(160.0), logf(10e3),  2.5e3 },
42       { R_EQ2GN, "eq2gn", false, -15.0,        15.0,        0.0   },
43 
44       { R_OPMIX, "opmix", false,  0.0,         1.0,         0.5   }
45       };
46 
47 //---------------------------------------------------------
48 //   Pareq
49 //---------------------------------------------------------
50 
Pareq()51 Pareq::Pareq()
52       {
53       }
54 
55 //---------------------------------------------------------
56 //   setfsamp
57 //---------------------------------------------------------
58 
setfsamp(float fsamp)59 void Pareq::setfsamp(float fsamp)
60       {
61       _fsamp = fsamp;
62       reset();
63       }
64 
65 //---------------------------------------------------------
66 //   reset
67 //---------------------------------------------------------
68 
reset()69 void Pareq::reset()
70       {
71       memset (_z1, 0, sizeof (float) * MAXCH);
72       memset (_z2, 0, sizeof (float) * MAXCH);
73       }
74 
75 //---------------------------------------------------------
76 //   prepare
77 //---------------------------------------------------------
78 
prepare(int nsamp)79 void Pareq::prepare(int nsamp)
80       {
81       bool  upd = false;
82 
83       uint16_t touchValue = _touch0;
84       if (_touch1 != touchValue) {
85             float g = _g0;
86             float f = _f0;
87             if (g != _g1) {
88                   upd = true;
89                   if (g > 2 * _g1)
90                        _g1 *= 2;
91                   else if (_g1 > 2 * g)
92                         _g1 /= 2;
93                   else
94                         _g1 = g;
95                   }
96             if (f != _f1) {
97                   upd = true;
98                   if (f > 2 * _f1)
99                         _f1 *= 2;
100                   else if (_f1 > 2 * f)
101                         _f1 /= 2;
102                   else
103                         _f1 = f;
104                   }
105             if (upd)  {
106                   if ((_state == BYPASS) && (_g1 == 1)) {
107                         calcpar1 (0, _g1, _f1);
108                         }
109                   else {
110                         _state = SMOOTH;
111                         calcpar1 (nsamp, _g1, _f1);
112                         }
113                   }
114             else {
115                   _touch1 = touchValue;
116                   if (fabs (_g1 - 1) < 0.001f) {
117                         _state = BYPASS;
118                         reset ();
119                         }
120                   else {
121                         _state = STATIC;
122                         }
123                   }
124             }
125       }
126 
127 //---------------------------------------------------------
128 //   calcpar1
129 //---------------------------------------------------------
130 
calcpar1(int nsamp,float g,float f)131 void Pareq::calcpar1(int nsamp, float g, float f)
132       {
133       f *= float (M_PI) / _fsamp;
134 
135       float b = 2 * f / sqrtf (g);
136       float gg = 0.5f * (g - 1);
137       float c1 = -cosf (2 * f);
138       float c2 = (1 - b) / (1 + b);
139       if (nsamp) {
140             _dc1 = (c1 - _c1) / nsamp + 1e-30f;
141             _dc2 = (c2 - _c2) / nsamp + 1e-30f;
142             _dgg = (gg - _gg) / nsamp + 1e-30f;
143             }
144       else {
145             _c1 = c1;
146             _c2 = c2;
147             _gg = gg;
148             }
149       }
150 
151 //---------------------------------------------------------
152 //   process1
153 //---------------------------------------------------------
154 
process1(int nsamp,float * data)155 void Pareq::process1(int nsamp, float* data)
156       {
157       int nchan = 2;
158       float c1 = _c1;
159       float c2 = _c2;
160       float gg = _gg;
161 
162       if (_state == SMOOTH) {
163             for (int i = 0; i < nchan; i++) {
164                   float z1 = _z1 [i];
165                   float z2 = _z2 [i];
166                   c1       = _c1;
167                   c2       = _c2;
168                   gg = _gg;
169                   for (int j = 0; j < nsamp; j++) {
170                         c1 += _dc1;
171                         c2 += _dc2;
172                         gg += _dgg;
173                         float* p = data + j * 2 + i;
174                         float x = *p;
175                         float y = x - c2 * z2;
176                         *p = x - gg * (z2 + c2 * y - x);
177                         y -= c1 * z1;
178                         z2 = z1 + c1 * y;
179                         z1 = y + 1e-20f;
180                         }
181                   _z1 [i] = z1;
182                   _z2 [i] = z2;
183                   }
184             _c1 = c1;
185             _c2 = c2;
186             _gg = gg;
187             }
188       else {
189             for (int i = 0; i < nchan; i++) {
190                   float z1 = _z1 [i];
191                   float z2 = _z2 [i];
192                   for (int j = 0; j < nsamp; j++) {
193                         float* p = data + j * 2 + i;
194                         float x = *p;
195                         float y = x - c2 * z2;
196                         *p = x - gg * (z2 + c2 * y - x);
197                         y -= c1 * z1;
198                         z2 = z1 + c1 * y;
199                         z1 = y + 1e-20f;
200                         }
201                   _z1 [i] = z1;
202                   _z2 [i] = z2;
203                   }
204             }
205       }
206 
~Diff1()207 Diff1::~Diff1()
208       {
209       fini();
210       }
211 
init(int size,float c)212 void Diff1::init (int size, float c)
213       {
214       _size = size;
215       _line = new float [size];
216       memset (_line, 0, size * sizeof (float));
217       _i = 0;
218       _c = c;
219       }
220 
fini()221 void Diff1::fini()
222       {
223       delete[] _line;
224       _size = 0;
225       _line = 0;
226       }
227 
Delay()228 Delay::Delay()
229    : _i(0), _size (0), _line (0)
230       {
231       }
232 
~Delay()233 Delay::~Delay()
234       {
235       fini();
236       }
237 
init(int size)238 void Delay::init (int size)
239       {
240       _size = size;
241       _line = new float [size];
242       memset (_line, 0, size * sizeof (float));
243       _i = 0;
244       }
245 
fini()246 void Delay::fini ()
247       {
248       delete[] _line;
249       _size = 0;
250       _line = 0;
251       }
252 
Vdelay()253 Vdelay::Vdelay ()
254    : _ir(0), _iw(0), _size (0), _line (0)
255       {
256       }
257 
~Vdelay()258 Vdelay::~Vdelay ()
259       {
260       fini();
261       }
262 
init(int size)263 void Vdelay::init (int size)
264       {
265       _size = size;
266       _line = new float [size];
267       memset (_line, 0, size * sizeof (float));
268       _ir = 0;
269       _iw = 0;
270       }
271 
fini()272 void Vdelay::fini ()
273       {
274       delete[] _line;
275       _size = 0;
276       _line = 0;
277       }
278 
set_delay(int del)279 void Vdelay::set_delay (int del)
280       {
281       _ir = _iw - del;
282       if (_ir < 0)
283             _ir += _size;
284       }
285 
set_params(float del,float tmf,float tlo,float wlo,float thi,float chi)286 void Filt1::set_params (float del, float tmf, float tlo, float wlo, float thi, float chi)
287       {
288       _gmf = powf (0.001f, del / tmf);
289       _glo = powf (0.001f, del / tlo) / _gmf - 1.0f;
290       _wlo = wlo;
291       float g    = powf (0.001f, del / thi) / _gmf;
292       float t    = (1 - g * g) / (2 * g * g * chi);
293       _whi = (sqrtf (1 + 4 * t) - 1) / (2 * t);
294       }
295 
296 float ZitaReverb::_tdiff1 [8] = {
297       20346e-6f,
298       24421e-6f,
299       31604e-6f,
300       27333e-6f,
301       22904e-6f,
302       29291e-6f,
303       13458e-6f,
304       19123e-6f
305       };
306 
307 float ZitaReverb::_tdelay [8] = {
308       153129e-6f,
309       210389e-6f,
310       127837e-6f,
311       256891e-6f,
312       174713e-6f,
313       192303e-6f,
314       125000e-6f,
315       219991e-6f
316       };
317 
~ZitaReverb()318 ZitaReverb::~ZitaReverb ()
319       {
320       fini();
321       }
322 
init(float fsamp)323 void ZitaReverb::init(float fsamp)
324       {
325       _fsamp = fsamp;
326       _cntA1 = 1;
327       _cntA2 = 0;
328       _cntB1 = 1;
329       _cntB2 = 0;
330       _cntC1 = 1;
331       _cntC2 = 0;
332 
333       _ipdel = 0.04f;
334       _xover = 200.0f;
335       _rtlow = 1.4f;
336       _rtmid = 2.0f;
337       _fdamp = 3e3f;
338       _opmix = 0.33f;
339 
340       _g0 = _d0 = 0;
341       _g1 = _d1 = 0;
342 
343       _vdelay0.init ((int)(0.1f * _fsamp));
344       _vdelay1.init ((int)(0.1f * _fsamp));
345       for (int i = 0; i < 8; i++) {
346             int k1 = (int)(floorf (_tdiff1 [i] * _fsamp + 0.5f));
347             int k2 = (int)(floorf (_tdelay [i] * _fsamp + 0.5f));
348             _diff1 [i].init (k1, (i & 1) ? -0.6f : 0.6f);
349             _delay [i].init (k2 - k1);
350             }
351 
352       _pareq1.setfsamp(fsamp);
353       _pareq2.setfsamp(fsamp);
354       _pareq1.setparam(160.0, 0.0);
355       _pareq2.setparam(2.5e3, 0.0);
356 
357       _fragm = 1024;
358       _nsamp = 0;
359       }
360 
361 
fini()362 void ZitaReverb::fini ()
363       {
364       for (int i = 0; i < 8; i++)
365             _delay [i].fini ();
366       }
367 
368 //---------------------------------------------------------
369 //   prepare
370 //---------------------------------------------------------
371 
prepare(int nfram)372 void ZitaReverb::prepare (int nfram)
373       {
374       int a = _cntA1;
375       int b = _cntB1;
376       int c = _cntC1;
377 
378       _d0 = _d1 = 0;
379 
380       if (a != _cntA2) {
381             int k = (int)(floorf ((_ipdel - 0.020f) * _fsamp + 0.5f));
382             _vdelay0.set_delay (k);
383             _vdelay1.set_delay (k);
384             _cntA2 = a;
385             }
386 
387       if (b != _cntB2) {
388             float wlo = 6.2832f * _xover / _fsamp;
389             float chi = 0.f;
390             if (_fdamp > 0.49f * _fsamp)
391                   chi = 2;
392             else
393                   chi = 1 - cosf (6.2832f * _fdamp / _fsamp);
394             for (int i = 0; i < 8; i++) {
395                   _filt1 [i].set_params (_tdelay [i], _rtmid, _rtlow, wlo, 0.5f * _rtmid, chi);
396                   }
397             _cntB2 = b;
398             }
399 
400       if (c != _cntC2) {
401             float t0 = (1 - _opmix) * (1 + _opmix);
402             float t1 = 0.7f * _opmix * (2 - _opmix) / sqrtf (_rtmid);
403             _d0 = (t0 - _g0) / nfram;
404             _d1 = (t1 - _g1) / nfram;
405             _cntC2 = c;
406             }
407 
408       _pareq1.prepare (nfram);
409       _pareq2.prepare (nfram);
410       }
411 
412 //---------------------------------------------------------
413 //   process
414 //---------------------------------------------------------
415 
process(int nfram,float * inp,float * out)416 void ZitaReverb::process (int nfram, float* inp, float* out)
417       {
418       float t = 0.f, g = 0.f, x0 = 0.f, x1 = 0.f, x2 = 0.f, x3 = 0.f, x4 = 0.f, x5 = 0.f, x6 = 0.f, x7 = 0.f;
419       g = sqrtf (0.125f);
420 
421 
422       while (nfram) {
423             if (!_nsamp) {
424                   prepare(_fragm);
425                   _nsamp = _fragm;
426                   }
427 
428             int k = _nsamp < nfram ? _nsamp : nfram;
429 
430             float* p0 = inp;
431             float* p1 = inp + 1;
432             float* q0 = out;
433             float* q1 = out + 1;
434 
435             for (int i = 0; i < k * 2; i += 2) {
436                   _vdelay0.write (p0 [i]);
437                   _vdelay1.write (p1 [i]);
438 
439                   t = 0.3f * _vdelay0.read ();
440                   x0 = _diff1 [0].process (_delay [0].read () + t);
441                   x1 = _diff1 [1].process (_delay [1].read () + t);
442                   x2 = _diff1 [2].process (_delay [2].read () - t);
443                   x3 = _diff1 [3].process (_delay [3].read () - t);
444                   t = 0.3f * _vdelay1.read ();
445                   x4 = _diff1 [4].process (_delay [4].read () + t);
446                   x5 = _diff1 [5].process (_delay [5].read () + t);
447                   x6 = _diff1 [6].process (_delay [6].read () - t);
448                   x7 = _diff1 [7].process (_delay [7].read () - t);
449 
450                   t = x0 - x1; x0 += x1;  x1 = t;
451                   t = x2 - x3; x2 += x3;  x3 = t;
452                   t = x4 - x5; x4 += x5;  x5 = t;
453                   t = x6 - x7; x6 += x7;  x7 = t;
454                   t = x0 - x2; x0 += x2;  x2 = t;
455                   t = x1 - x3; x1 += x3;  x3 = t;
456                   t = x4 - x6; x4 += x6;  x6 = t;
457                   t = x5 - x7; x5 += x7;  x7 = t;
458                   t = x0 - x4; x0 += x4;  x4 = t;
459                   t = x1 - x5; x1 += x5;  x5 = t;
460                   t = x2 - x6; x2 += x6;  x6 = t;
461                   t = x3 - x7; x3 += x7;  x7 = t;
462 
463                   _g1 += _d1;
464 
465                   q0 [i] = _g1 * (x1 + x2);
466                   q1 [i] = _g1 * (x1 - x2);
467 
468                   _delay [0].write (_filt1 [0].process (g * x0));
469                   _delay [1].write (_filt1 [1].process (g * x1));
470                   _delay [2].write (_filt1 [2].process (g * x2));
471                   _delay [3].write (_filt1 [3].process (g * x3));
472                   _delay [4].write (_filt1 [4].process (g * x4));
473                   _delay [5].write (_filt1 [5].process (g * x5));
474                   _delay [6].write (_filt1 [6].process (g * x6));
475                   _delay [7].write (_filt1 [7].process (g * x7));
476                   }
477             _pareq1.process (k, out);
478             _pareq2.process (k, out);
479 
480             for (int i = 0; i < k; i++) {
481                   *out++ += _g0 * *inp++;
482                   *out++ += _g0 * *inp++;
483                   _g0 += _d0;
484                   }
485             nfram  -= k;
486             _nsamp -= k;
487             }
488       }
489 
setNValue(int idx,double value)490 void ZitaReverb::setNValue(int idx, double value)
491       {
492       float fvalue = static_cast<float>(value);
493       if (!std::isfinite(fvalue)) {
494             fvalue = pd[idx].init;
495       }
496       float min = (pd[idx].log ? pow(M_E, pd[idx].min) : pd[idx].min);
497       float max = (pd[idx].log ? pow(M_E, pd[idx].max) : pd[idx].max);
498       fvalue = std::max(min, fvalue);
499       fvalue = std::min(max, fvalue);
500       switch (idx) {
501             case R_DELAY: set_delay(fvalue); break;
502             case R_XOVER: set_xover(fvalue); break;
503             case R_RTLOW: set_rtlow(fvalue); break;
504             case R_RTMID: set_rtmid(fvalue); break;
505             case R_FDAMP: set_fdamp(fvalue); break;
506             case R_EQ1FR: set_eq1fr(fvalue); break;
507             case R_EQ1GN: set_eq1gn(fvalue); break;
508             case R_EQ2FR: set_eq2fr(fvalue); break;
509             case R_EQ2GN: set_eq2gn(fvalue); break;
510             case R_OPMIX: set_opmix(fvalue); break;
511             }
512       }
513 
514 //---------------------------------------------------------
515 //   value
516 //    return normalized value 0-1.0
517 //---------------------------------------------------------
518 
nvalue(int idx) const519 double ZitaReverb::nvalue(int idx) const
520       {
521       double v = 0.0;
522       switch (idx) {
523             case R_DELAY: v = delay(); break;
524             case R_XOVER: v = xover(); break;
525             case R_RTLOW: v = rtlow(); break;
526             case R_RTMID: v = rtmid(); break;
527             case R_FDAMP: v = fdamp(); break;
528             case R_EQ1FR: v = eq1fr(); break;
529             case R_EQ1GN: v = eq1gn(); break;
530             case R_EQ2FR: v = eq2fr(); break;
531             case R_EQ2GN: v = eq2gn(); break;
532             case R_OPMIX: v = opmix(); break;
533             }
534       return v;
535       }
536 
537 //---------------------------------------------------------
538 //   parDescr
539 //---------------------------------------------------------
540 
parDescr() const541 const std::vector<ParDescr>& ZitaReverb::parDescr() const
542       {
543       return pd;
544       }
545 
546 //---------------------------------------------------------
547 //   state
548 //---------------------------------------------------------
549 
state() const550 SynthesizerGroup ZitaReverb::state() const
551       {
552       SynthesizerGroup g;
553       g.setName(name());
554 
555       for (const ParDescr& d : pd)
556             g.push_back(IdValue(d.id, QString("%1").arg(value(d.id))));
557       return g;
558       }
559 
560 //---------------------------------------------------------
561 //   setState
562 //---------------------------------------------------------
563 
setState(const SynthesizerGroup & g)564 void ZitaReverb::setState(const SynthesizerGroup& g)
565       {
566       for (const IdValue& v : g)
567             setValue(v.id, v.data.toDouble());
568       }
569 }
570 
571