basics.lib

//#################################### basics.lib ########################################
// A library of basic elements. Its official prefix is `ba`.
//
// #### References
// * <https://github.com/grame-cncm/faustlibraries/blob/master/basics.lib>
//########################################################################################
// A library of basic elements for Faust organized in 5 sections:
//
// * Conversion Tools
// * Counters and Time/Tempo Tools
// * Array Processing/Pattern Matching
// * Selectors (Conditions)
// * Other Tools (Misc)

//########################################################################################

/************************************************************************
************************************************************************
FAUST library file, GRAME section

Except where noted otherwise, Copyright (C) 2003-2017 by GRAME,
Centre National de Creation Musicale.
----------------------------------------------------------------------
GRAME LICENSE

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2.1 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA.

EXCEPTION TO THE LGPL LICENSE : As a special exception, you may create a
larger FAUST program which directly or indirectly imports this library
file and still distribute the compiled code generated by the FAUST
compiler, or a modified version of this compiled code, under your own
copyright and license. This EXCEPTION TO THE LGPL LICENSE explicitly
grants you the right to freely choose the license for the resulting
compiled code. In particular the resulting compiled code has no obligation
to be LGPL or GPL. For example you are free to choose a commercial or
closed source license or any other license if you decide so.
************************************************************************
************************************************************************/

ma = library("maths.lib");
ro = library("routes.lib");
ba = library("basics.lib"); // for compatible copy/paste out of this file
fi = library("filters.lib");
it = library("interpolators.lib");
si = library("signals.lib");

declare name "Faust Basic Element Library";
declare version "1.21.0";

//=============================Conversion Tools===========================================
//========================================================================================

//-------`(ba.)samp2sec`----------
// Converts a number of samples to a duration in seconds at the current sampling rate (see `ma.SR`).
// `samp2sec` is a standard Faust function.
//
// #### Usage
//
// ```
// samp2sec(n) : _
// ```
//
// Where:
//
// * `n`: number of samples
//----------------------------
samp2sec(n) = n/ma.SR;


//-------`(ba.)sec2samp`----------
// Converts a duration in seconds to a number of samples at the current sampling rate (see `ma.SR`).
// `samp2sec` is a standard Faust function.
//
// #### Usage
//
// ```
// sec2samp(d) : _
// ```
//
// Where:
//
// * `d`: duration in seconds
//----------------------------
sec2samp(d) = d*ma.SR;


//-------`(ba.)db2linear`----------
// dB-to-linear value converter. It can be used to convert an amplitude in dB to a linear gain ]0-N].
// `db2linear` is a standard Faust function.
//
// #### Usage
//
// ```
// db2linear(l) : _
// ```
//
// Where:
//
// * `l`: amplitude in dB
//-----------------------------
db2linear(l) = pow(10.0, l/20.0);


//-------`(ba.)linear2db`----------
// linea-to-dB value converter. It can be used to convert a linear gain ]0-N] to an amplitude in dB.
// `linear2db` is a standard Faust function.
//
// #### Usage
//
// ```
// linear2db(g) : _
// ```
//
// Where:
//
// * `g`: a linear gain
//-----------------------------
linear2db(g) = 20.0*log10(max(ma.MIN, g));


//----------`(ba.)lin2LogGain`------------------
// Converts a linear gain (0-1) to a log gain (0-1).
//
// #### Usage
//
// ```
// lin2LogGain(n) : _
// ```
//
// Where:
//
// * `n`: the linear gain
//---------------------------------------------
lin2LogGain(n) = n*n;


//----------`(ba.)log2LinGain`------------------
// Converts a log gain (0-1) to a linear gain (0-1).
//
// #### Usage
//
// ```
// log2LinGain(n) : _
// ```
//
// Where:
//
// * `n`: the log gain
//---------------------------------------------
log2LinGain(n) = sqrt(n);


// end GRAME section
//########################################################################################
/************************************************************************
FAUST library file, jos section

Except where noted otherwise, The Faust functions below in this
section are Copyright (C) 2003-2017 by Julius O. Smith III <jos@ccrma.stanford.edu>
([jos](http://ccrma.stanford.edu/~jos/)), and released under the
(MIT-style) [STK-4.3](#stk-4.3-license) license.

The MarkDown comments in this section are Copyright 2016-2017 by Romain
Michon and Julius O. Smith III, and are released under the
[CCA4I](https://creativecommons.org/licenses/by/4.0/) license (TODO: if/when Romain agrees)

************************************************************************/

//-------`(ba.)tau2pole`----------
// Returns a real pole giving exponential decay.
// Note that t60 (time to decay 60 dB) is ~6.91 time constants.
// `tau2pole` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : smooth(tau2pole(tau)) : _
// ```
//
// Where:
//
// * `tau`: time-constant in seconds
//-----------------------------
// tau2pole(tau) = exp(-1.0/(tau*ma.SR));

tau2pole(tau) = ba.if(clipCond, 0.0, exp(-1.0/(tauCenterClipped*float(ma.SR))))
with {
    clipCond = abs(tau)<ma.EPSILON;
    tauCenterClipped = ba.if(clipCond, 1.0, tau); // 1.0 can be any nonzero value (not used)
};


//-------`(ba.)pole2tau`----------
// Returns the time-constant, in seconds, corresponding to the given real,
// positive pole in (0-1).
// `pole2tau` is a standard Faust function.
//
// #### Usage
//
// ```
// pole2tau(pole) : _
// ```
//
// Where:
//
// * `pole`: the pole
//-----------------------------
pole2tau(pole) = -1.0/(log(max(ma.MIN, pole))*ma.SR);


//-------`(ba.)midikey2hz`----------
// Converts a MIDI key number to a frequency in Hz (MIDI key 69 = A440).
// `midikey2hz` is a standard Faust function.
//
// #### Usage
//
// ```
// midikey2hz(mk) : _
// ```
//
// Where:
//
// * `mk`: the MIDI key number
//-----------------------------
midikey2hz(mk) = 440.0*pow(2.0, (mk-69.0)/12.0);


//-------`(ba.)hz2midikey`----------
// Converts a frequency in Hz to a MIDI key number (MIDI key 69 = A440).
// `hz2midikey` is a standard Faust function.
//
// #### Usage
//
// ```
// hz2midikey(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//-----------------------------
hz2midikey(freq) = 12.0*ma.log2(freq/440.0) + 69.0;


//-------`(ba.)semi2ratio`----------
// Converts semitones in a frequency multiplicative ratio.
// `semi2ratio` is a standard Faust function.
//
// #### Usage
//
// ```
// semi2ratio(semi) : _
// ```
//
// Where:
//
// * `semi`: number of semitone
//-----------------------------
semi2ratio(semi) = pow(2.0, semi/12.0);


//-------`(ba.)ratio2semi`----------
// Converts a frequency multiplicative ratio in semitones.
// `ratio2semi` is a standard Faust function.
//
// #### Usage
//
// ```
// ratio2semi(ratio) : _
// ```
//
// Where:
//
// * `ratio`: frequency multiplicative ratio
//-----------------------------
ratio2semi(ratio) = 12.0*log(ratio)/log(2.0);


//-------`(ba.)cent2ratio`----------
// Converts cents in a frequency multiplicative ratio.
//
// #### Usage
//
// ```
// cent2ratio(cent) : _
// ```
//
// Where:
//
// * `cent`: number of cents
//-----------------------------
cent2ratio(cent) = pow(2.0, cent/1200.0);


//-------`(ba.)ratio2cent`----------
// Converts a frequency multiplicative ratio in cents.
//
// #### Usage
//
// ```
// ratio2cent(ratio) : _
// ```
//
// Where:
//
// * `ratio`: frequency multiplicative ratio
//-----------------------------
ratio2cent(ratio) = 1200.0*log(ratio)/log(2.0);


//-------`(ba.)pianokey2hz`----------
// Converts a piano key number to a frequency in Hz (piano key 49 = A440).
//
// #### Usage
//
// ```
// pianokey2hz(pk) : _
// ```
//
// Where:
//
// * `pk`: the piano key number
//-----------------------------
pianokey2hz(pk) = 440.0*pow(2.0, (pk-49.0)/12.0);


//-------`(ba.)hz2pianokey`----------
// Converts a frequency in Hz to a piano key number (piano key 49 = A440).
//
// #### Usage
//
// ```
// hz2pianokey(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//-----------------------------
hz2pianokey(freq) = 12.0*ma.log2(freq/440.0) + 49.0;


// end jos section
//########################################################################################
/************************************************************************
FAUST library file, GRAME section 2
************************************************************************/

//==============================Counters and Time/Tempo Tools=============================
//========================================================================================

//----------------------------`(ba.)counter`------------------------------
// Starts counting 0, 1, 2, 3..., and raise the current integer value
// at each upfront of the trigger.
//
// #### Usage
//
// ```
// counter(trig) : _
// ```
//
// Where:
//
// * `trig`: the trigger signal, each upfront will move the counter to the next integer
//-----------------------------------------------------------------------------
declare counter author "Stephane Letz";

counter(trig) = upfront(trig) : + ~ _ with { upfront(x) = x > x'; };


//----------------------------`(ba.)countdown`------------------------------
// Starts counting down from n included to 0. While trig is 1 the output is n.
// The countdown starts with the transition of trig from 1 to 0. At the end
// of the countdown the output value will remain at 0 until the next trig.
// `countdown` is a standard Faust function.
//
// #### Usage
//
// ```
// countdown(n,trig) : _
// ```
//
// Where:
//
// * `n`: the starting point of the countdown
// * `trig`: the trigger signal (1: start at `n`; 0: decrease until 0)
//-----------------------------------------------------------------------------
countdown(n, trig) = \(c).(if(trig>0, n, max(0, c-1))) ~ _;


//----------------------------`(ba.)countup`--------------------------------
// Starts counting up from 0 to n included. While trig is 1 the output is 0.
// The countup starts with the transition of trig from 1 to 0. At the end
// of the countup the output value will remain at n until the next trig.
// `countup` is a standard Faust function.
//
// #### Usage
//
// ```
// countup(n,trig) : _
// ```
//
// Where:
//
// * `n`: the maximum count value
// * `trig`: the trigger signal (1: start at 0; 0: increase until `n`)
//-----------------------------------------------------------------------------
countup(n, trig) = \(c).(if(trig>0, 0, min(n, c+1))) ~ _;


//--------------------`(ba.)sweep`--------------------------
// Counts from 0 to `period-1` repeatedly, generating a
// sawtooth waveform, like `os.lf_rawsaw`,
// starting at 1 when `run` transitions from 0 to 1.
// Outputs zero while `run` is 0.
//
// #### Usage
//
// ```
// sweep(period,run) : _
// ```
//-----------------------------------------------------------------
declare sweep author "Jonatan Liljedahl";

sweep = %(int(*:max(1)))~+(1);


//-------`(ba.)time`----------
// A simple timer that counts every samples from the beginning of the process.
// `time` is a standard Faust function.
//
// #### Usage
//
// ```
// time : _
// ```
//------------------------
time = +(1)~_ : mem;


//-------`(ba.)ramp`----------
// A linear ramp with a slope of '(+/-)1/n' samples to reach the next target value.
//
// #### Usage
//
// ```
// _ : ramp(n) : _
// ```
// Where:
//
// * `n`: number of samples to increment/decrement the value by one
//------------------------
ramp = case {
    (0) => _;
    (n) => \(y,x).(if(y+1.0/n < x, y+1.0/n, if(y-1.0/n > x, y-1.0/n, x))) ~ _;
};


//-------`(ba.)line`----------
// A ramp interpolator that generates a linear transition to reach a target value:
//
//  - the interpolation process restarts each time a new and distinct input value is received
//  - it utilizes 'n' samples to achieve the transition to the target value
//  - after reaching the target value, the output value is maintained.
//
// #### Usage
//
// ```
// _ : line(n) : _
// ```
// Where:
//
// * `n`: number of samples to reach the new target received at its input
//------------------------
line(n, x) = state ~ (_,_) : !,_
with {
    state(t, c) = nt,nc
    with {
        nt = ba.if(x != x', n, t-1);
        nc = ba.if(nt > 0, c + (x - c)/nt, x);
    };
};


//-------`(ba.)tempo`----------
// Converts a tempo in BPM into a number of samples.
//
// #### Usage
//
// ```
// tempo(t) : _
// ```
//
// Where:
//
// * `t`: tempo in BPM
//------------------------
tempo(t) = (60*ma.SR)/t;


//-------`(ba.)period`----------
// Basic sawtooth wave of period `p`.
//
// #### Usage
//
// ```
// period(p) : _
// ```
//
// Where:
//
// * `p`: period as a number of samples
//------------------------
// NOTE: may be this should go in oscillators.lib
period(p) = %(int(p))~+(1');


//-------`(ba.)spulse`----------
// Produces a single pulse of n samples when trig goes from 0 to 1.
//
// #### Usage
//
// ```
// spulse(n,trig) : _
// ```
//
// Where:
//
// * `n`: pulse length as a number of samples
// * `trig`: the trigger signal (1: start the pulse)
//------------------------
spulse(n, trig) = trig : trigger(n)
with {
    upfront(x) = x > x';    // detect rising edge
    reset = upfront(trig);  // reset signal
    decay(n, x) = ba.if(reset, 0, x - (x>0.0)/n);   // decay the pulse, possibly using reset to restart the pulse
    release(n) = + ~ decay(n);                      
    trigger(n) = upfront : release(n) : >(0.0);    
};

//-------`(ba.)pulse`----------
// Pulses (like 10000) generated at period `p`.
//
// #### Usage
//
// ```
// pulse(p) : _
// ```
//
// Where:
//
// * `p`: period as a number of samples
//------------------------
// NOTE: may be this should go in oscillators.lib
pulse(p) = period(p) : \(x).(x <= x');


//-------`(ba.)pulsen`----------
// Pulses (like 11110000) of length `n` generated at period `p`.
//
// #### Usage
//
// ```
// pulsen(n,p) : _
// ```
//
// Where:
//
// * `n`: pulse length as a number of samples
// * `p`: period as a number of samples
//------------------------
// NOTE: may be this should go in oscillators.lib
pulsen(n,p) = period(p)<n;


//-----------------------`(ba.)cycle`---------------------------
// Split nonzero input values into `n` cycles.
//
// #### Usage
//
// ```
// _ : cycle(n) : si.bus(n)
// ```
//
// Where:
//
// * `n`: the number of cycles/output signals
//---------------------------------------------------------
declare cycle author "Mike Olsen";

cycle(n) = _ <: par(i,n,resetCtr(n,(i+1)));


//-------`(ba.)beat`----------
// Pulses at tempo `t` in BPM.
// `beat` is a standard Faust function.
//
// #### Usage
//
// ```
// beat(t) : _
// ```
//
// Where:
//
// * `t`: tempo in BPM
//------------------------
beat(t) = pulse(tempo(t));


//----------------------------`(ba.)pulse_countup`-----------------------------------
// Starts counting up pulses. While trig is 1 the output is
// counting up, while trig is 0 the counter is reset to 0.
//
// #### Usage
//
// ```
// _ : pulse_countup(trig) : _
// ```
//
// Where:
//
// * `trig`: the trigger signal (1: start at next pulse; 0: reset to 0)
//------------------------------------------------------------------------------
declare pulse_countup author "Vince";

pulse_countup(trig) = + ~ _ * trig;


//----------------------------`(ba.)pulse_countdown`---------------------------------
// Starts counting down pulses. While trig is 1 the output is
// counting down, while trig is 0 the counter is reset to 0.
//
// #### Usage
//
// ```
// _ : pulse_countdown(trig) : _
// ```
//
// Where:
//
// * `trig`: the trigger signal (1: start at next pulse; 0: reset to 0)
//------------------------------------------------------------------------------
declare pulse_countdown author "Vince";

pulse_countdown(trig) = - ~ _ * trig;


//----------------------------`(ba.)pulse_countup_loop`------------------------------
// Starts counting up pulses from 0 to n included. While trig is 1 the output is
// counting up, while trig is 0 the counter is reset to 0. At the end
// of the countup (n) the output value will be reset to 0.
//
// #### Usage
//
// ```
// _ : pulse_countup_loop(n,trig) : _
// ```
//
// Where:
//
// * `n`: the highest number of the countup (included) before reset to 0
// * `trig`: the trigger signal (1: start at next pulse; 0: reset to 0)
//------------------------------------------------------------------------------
declare pulse_countup_loop author "Vince";

pulse_countup_loop(n, trig) = + ~ cond(n)*trig
with {
    cond(n, x) = x * (x <= n);
};


//----------------------------`(ba.)pulse_countdown_loop`----------------------------
// Starts counting down pulses from 0 to n included. While trig is 1 the output
// is counting down, while trig is 0 the counter is reset to 0. At the end
// of the countdown (n) the output value will be reset to 0.
//
// #### Usage
//
// ```
// _ : pulse_countdown_loop(n,trig) : _
// ```
//
// Where:
//
// * `n`: the highest number of the countup (included) before reset to 0
// * `trig`: the trigger signal (1: start at next pulse; 0: reset to 0)
//------------------------------------------------------------------------------
declare pulse_countdown_loop author "Vince";

pulse_countdown_loop(n, trig) = - ~ cond(n)*trig
with {
    cond(n, x) = x * (x >= n);
};


//-----------------------`(ba.)resetCtr`------------------------
// Function that lets through the mth impulse out of
// each consecutive group of `n` impulses.
//
// #### Usage
//
// ```
// _ : resetCtr(n,m) : _
// ```
//
// Where:
//
// * `n`: the total number of impulses being split
// * `m`: index of impulse to allow to be output
//---------------------------------------------------------
declare resetCtr author "Mike Olsen";

resetCtr(n,m) = _ <: (_,pulse_countup_loop(n-1,1)) : (_,(_==m)) : *;


//===============================Array Processing/Pattern Matching========================
//========================================================================================

//---------------------------------`(ba.)count`---------------------------------
// Count the number of elements of list l.
// `count` is a standard Faust function.
//
// #### Usage
//
// ```
// count(l)
// count((10,20,30,40)) -> 4
// ```
//
// Where:
//
// * `l`: list of elements
//-----------------------------------------------------------------------------
count((xs, xxs)) = 1 + count(xxs);
count(xx) = 1;


//-------------------------------`(ba.)take`-----------------------------------
// Take an element from a list.
// `take` is a standard Faust function.
//
// #### Usage
//
// ```
// take(P,l)
// take(3,(10,20,30,40)) -> 30
// ```
//
// Where:
//
// * `P`: position (int, known at compile time, P > 0)
// * `l`: list of elements
//-----------------------------------------------------------------------------
take(1, (xs, xxs))  = xs;
take(1, xs)         = xs;
take(N, (xs, xxs)) = take(N-1, xxs);


//----------------------------`(ba.)subseq`--------------------------------
// Extract a part of a list.
//
// #### Usage
//
// ```
// subseq(l, P, N)
// subseq((10,20,30,40,50,60), 1, 3) -> (20,30,40)
// subseq((10,20,30,40,50,60), 4, 1) -> 50
// ```
//
// Where:
//
// * `l`: list
// * `P`: start point (int, known at compile time, 0: begin of list)
// * `N`: number of elements (int, known at compile time)
//
// #### Note:
//
// Faust doesn't have proper lists. Lists are simulated with parallel
// compositions and there is no empty list.
//-----------------------------------------------------------------------------
subseq((head, tail), 0, 1) = head;
subseq((head, tail), 0, N) = head, subseq(tail, 0, N-1);
subseq((head, tail), P, N) = subseq(tail, P-1, N);
subseq(head, 0, N)         = head;


//============================Function tabulation=========================================
// The purpose of function tabulation is to speed up the computation of heavy functions over an interval, 
// so that the computation at runtime can be faster than directly using the function. 
// Two techniques are implemented: 
//
// * `tabulate` computes the function in a table and read the points using interpolation. `tabulateNd` is the N dimensions version of `tabulate`
//
// * `tabulate_chebychev` uses Chebyshev polynomial approximation
//
// #### Comparison program example 
// ```
///* Both tabulate() and tabulate_chebychev() create rdtable of size = 200, both use */
///* cubic polynomials, so this comparison is more or less fair. */
// process = line(50000, r0, r1) <: FX-tb,FX-ch : par(i, 2, maxerr)
// with {
//    C = 0;
//    FX = sin; 
//    NX = 50; 
//    CD = 3;
//    r0 = 0;
//    r1 = ma.PI;
//    tb(x) = ba.tabulate(C, FX, NX*(CD+1), r0, r1, x).cub;
//    ch(x) = ba.tabulate_chebychev(C, FX, NX, CD, r0, r1, x);
//    maxerr = abs : max ~ _;
//    line(n, x0, x1) = x0 + (ba.time%n)/n * (x1-x0);
// };
// ```


//-------`(ba.)tabulate`----------
// Tabulate a 1D function over the range [r0, r1] for access via nearest-value, linear, cubic interpolation.
// In other words, the uniformly tabulated function can be evaluated using interpolation of order 0 (none), 1 (linear), or 3 (cubic).
//
// #### Usage
//
// ```
// tabulate(C, FX, S, r0, r1, x).(val|lin|cub) : _
// ```
//
// * `C`: whether to dynamically force the `x` value to the range [r0, r1]: 1 forces the check, 0 deactivates it (constant numerical expression)
// * `FX`: unary function Y=F(X) with one output (scalar function of one variable) 
// * `S`: size of the table in samples (constant numerical expression)
// * `r0`: minimum value of argument x
// * `r1`: maximum value of argument x
//
// ```
// tabulate(C, FX, S, r0, r1, x).val uses the value in the table closest to x
// ```
//
// ```
// tabulate(C, FX, S, r0, r1, x).lin evaluates at x using linear interpolation between the closest stored values
// ```
//
// ```
// tabulate(C, FX, S, r0, r1, x).cub evaluates at x using cubic interpolation between the closest stored values
// ```
//
// #### Example test program
//
// ```
// midikey2hz(mk) = ba.tabulate(1, ba.midikey2hz, 512, 0, 127, mk).lin;
// process = midikey2hz(ba.time), ba.midikey2hz(ba.time);
// ```
//
//--------------------------------------------
tabulate(C, FX, S, r0, r1, x) = environment {

    // Maximum index to access
    mid = S-1;

    // Create the table
    wf = r0 + float(rid(ba.time, 1))*(r1-r0)/float(mid) : FX;

    // Prepare the 'float' table read index
    id = (x-r0)/(r1-r0)*mid;

    // Limit the table read index in [0, mid] if C = 1
    rid(x, 0) = x;
    rid(x, 1) = max(0, min(x, mid));

    // Tabulate an unary 'FX' function on a range [r0, r1]
    val = y0 with { y0 = rdtable(S, wf, rid(int(id+0.5), C)); };

    // Tabulate an unary 'FX' function over the range [r0, r1] with linear interpolation
    lin = it.interpolate_linear(d,y0,y1)
    with {
        x0 = int(id);
        x1 = x0+1;
        d  = id-x0;
        y0 = rdtable(S, wf, rid(x0, C));
        y1 = rdtable(S, wf, rid(x1, C));
    };

    // Tabulate an unary 'FX' function over the range [r0, r1] with cubic interpolation
    cub = it.interpolate_cubic(d,y0,y1,y2,y3)
    with {
        x0 = x1-1;
        x1 = int(id);
        x2 = x1+1;
        x3 = x2+1;
        d  = id-x1;
        y0 = rdtable(S, wf, rid(x0, C));
        y1 = rdtable(S, wf, rid(x1, C));
        y2 = rdtable(S, wf, rid(x2, C));
        y3 = rdtable(S, wf, rid(x3, C));
    };
};

declare tabulate author "Stephane Letz";


//-------`(ba.)tabulate_chebychev`----------
// Tabulate a 1D function over the range [r0, r1] for access via Chebyshev polynomial approximation.
// In contrast to `(ba.)tabulate`, which interpolates only between tabulated samples, `(ba.)tabulate_chebychev`
// stores coefficients of Chebyshev polynomials that are evaluated to provide better approximations in many cases.
// Two new arguments controlling this are NX, the number of segments into which [r0, r1] is divided, and CD,
// the maximum Chebyshev polynomial degree to use for each segment. A `rdtable` of size NX*(CD+1) is internally used.
//
// Note that processing `r1` the last point in the interval is not safe. So either be sure the input stays in [r0, r1[ 
// or use `C = 1`.
//
// #### Usage
//
// ```
// _ : tabulate_chebychev(C, FX, NX, CD, r0, r1) : _
// ```
//
// * `C`: whether to dynamically force the value to the range [r0, r1]: 1 forces the check, 0 deactivates it (constant numerical expression)
// * `FX`: unary function Y=F(X) with one output (scalar function of one variable)
// * `NX`: number of segments for uniformly partitioning [r0, r1] (constant numerical expression)
// * `CD`: maximum polynomial degree for each Chebyshev polynomial (constant numerical expression)
// * `r0`: minimum value of argument x
// * `r1`: maximum value of argument x
//
// #### Example test program
//
// ```
// midikey2hz_chebychev(mk) = ba.tabulate_chebychev(1, ba.midikey2hz, 100, 4, 0, 127, mk);
// process = midikey2hz_chebychev(ba.time), ba.midikey2hz(ba.time);
// ```
//
//--------------------------------------------
tabulate_chebychev(C, FX, NX, CD, r0, r1, x) = y with {
    ck(0) = _;
    ck(1) = max(0) : min(NX-1);

    // number of chebyshev coefficients
    NC = CD + 1;
    // length of the segments
    DX = (r1 - r0) / NX;
    // number of segment 'x' falls in
    nx = (x  - r0) / DX : int : ck(C);
    // center of n's segment
    xc(n) = r0 + DX * (n + 1/2);
    // so ch(0) .. ch(NC) are the coeffs we use for approximation
    // on nx's segment
    ch(i) = chtab(NC * nx + i);

    // map the input in segment [nx*DX, (nx+1)*DX] to [-1,1]
    y = (x - xc(nx)) * 2/DX <: sum(i, NC, ch(i) * ma.chebychev(i));

    // map [-1,1] to the segment [nx*DX, (nx+1)*DX] so mapfx(nx)
    // is simply the "renormalized" FX defined on [-1,1]
    mapfx(nx, x) = FX(xc(nx) + DX/2 * x);

    // calculate the nc's chebyshev coefficient we use on nx's segment
    gench(nx, nc) = (1+(nc!=0))/NC * sum(k,NC,
        mapfx(nx, cos(ma.PI*(k+1/2)/NC)) * cos(ma.PI*nc*(k+1/2)/NC));

    // record gench(nx, nc) in rdtable() to avoid the run-time calculations
    chtab = rdtable(NX*NC, (ba.time <: int(/(NC)), %(NC) : gench));
};

declare tabulate_chebychev author "Oleg Nesterov";
declare tabulate_chebychev copyright "Copyright (C) 2022 Oleg Nesterov <oleg@redhat.com>";
declare tabulate_chebychev license "MIT-style STK-4.3 license";


//-------`(ba.)tabulateNd`----------
// Tabulate an nD function for access via nearest-value or linear or cubic interpolation. In other words, the tabulated function can be evaluated using interpolation of order 0 (none), 1 (linear), or 3 (cubic).  
//
// The table size and parameter range of each dimension can and must be separately specified. You can use it anywhere you have an expensive function with multiple parameters with known ranges. You could use it to build a wavetable synth, for example.
//
// The number of dimensions is deduced from the number of parameters you give, see below.
//
// Note that processing the last point in each interval is not safe. So either be sure the inputs stay in their respective ranges, or use `C = 1`. Similarly for the first point when doing cubic interpolation.
//
// #### Usage
//
// ```
// tabulateNd(C, function, (parameters) ).(val|lin|cub) : _
// ```
//
// * `C`: whether to dynamically force the parameter values for each dimension to the ranges specified in parameters: 1 forces the check, 0 deactivates it (constant numerical expression)
// * `function`: the function we want to tabulate. Can have any number of inputs, but needs to have just one output.
// * `(parameters)`: sizes, ranges and read values. Note: these need to be in brackets, to make them one entity.  
//
//   If N is the number of dimensions, we need:
// 
//   * N times `S`: number of values to store for this dimension (constant numerical expression)
//   * N times `r0`: minimum value of this dimension
//   * N times `r1`: maximum value of this dimension
//   * N times `x`: read value of this dimension
//
// By providing these parameters, you indirectly specify the number of dimensions; it's the number of parameters divided by 4.
//
// The user facing functions are:
// ```
// tabulateNd(C, function, S, parameters).val
// ```
//  - Uses the value in the table closest to x.
// ```
// tabulateNd(C, function, S, parameters).lin
// ```
//  - Evaluates at x using linear interpolation between the closest stored values.
// ```
// tabulateNd(C, function, S, parameters).cub
// ```
//  - Evaluates at x using cubic interpolation between the closest stored values.
//
//
// #### Example test program
//
// ```faust
// powSin(x,y) = sin(pow(x,y)); // The function we want to tabulate
// powSinTable(x,y) = ba.tabulateNd(1, powSin, (sizeX,sizeY, rx0,ry0, rx1,ry1, x,y) ).lin;
// sizeX = 512; // table size of the first parameter
// sizeY = 512; // table size of the second parameter
// rx0 = 2; // start of the range of the first parameter
// ry0 = 2; // start of the range of the second parameter
// rx1 = 10; // end of the range of the first parameter
// ry1 = 10; // end of the range of the second parameter
// x = hslider("x", rx0, rx0, rx1, 0.001):si.smoo;
// y = hslider("y", ry0, ry0, ry1, 0.001):si.smoo;
// process = powSinTable(x,y), powSin(x,y);
// ```
//
// #### Working principle
//
// The ``.val`` function just outputs the closest stored value.
// The ``.lin`` and ``.cub`` functions interpolate in N dimensions.
//
// ##### Multi dimensional interpolation
//
// To understand what it means to interpolate in N dimensions, here's a quick reminder on the general principle of 2D linear interpolation:
//
// * We have a grid of values, and we want to find the value at a point (x, y) within this grid.  
// * We first find the four closest points (A, B, C, D) in the grid surrounding the point (x, y).  
// 
// Then, we perform linear interpolation in the x-direction between points A and B, and between points C and D. This gives us two new points E and F. Finally, we perform linear interpolation in the y-direction between points E and F to get our value.
//
// To implement this in Faust, we need N sequential groups of interpolators, where N is the number of dimensions.  
// Each group feeds into the next, with the last "group" being a single interpolator, and the group before it containing one interpolator for each input of the group it's feeding.
//
// Some examples:
//
// * Our 2D linear example has two interpolators feeding into one.
// * A 3D linear interpolator has four interpolators feeding into two, feeding into one.
// * A 2D cubic interpolater has four interpolators feeding into one.
// * A 3D cubic interpolator has sixteen interpolators feeding into four, feeding into one.
//
// To understand which values we need to look up, let's consider the 2D linear example again.
// The four values going into the first group represent the four closest points (A, B, C, D) mentioned above.
//
// 1) The first interpolator gets:
//
// * The closest value that is stored (A)
// * The next value in the x dimension, keeping y fixed (B)
//
// 2) The second interpolator gets:
//
// * One step over in the y dimension, keeping x fixed (C)
// * One step over in both the x dimension and the y dimension (D)
//
// The outputs of these two interpolators are points E and F.
// In other words: the interpolated x values and, respectively, the following y values:
//
// * The closest stored value of the y dimension
// * One step forward in the y dimension
//
// The last interpolator takes these two values and interpolates them in the y dimension.
//
// To generalize for N dimensions and linear interpolation:
//
// * The first group has 2^(n-1) parallel interpolators interpolating in the first dimension.
// * The second group has 2^(n-2) parallel interpolators interpolating in the second dimension.
// * The process continues until the n-th group, which has a single interpolator interpolating in the n-th dimension.
//
// The same principle applies to the cubic interpolation in nD. The only difference is that there would be 4^(n-1) parallel interpolators in the first group, compared to 2^(n-1) for linear interpolation.
//
// This is what the ``mixers`` function does.
//
// Besides the values, each interpolator also needs to know the weight of each value in it's output.  
// Let's call this `d`, like in ``ba.interpolate``. It is the same for each group of interpolators, since it correlates to a dimension.  
// It's value is calculated the similarly to ``ba.interpolate``:
// 
// * First we prepare a "float table read-index" for that dimension (``id`` in ``ba.tabulate``)
// * If the table only had that dimension and it could read a float index, what would it be.
// * Then we ``int`` the float index to get the value we have stored that is closest to, but lower than the input value; the actual index for that dimension.
// Our ``d`` is the difference between the float index and the actual index.
//
// The ``ids`` function calculates the ``id`` for each dimension and inside the ``mixer`` function they get turned into ``d``s.
//
// ##### Storage method
//
// The elephant in the room is: how do we get these indexes? For that we need to know how the values are stored.
// We use one big table to store everything.
//
// To understand the concept, let's look at the 2D example again, and then we'll extend it to 3d and the general nD case.
//
// Let's say we have a 2D table with dimensions A and B where:
// A has 3 values between 0 and 5 and B has 4 values between 0 and 1.
// The 1D array representation of this 2D table will have a size of 3 * 4 = 12.
//
// The values are stored in the following way:
//
// * First 3 values: A is 0, then 3, then 5 while B is at 0.
// * Next 3 values: A changes from 0 to 5 while B is at 1/3.
// * Next 3 values: A changes from 0 to 5 while B is at 2/3.
// * Last 3 values: A changes from 0 to 5 while B is at 1.
//
// For the 3D example, let's extend the 2D example with an additional dimension C having 2 values between 0 and 2.
// The total size will be 3 * 4 * 2 = 24.
//
// The values are stored like so:
//
// * First 3 values: A changes from 0 to 5, B is at 0, and C is at 0.
// * Next 3 values: A changes from 0 to 5, B is at 1/3, and C is at 0.
// * Next 3 values: A changes from 0 to 5, B is at 2/3, and C is at 0.
// * Next 3 values: A changes from 0 to 5, B is at 1, and C is at 0.
//
// The last 12 values are the same as the first 12, but with C at 2.
//
// For the general n-dimensional case, we iterate through all dimensions, changing the values of the innermost dimension first, then moving towards the outer dimensions.
//
// ##### Read indexes
//
// To get the float read index (``id``) corresponding to a particular dimension, we scale the function input value to be between 0 and 1, and multiply it by the size of that dimension minus one.
//
// To understand how we get the ``readIndex``for ``.val``, let's work trough how we'd do it in our 2D linear example.  
// For simplicity's sake, the ranges of the inputs to our ``function`` are both 0 to 1.  
// Say we wanted to read the value closest to ``x=0.5`` and ``y=0``, so the ``id`` of ``x`` is ``1`` (the second value) and the ``id`` of ``y`` is 0 (first value). In this case, the read index is just the ``id`` of ``x``, rounded to the nearest integer, just like in ``ba.tabulate``.
// 
// If we want to read the value belonging to ``x=0.5`` and ``y=2/3``, things get more complicated. The ``id`` for ``y`` is now ``2``, the third value. For each step in the ``y`` direction, we need to increase the index by ``3``, the number of values that are stored for ``x``. So the influence of the ``y`` is:  the size of ``x`` times the rounded ``id`` of ``y``. The final read index is the rounded ``id`` of ``x`` plus the influence of ``y``.
//
// For the general nD case, we need to do the same operation N times, each feeding into the next. This operation is the ``riN`` function. We take four parameters: the size of the dimension before it ``prevSize``, the index of the previous dimension ``prevIX``, the current size ``sizeX`` and the current id ``idX``.  ``riN`` has 2 outputs, the size, for feeding into the next dimension's ``prevSize``, and the read index feeding into the next dimension's ``prevIX``.  
// The size is the ``sizeX`` times ``prevSize``. The read index is the rounded ``idX`` times ``prevSize`` added to the ``prevIX``. Our final ``readIndex`` is the read index output of the last dimension.
//
// To get the read values for the  interpolators need a pattern of offsets in each dimension, since we are looking for the read indexes surrounding the point of interest. These offsets are best explained by looking at the code of ``tabulate2d``, the hardcoded 2D version:
//
// ```faust
// tabulate2d(C,function, sizeX,sizeY, rx0,ry0, rx1,ry1, x,y) =
//   environment {
//     size = sizeX*sizeY;
//     // Maximum X index to access
//     midX = sizeX-1;
//     // Maximum Y index to access
//     midY = sizeY-1;
//     // Maximum total index to access
//     mid = size-1;
//     // Create the table
//     wf = function(wfX,wfY);
//     // Prepare the 'float' table read index for X
//     idX = (x-rx0)/(rx1-rx0)*midX;
//     // Prepare the 'float' table read index for Y
//     idY = ((y-ry0)/(ry1-ry0))*midY;
//     // table creation X:
//     wfX =
//       rx0+float(ba.time%sizeX)*(rx1-rx0)
//       /float(midX);
//     // table creation Y:
//     wfY =
//       ry0+
//       ((float(ba.time-(ba.time%sizeX))
//         /float(sizeX))
//        *(ry1-ry0))
//       /float(midY);
//
//     // Limit the table read index in [0, mid] if C = 1
//     rid(x,mid, 0) = x;
//     rid(x,mid, 1) = max(0, min(x, mid));
//
//     // Tabulate a binary 'FX' function on a range [rx0, rx1] [ry0, ry1]
//     val(x,y) =
//       rdtable(size, wf, readIndex);
//     readIndex =
//       rid(
//         rid(int(idX+0.5),midX, C)
//         +yOffset
//       , mid, C);
//     yOffset = sizeX*rid(int(idY),midY,C);
//
//     // Tabulate a binary 'FX' function over the range [rx0, rx1] [ry0, ry1] with linear interpolation
//     lin =
//       it.interpolate_linear(
//         dy
//       , it.interpolate_linear(dx,v0,v1)
//       , it.interpolate_linear(dx,v2,v3))
//     with {
//       i0 = rid(int(idX), midX, C)+yOffset;
//       i1 = i0+1;
//       i2 = i0+sizeX;
//       i3 = i1+sizeX;
//       dx  = idX-int(idX);
//       dy  = idY-int(idY);
//       v0 = rdtable(size, wf, rid(i0, mid, C));
//       v1 = rdtable(size, wf, rid(i1, mid, C));
//       v2 = rdtable(size, wf, rid(i2, mid, C));
//       v3 = rdtable(size, wf, rid(i3, mid, C));
//     };
// 
//     // Tabulate a binary 'FX' function over the range [rx0, rx1] [ry0, ry1] with cubic interpolation
//     cub =
//       it.interpolate_cubic(
//         dy
//       , it.interpolate_cubic(dx,v0,v1,v2,v3)
//       , it.interpolate_cubic(dx,v4,v5,v6,v7)
//       , it.interpolate_cubic(dx,v8,v9,v10,v11)
//       , it.interpolate_cubic(dx,v12,v13,v14,v15)
//       )
//     with {
//       i0  = i4-sizeX;
//       i1  = i5-sizeX;
//       i2  = i6-sizeX;
//       i3  = i7-sizeX;
//
//       i4  = i5-1;
//       i5  = rid(int(idX), midX, C)+yOffset;
//       i6  = i5+1;
//       i7  = i6+1;
//
//       i8  = i4+sizeX;
//       i9  = i5+sizeX;
//       i10 = i6+sizeX;
//       i11 = i7+sizeX;
//
//       i12 = i4+(2*sizeX);
//       i13 = i5+(2*sizeX);
//       i14 = i6+(2*sizeX);
//       i15 = i7+(2*sizeX);
//
//       dx  = idX-int(idX);
//       dy  = idY-int(idY);
//       v0  = rdtable(size, wf, rid(i0 , mid, C));
//       v1  = rdtable(size, wf, rid(i1 , mid, C));
//       v2  = rdtable(size, wf, rid(i2 , mid, C));
//       v3  = rdtable(size, wf, rid(i3 , mid, C));
//       v4  = rdtable(size, wf, rid(i4 , mid, C));
//       v5  = rdtable(size, wf, rid(i5 , mid, C));
//       v6  = rdtable(size, wf, rid(i6 , mid, C));
//       v7  = rdtable(size, wf, rid(i7 , mid, C));
//       v8  = rdtable(size, wf, rid(i8 , mid, C));
//       v9  = rdtable(size, wf, rid(i9 , mid, C));
//       v10 = rdtable(size, wf, rid(i10, mid, C));
//       v11 = rdtable(size, wf, rid(i11, mid, C));
//       v12 = rdtable(size, wf, rid(i12, mid, C));
//       v13 = rdtable(size, wf, rid(i13, mid, C));
//       v14 = rdtable(size, wf, rid(i14, mid, C));
//       v15 = rdtable(size, wf, rid(i15, mid, C));
//     };
//   };
// ```
//
// In the interest of brevity, we'll stop explaining here. If you have any more questions, feel free to open an issue on [faustlibraries](https://github.com/grame-cncm/faustlibraries) and tag @magnetophon.
//
//--------------------------------------------
tabulateNd(C,function,parameters) =
    environment {
        val =
            parameters
            // our sizes can be int, should be faster to calculate
            : (par(i, N, int),si.bus(N*3))
              // table size, waveform, read index
              <: (tableSize,wf,readIndex(idsGetRoundedNotFloored))
            : rdtable
        with {
        // for val we want the rounded version of id, the interpolators need the rounded down version
        // see: https://github.com/grame-cncm/faustlibraries/pull/152
        idsGetRoundedNotFloored =
            // table sizes and ids for each dimension
            sizesIds:
            (bs,par(i, N, _+0.5));
        };

        lin =
            parameters
            // our sizes can be int, should be faster to calculate
            : (par(i, N, int),si.bus(N*3))
              // the mixers need the float indexes and the values read from the tables
              <: (ids,tables(nrReadIndexes,readIndexes))
                 // the actual interpolation
            :mixers(0,nrReadIndexes)
        with {
            // the read indexes to form the closest points in the grid surrounding the point of interest.
            readIndexes =
                // for cleaner block diagrams
                si.bus(nParams) <:
                // the closest stored value at or below the point of interest
                ((readIndex(sizesIds) <:si.bus(nrReadIndexes))
                 // the offsets to get to the points around it
                , offsets)
                // add them up
                : ro.interleave(nrReadIndexes,2) : par(i, nrReadIndexes, +) ;
            offsets =
                // the size of a step for each dimension
                stepSizes
                // for each read index
                <: par(i, nrReadIndexes,
                       // sum up the N components that determine the final offset
                       par(j, N, switch(i,j)):>_)
            with {
                // the truth table of which step size to use is a binary counting table
                // so we use a variant on a function I wrote for slidingReduce
                switch(i,j) = _*int2bin(j,i,nrReadIndexes);
            };
            // since each interpolator has 2 inputs
            nrReadIndexes = pow(2,N);
        };

        cub =
            // same as .lin so far
            parameters
            : (par(i, N, int),si.bus(N*3))
              <: (ids,tables(nrReadIndexes,readIndexes))
            :mixers(1,nrReadIndexes)
        with {
            readIndexes =
                si.bus(nParams) <:
                ((stepSizes:ro.cross(N))
                , readIndex(sizesIds))
                // calculate the offsets for each interpolation read index and add them to the "base" readindex
                :cubifiers;
            // since each interpolator has 4 inputs
            nrReadIndexes = pow(4,N);
            // calculate the offsets and add them to the "base" readindex
            // we do this by creating a tree, where each branch has 4 sub-branches
            // there are N splits, so we end up with nrReadIndexes "leaves"
            cubifiers =
                // iterate trough the dimensions, each feeding into the next
                seq(i, N,
                    si.bus(N-i-1),cubifier(pow(4,i)));
            // take a step size and ``len`` input indexes and create 4*len output indexes
            // offset the input index(es) with the amount needed for that dimension
            cubifier(len) =
                ((_<:si.bus(len*4))
                , (si.bus(len)<:si.bus(len*4)))
                : ro.interleave(len*4,2)
                :par(i, 4,
                     par(j, len, off(i)+_))
            with {
                // the hardcoded numbers are the same offsets as in ``ba.tabulate``
                // but they get multiplied by the step size of the current dimension
                off(0,stepSize) = -1*stepSize;
                off(1,stepSize) = 0;
                off(2,stepSize) = 1*stepSize;
                off(3,stepSize) = 2*stepSize;
            };
        };

        // how many parameters have we been given?
        nParams = outputs(parameters);
        // the number of dimensions
        N = int(nParams/4);
        // the values we read from the tables
        tables(nrReadIndexes,readIndexes) =
            si.bus(nParams)<:
            // table size, waveform, read index for each table
            ((tableSize<:si.bus(nrReadIndexes))
            , (wf<:si.bus(nrReadIndexes))
            , readIndexes)
            :ro.interleave(nrReadIndexes,3)
             // the actual tables
            :par(i, nrReadIndexes, rdtable);

        // the interpolators
        mixers(linCub,nrReadIndexes)=
            (ro.cross(N),si.bus(nrReadIndexes))
            : seq(i, N, mixer(linCub,i));

        // the interpolator for a single dimension
        // linear version
        mixer(0,i) =
            mixerUniversal(i,2,(_,!,_),it.interpolate_linear) ;
        // cubic version
        mixer(1,i) =
            mixerUniversal(i,4,(_,!,_,!,_,!,_),it.interpolate_cubic) ;

        // i is the current dimension
        // mult is the number of inputs of each interpolator
        // sieve tells us which outputs to let trough and which to block, more on this soon
        mixerUniversal(i,mult,sieve,it) =
            // bypass the weights needed for the next dimension
            si.bus(N-i-1),
            // split our own weight
            // in the end we need one weight per interpolator
            // but since we need to interleave these with nrInterpolators*mult read values, we split it into as many busses as we have read values
            (((_<:si.bus(nrInterpolators(i)*mult))
              // the read values bypass this step
             , (si.bus(nrInterpolators(i)*mult)))
             // interleave the weights with the read values
             : ro.interleave(nrInterpolators(i)*mult,2)
               // the actual interpolators
             : par(i, nrInterpolators(i),
                   // take the id and turn it into the weight for this dimension
                   ((_<:(_-int(_)))
                    // throw away the extra weights we created for the interleave
                   ,sieve)
                   // a single interpolator
                   :it))
        with {
            // the number of interpolators for this dimension
            nrInterpolators(i) = pow(mult,N-i-1);
        };

        // total size of the table: s(0) * s(1)  ...  * s(N-2) * s(N-1)
        // N in, 1 out
        size(1) = _;
        size(N) = _*size(N-1);
        tableSize = size(N),par(i, N*3, !);
        // the size of a step for each dimension.
        // the first is 1, the second is the size of the first dimension times one
        // the third is the second number, times the size of the second dimension
        stepSizes =
            (int(1),si.bus(N-1),par(i, 3*N+1, !))
            : seq(i, N-1,
                  ((si.bus(i),(_<:(_,_)), si.bus(N-i-1))
                   :(si.bus(i+1),*,si.bus(N-i-2))));
        // Prepare the 'float' table read index for one parameter
        id(sizeX,r0,r1,x) = (x-r0)/(r1-r0)*(sizeX-1);
        // Prepare the 'float' table read index for all parameters
        ids =
            ro.interleave(N,4)
            : par(i, N, id) ;

        // one waveform parameter write value:
        wfp(prevSize,sizeX,r0,r1) =
            r0+
            ((float(int(ba.time%(prevSize*sizeX)/prevSize))*(r1-r0))
            /float(sizeX-1))
           ,(prevSize*sizeX);

        // all waveform parameters write values:
        wfps =
            ro.interleave(N,3)
            : (1,si.bus(3*N))
            : seq(i, N, si.bus(i),wfp, si.bus(3*N-(3*(i+1))))
            : (si.bus(N),!);

        // Create the table
        wf = (wfps,par(i, N, !)):function;

        // Limit the table read index in [0, mid] if C = 1
        // we do this both for the total table and per dimension
        rid(x,mid, 0) = int(x);
        rid(x,mid, 1) = max(int(0), min(int(x), mid));

        // for ``.val`` this is the stored value closest to the point of interest
        // for  ``.lin`` and ``cub`` it's the closest stored value at or below the point of interest
        readIndex(sizesIds) =
            // table sizes and ids for each dimension
            sizesIds
            // get the raw read index
            : ri
              // limit it and make it an int
            : riPost ;
        // helper function to route the arguments of rid
        riPost(size,ri) =
            rid(ri,size-int(1),C);
        // the raw read index
        ri =
            // interleave the sizes and the ids
            ro.interleave(N,2)
            // the first iteration gets:
            // 1 as the size of the previous dimension and 0 as the read index of the previous dimension
            : (1,0
                 // pass trough the sizes and ids
               ,si.bus(2*N))
              // for each dimension
            : seq(i, N,
                  // get the step size and the partial index for the next dimension
                  riN
                  // pass trough the sizes and ids for the next dimensions
                  , si.bus(2*(N-i-1))) ;

        // get the step size and the partial index for the next dimension
        riN(prevSize,prevIX,sizeX,idX) =
            // each step size is the previous step size times the current
            (prevSize*sizeX)
            // but the step we actually take in this dimension is the stepsize we calculated in the previous dimension
          , ( (prevSize*
               // round down and limit the id for this dimension to get the number of steps
               rid(int(idX),(sizeX-int(1)),C))
              // add the result of that to the one we calculated in the previous dimension
              +prevIX) ;

        // table sizes and ids for each dimension
        // just a helper function to simplify routing
        sizesIds =
            (
                // from sizes
                ( bs<:si.bus(N*2) )
                // from r0s,r1s,xs
              , si.bus(N*3)
            ) :
            // from sizes
            (si.bus(N)
             // takes (midX,r0,r1,x)
            ,ids);

        // the value of a single bit when binary counting
        int2bin(i,n,maxN) = int(int((n)/(1<<i))%int(2));
        // shortcut
        bs = si.bus(N);
};

declare tabulateNd author "Bart Brouns";
declare tabulateNd "Copyright (C) 2023 Bart Brouns <bart@magnetophon.nl>";
declare tabulateNd license "AGPL-3.0";

//============================Selectors (Conditions)======================================
//========================================================================================

//-----------------------------`(ba.)if`-----------------------------------
// if-then-else implemented with a select2. WARNING: since `select2` is strict (always evaluating both branches),
// the resulting if does not have the usual "lazy" semantic of the C if form, and thus cannot be used to
// protect against forbidden computations like division-by-zero for instance.
//
// #### Usage
//
// *   `if(cond, then, else) : _`
//
// Where:
//
// * `cond`: condition
// * `then`: signal selected while cond is true
// * `else`: signal selected while cond is false
//-----------------------------------------------------------------------------
if(cond,then,else) = select2(cond,else,then);
// TODO: perhaps it would make more sense to have an if(a,b) and an ifelse(a,b,c)?

//-----------------------------`(ba.)ifNc`--------------------------------------
// if-then-elseif-then-...elsif-then-else implemented on top of `ba.if`.
//
// #### Usage
//
// ```
//    ifNc((cond1,then1, cond2,then2, ... condN,thenN, else)) : _
// or
//    ifNc(Nc, cond1,then1, cond2,then2, ... condN,thenN, else) : _
// or
//    cond1,then1, cond2,then2, ... condN,thenN, else : ifNc(Nc) : _
// ```
//
// Where:
//
// * `Nc` : number of branches/conditions (constant numerical expression)
// * `condX`: condition
// * `thenX`: signal selected if condX is the 1st true condition
// * `else`: signal selected if all the cond1-condN conditions are false
//
// #### Example test program
//
// ```
//    process(x,y) = ifNc((x<y,-1, x>y,+1, 0));
// or
//    process(x,y) = ifNc(2, x<y,-1, x>y,+1, 0);
// or
//    process(x,y) = x<y,-1, x>y,+1, 0 : ifNc(2);
// ```
//
// outputs `-1` if `x<y`, `+1` if `x>y`, `0` otherwise.
//-----------------------------------------------------------------------------
ifNc((c,r)) = c,r : ifNc(outputs(r)/2);
ifNc(1) = if;
ifNc(n) = _,_,ifNc(n-1) : ifNc(1);

declare ifNc author "Oleg Nesterov";
declare ifNc copyright "Copyright (C) 2023 Oleg Nesterov <oleg@redhat.com>";
declare ifNc license "MIT-style STK-4.3 license";

//-----------------------------`(ba.)ifNcNo`-------------------------------------
// `ifNcNo(Nc,No)` is similar to `ifNc(Nc)` above but then/else branches have `No` outputs.
//
// #### Usage
//
// ```
//    ifNcNo(Nc,No, cond1,then1, cond2,then2, ... condN,thenN, else) : sig.bus(No)
// ```
//
// Where:
//
// * `Nc` : number of branches/conditions (constant numerical expression)
// * `No` : number of outputs (constant numerical expression)
// * `condX`: condition
// * `thenX`: list of No signals selected if condX is the 1st true condition
// * `else`: list of No signals selected if all the cond1-condN conditions are false
//
// #### Example test program
//
// ```
//    process(x) = ifNcNo(2,3, x<0, -1,-1,-1, x>0, 1,1,1, 0,0,0);
// ```
//
// outputs `-1,-1,-1` if `x<0`, `1,1,1` if `x>0`, `0,0,0` otherwise.
//-----------------------------------------------------------------------------
ifNcNo(Nc,No) = si.bus(Ni) <: par(no,No, slice(no) : ifNc(Nc))
with {
    Ni = Nc*No + Nc+No;
    slice(no) = route(Ni,2*Nc+1, (par(nc,Nc, ct(nc,no)),else(no)));
    ct(nc,no) = nc*(No+1)+1,    2*nc+1,    // cond[nc]
                nc*(No+1)+2+no, 2*nc+2;    // then[nc][no]
    else(no)  = Nc*(No+1)+1+no, 2*Nc+1;    // else[no]
};

declare ifNcNo author "Oleg Nesterov";
declare ifNcNo copyright "Copyright (C) 2023 Oleg Nesterov <oleg@redhat.com>";
declare ifNcNo license "MIT-style STK-4.3 license";

//-----------------------------`(ba.)selector`---------------------------------
// Selects the ith input among n at compile time.
//
// #### Usage
//
// ```
// selector(I,N)
// _,_,_,_ : selector(2,4) : _ // selects the 3rd input among 4
// ```
//
// Where:
//
// * `I`: input to select (int, numbered from 0, known at compile time)
// * `N`: number of inputs (int, known at compile time, N > I)
//
// There is also `cselector` for selecting among complex input signals of the form (real,imag).
//
//-----------------------------------------------------------------------------
selector(i,n) = par(j, n, S(i, j))    with { S(i,i) = _; S(i,j) = !; };
cselector(i,n) = par(j, n, S(i, j))   with { S(i,i) = (_,_); S(i,j) = (!,!); }; // for complex numbers


//--------------------`(ba.)select2stereo`--------------------
// Select between 2 stereo signals.
//
// #### Usage
//
// ```
// _,_,_,_ : select2stereo(bpc) : _,_
// ```
//
// Where:
//
// * `bpc`: the selector switch (0/1)
//------------------------------------------------------------
select2stereo(bpc) = ro.cross2 : select2(bpc), select2(bpc) : _,_;


//-----------------------------`(ba.)selectn`---------------------------------
// Selects the ith input among N at run time.
//
// #### Usage
//
// ```
// selectn(N,i)
// _,_,_,_ : selectn(4,2) : _ // selects the 3rd input among 4
// ```
//
// Where:
//
// * `N`: number of inputs (int, known at compile time, N > 0)
// * `i`: input to select (int, numbered from 0)
//
// #### Example test program
//
// ```
// N = 64;
// process = par(n, N, (par(i,N,i) : selectn(N,n)));
// ```
//-----------------------------------------------------------------------------
selectn(N,i) = selectnX(N,i,selector)
with {
    selector(i,j,x,y) = select2((i >= j), x, y);
};


// The generic version with a 'sel' function to be applied on:
// - the channel index as a (possibly) fractional value
// - the next channel index as an integer value
// - the 2 signals to be selected between

selectnX(N,i,sel) = S(N,0)
with {
    S(1,offset) = _;
    S(n,offset) = S(left, offset), S(right, offset+left) : sel(i, offset+left)
    with {
        right = int(n/2);
        left  = n-right;
    };
};


//----------------------`(ba.)selectbus`-----------------------------------------
// Select a bus among `NUM_BUSES` buses, where each bus has `BUS_SIZE` outputs.
// The order of the signal inputs should be the signals of the first bus, the
// signals of the second bus, and so on.
//
// #### Usage
//
// ```
// process = si.bus(BUS_SIZE*NUM_BUSES) : selectbus(BUS_SIZE, NUM_BUSES, id) : si.bus(BUS_SIZE);
// ```
//
// Where:
//
// * `BUS_SIZE`: number of outputs from each bus (int, known at compile time).
// * `NUM_BUSES`: number of buses (int, known at compile time).
// * `id`: index of the bus to select (int, `0<=id<NUM_BUSES`)
//
//-------------------------------------------------------------------
selectbus(BUS_SIZE, NUM_BUSES, id) = ro.interleave(BUS_SIZE, NUM_BUSES) : par(j, BUS_SIZE, selectn(NUM_BUSES, id));
declare selectbus author "David Braun";
declare selectbus license "MIT";


//----------------------`(ba.)selectxbus`-----------------------------------------
// Like `ba.selectbus`, but with a cross-fade when selecting the bus using the same 
// technique than `ba.selectmulti`.
//
// #### Usage
//
// ```
// process = si.bus(BUS_SIZE*NUM_BUSES) : selectbus(BUS_SIZE, NUM_BUSES, FADE, id) : si.bus(BUS_SIZE);
// ```
//
// Where:
//
// * `BUS_SIZE`: number of outputs from each bus (int, known at compile time).
// * `NUM_BUSES`: number of buses (int, known at compile time).
// * `fade`: number of samples for the crossfade.
// * `id`: index of the bus to select (int, `0<=id<NUM_BUSES`)
//
//-------------------------------------------------------------------
selectxbus(BUS_SIZE, NUM_BUSES, fade, id) =
    ro.interleave(BUS_SIZE, NUM_BUSES):
    par(i, BUS_SIZE, selectnX(NUM_BUSES, id, xfade))
    with {
        xfade(i, j, x, y) = x*(1-xb) + y*xb with {
            xb = ramp(fade, (i >= j));
        };
    };
declare selectxbus author "Marc Lavallée";
declare selectxbus license "MIT";


//-----------------------------`(ba.)selectmulti`---------------------------------
// Selects the ith circuit among N at run time (all should have the same number of inputs and outputs)
// with a crossfade.
//
// #### Usage
//
// ```
// selectmulti(n,lgen,id)
// ```
//
// Where:
//
// * `n`: crossfade in samples
// * `lgen`: list of circuits
// * `id`: circuit to select (int, numbered from 0)
//
// #### Example test program
//
// ```
// process = selectmulti(ma.SR/10, ((3,9),(2,8),(5,7)), nentry("choice", 0, 0, 2, 1));
// process = selectmulti(ma.SR/10, ((_*3,_*9),(_*2,_*8),(_*5,_*7)), nentry("choice", 0, 0, 2, 1));
// ```
//-----------------------------------------------------------------------------
selectmulti(n, lgen, id) = selectmultiX(ins, lgen, id)
with {
    selectmultiX(0, lgen, id) = selector;                    // No inputs
    selectmultiX(N, lgen, id) = par(i, ins, _) <: selector;  // General case

    selector = lgen : ro.interleave(outs, N) : par(i, outs, selectnX(N, id, xfade))
    with {
        // crossfade of 'n' samples between 'x' and 'y' channels when the channel index changes
        xfade(i, j, x, y) = x*(1-xb) + y*xb with { xb = ramp(n, (i >= j)); };
    };

    outs = outputs(take(1, lgen));  // Number of outputs of the first circuit (all should have the same value)
    ins = inputs(take(1, lgen));    // Number of inputs of the first circuit (all should have the same value)
    N = outputs(lgen)/outs;         // Number of items in the list
};


//-----------------------------`(ba.)selectoutn`---------------------------------
// Route input to the output among N at run time.
//
// #### Usage
//
// ```
// _ : selectoutn(N, i) : si.bus(N)
// ```
//
// Where:
//
// * `N`: number of outputs (int, known at compile time, N > 0)
// * `i`: output number to route to (int, numbered from 0) (i.e. slider)
//
// #### Example test program
//
// ```
// process = 1 : selectoutn(3, sel) : par(i, 3, vbargraph("v.bargraph %i", 0, 1));
// sel = hslider("volume", 0, 0, 2, 1) : int;
// ```
//--------------------------------------------------------------------------
declare selectoutn author "Vince";

selectoutn(N, s) = _ <: par(i, N, *(s==i));


//=====================================Other==============================================
//========================================================================================

//----------------------------`(ba.)latch`--------------------------------
// Latch input on positive-going transition of trig: "records" the input when trig
// switches from 0 to 1, outputs a frozen values everytime else.
//
// #### Usage
//
// ```
// _ : latch(trig) : _
// ```
//
// Where:
//
// * `trig`: hold trigger (0 for hold, 1 for bypass)
//------------------------------------------------------------
latch(trig, x) = x * s : + ~ *(1-s) with { s = (trig' <= 0) & (trig > 0); };


//--------------------------`(ba.)sAndH`-------------------------------
// Sample And Hold: "records" the input when trig is 1, outputs a frozen value when trig is 0.
// `sAndH` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : sAndH(trig) : _
// ```
//
// Where:
//
// * `trig`: hold trigger (0 for hold, 1 for bypass)
//----------------------------------------------------------------
declare sAndH author "Romain Michon";

sAndH(trig) = select2(trig) ~ _;

//--------------------------`(ba.)tAndH`-------------------------------
// Test And Hold: "records" the input when pred(input) is true, outputs a frozen value otherwise.
//
// #### Usage
//
// ```
// _ : tAndH(pred) : _
// ```
//
// Where:
//
// * `pred`: predicate to test the input
//----------------------------------------------------------------
tAndH(pred) = _ <: pred,_ : sAndH;

//--------------------------`(ba.)downSample`-------------------------------
// Down sample a signal. WARNING: this function doesn't change the
// rate of a signal, it just holds samples...
// `downSample` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : downSample(freq) : _
// ```
//
// Where:
//
// * `freq`: new rate in Hz
//----------------------------------------------------------------
declare downSample author "Romain Michon";

downSample(freq) = sAndH(hold)
with {
    hold = time%int(ma.SR/freq) == 0;
};


//--------------------------`(ba.)downSampleCV`---------------------------
//
// A version of `ba.downSample` where the frequency parameter has
// been replaced by an `amount` parameter that is in the range zero
// to one. WARNING: this function doesn't change the rate of a
// signal, it just holds samples...
//
// #### Usage
// ```
// _ : downSampleCV(amount) : _
// ```
// Where:
//
// * `amount`: The amount of down-sampling to perform [0..1]
//-----------------------------------------------------------------------
downSampleCV(amt) = sAndH(hold)
with {
    // If amt is 0, then freq should be ma.SR/2.
    // If amt is 1, then freq should be 491.
    // 491 is derived from listening to the Serum synthesizer.
    freq = amt : it.remap(0, 1, hz2midikey(ma.SR/2), hz2midikey(491)) : midikey2hz;
    hold = time%int(ma.SR/freq) == 0;
};

declare downSampleCV author "Romain Michon/David Braun";


//------------------`(ba.)peakhold`---------------------------
// Outputs current max value above zero.
//
// #### Usage
//
// ```
// _ : peakhold(mode) : _
// ```
//
// Where:
//
// `mode` means:
//
//  0 - Pass through. A single sample 0 trigger will work as a reset.
//
//  1 - Track and hold max value.
//----------------------------------------------------------------
declare peakhold author "Jonatan Liljedahl, revised by Romain Michon";

peakhold = (*,_ : max) ~ _;

//------------------`(ba.)peakholder`-------------------------------------------
//
// While peak-holder functions are scarcely discussed in the literature
// (please do send me an email if you know otherwise), common sense
// tells that the expected behaviour should be as follows: the absolute
// value of the input signal is compared with the output of the peak-holder;
// if the input is greater or equal to the output, a new peak is detected
// and sent to the output; otherwise, a timer starts and the current peak
// is held for N samples; once the timer is out and no new peaks have been
// detected, the absolute value of the current input becomes the new peak.
//
// #### Usage
//
// ```
// _ : peakholder(holdTime) : _
// ```
//
// Where:
//
// * `holdTime`: hold time in samples
//------------------------------------------------------------------------------
declare peakholder author "Dario Sanfilippo";
declare peakholder copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare peakholder license "MIT-style STK-4.3 license";
peakholder(holdTime, x) = loop ~ si.bus(2) : ! , _
    with {
        loop(timerState, outState) = timer , output
            with {
                isNewPeak = abs(x) >= outState;
                isTimeOut = timerState >= holdTime;
                bypass = isNewPeak | isTimeOut;
                timer = ba.if(bypass, 0, timerState + 1);
                output = ba.if(bypass, abs(x), outState);
            };
    };


/*
// Alternate version with branchless code
//----------------------------------------
peakholder(holdTime, x) = loop ~ si.bus(2) : ! , _
    with {
        loop(timerState, outState) = timer , output
            with {
                isNewPeak = abs(x) >= outState;
                isTimeOut = timerState >= holdTime;
                bypass = isNewPeak | isTimeOut;
                timer = (1 - bypass) * (timerState + 1);
                output = bypass * (abs(x) - outState) + outState;
            };
    };
*/

/*
// The function below is kept for back-compatibility in case any user relies
// on it for their software. However, the function behaves differently than
// expected: currently, the timer is not reset when a new peak is detected.
//------------------------------------------------------------------------------
declare peakholder author "Jonatan Liljedahl";
peakholder(n) = peakhold2 ~ reset : (!,_) with {
    reset = sweep(n) > 0;
    // first out is gate that is 1 while holding last peak
    peakhold2 = _,abs <: peakhold,!,_ <: >=,_,!;
};
*/

//--------------------------`(ba.)kr2ar`---------------------------
// Force a control rate signal to be used as an audio rate signal.
//
// #### Usage
//
// ```
// hslider("freq", 200, 200, 2000, 0.1) : kr2ar;
// ```
//----------------------------------------------------------------
kr2ar = + ~ *(0); 


//--------------------------`(ba.)impulsify`---------------------------
// Turns a signal into an impulse with the value of the current sample
// (0.3,0.2,0.1 becomes 0.3,0.0,0.0). This function is typically used with a
// `button` to turn its output into an impulse. `impulsify` is a standard Faust
// function.
//
// #### Usage
//
// ```
// button("gate") : impulsify;
// ```
//----------------------------------------------------------------
impulsify = _ <: _,mem : - <: >(0)*_;


//-----------------------`(ba.)automat`------------------------------
// Record and replay in a loop the successives values of the input signal.
//
// #### Usage
//
// ```
// hslider(...) : automat(t, size, init) : _
// ```
//
// Where:
//
// * `t`: tempo in BPM
// * `size`: number of items in the loop
// * `init`: init value in the loop
//-----------------------------------------------------------------------
automat(t, size, init, input) = rwtable(size+1, init, windex, input, rindex)
with {
    clock = beat(t);
    rindex = int(clock) : (+ : %(size)) ~ _;    // each clock read the next entry of the table
    windex = if(timeToRenew, rindex, size);     // we ignore input unless it is time to renew
    timeToRenew = int(clock) & (inputHasMoved | (input <= init));
    inputHasMoved = abs(input-input') : countfrom(int(clock)') : >(0);
    countfrom(reset) = (+ : if(reset, 0, _)) ~ _;
};


//-----------------`(ba.)bpf`-------------------
// bpf is an environment (a group of related definitions) that can be used to
// create break-point functions. It contains three functions:
//
// * `start(x,y)` to start a break-point function
// * `end(x,y)` to end a break-point function
// * `point(x,y)` to add intermediate points to a break-point function, using linear interpolation
//
// A minimal break-point function must contain at least a start and an end point:
//
// ```
// f = bpf.start(x0,y0) : bpf.end(x1,y1);
// ```
//
// A more involved break-point function can contains any number of intermediate
// points:
//
// ```
// f = bpf.start(x0,y0) : bpf.point(x1,y1) : bpf.point(x2,y2) : bpf.end(x3,y3);
// ```
//
// In any case the `x_{i}` must be in increasing order (for all `i`, `x_{i} < x_{i+1}`).
// For example the following definition:
//
// ```
// f = bpf.start(x0,y0) : ... : bpf.point(xi,yi) : ... : bpf.end(xn,yn);
// ```
//
// implements a break-point function f such that:
//
// * `f(x) = y_{0}` when `x < x_{0}`
// * `f(x) = y_{n}` when `x > x_{n}`
// * `f(x) = y_{i} + (y_{i+1}-y_{i})*(x-x_{i})/(x_{i+1}-x_{i})` when `x_{i} <= x`
// and `x < x_{i+1}`
//
// In addition to `bpf.point`, there are also `step` and `curve` functions:
//
// * `step(x,y)` to add a flat section
// * `step_end(x,y)` to end with a flat section
// * `curve(B,x,y)` to add a curved section
// * `curve_end(B,x,y)` to end with a curved section
//
// These functions can be combined with the other `bpf` functions.
//
// Here's an example using `bpf.step`:
//
// `f(x) = x : bpf.start(0,0) : bpf.step(.2,.3) : bpf.step(.4,.6) : bpf.step_end(1,1);`
//
// For `x < 0.0`, the output is 0.0.
// For `0.0 <= x < 0.2`, the output is 0.0.
// For `0.2 <= x < 0.4`, the output is 0.3.
// For `0.4 <= x < 1.0`, the output is 0.6.
// For `1.0 <= x`, the output is 1.0
//
// For the `curve` functions, `B` (compile-time constant)
// is a "bias" value strictly greater than zero and less than or equal to 1. When `B` is 0.5, the
// output curve is exactly linear and equivalent to `bpf.point`. When `B` is less than 0.5, the
// output is biased towards the `y` value of the previous breakpoint. When `B` is greater than 0.5,
// the output is biased towards the `y` value of the curve breakpoint. Here's an example:
//
// `f = bpf.start(0,0) : bpf.curve(.15,.5,.5) : bpf.curve_end(.85,1,1);`
//
// In the following example, the output is biased towards zero (the latter y value) instead of
// being a linear ramp from 1 to 0.
//
// `f = bpf.start(0,1) : bpf.curve_end(.9,1,0);`
//
// `bpf` is a standard Faust function.
//--------------------------------------------------------
bpf = environment
{
    // Start a break-point function
    start(x0,y0) = \(x).(x0,y0,x,y0);
    // Add a break-point
    point(x1,y1) = \(x0,y0,x,y).(x1, y1, x, if(x < x0, y, if(x < x1, y0 + (x-x0)*(y1-y0)/(x1-x0), y1)));
    // End a break-point function
    end(x1,y1) = \(x0,y0,x,y).(if(x < x0, y, if(x < x1, y0 + (x-x0)*(y1-y0)/(x1-x0), y1)));

    // Add a stepped "flat" section
    step(x1,y1) = \(x0,y0,x,y).(x1, y1, x, if(x < x0, y, if(x < x1, y0, y1)));
    // End with a stepped "flat" section
    step_end(x1,y1) = \(x0,y0,x,y).(if(x < x0, y, if(x < x1, y0, y1)));

    // Add a curve-point. `B` (compile-time constant) must be in range (0,1]
    curve(B,x1,y1) = \(x0,y0,x,y).(x1, y1, x, if(x < x0, y, if(x < x1, _curve_util(x0,x1,y0,y1,B,x), y1)));
    // End with a curve-point. `B` (compile-time constant) must be in range (0,1]
    curve_end(B,x1,y1) = \(x0,y0,x,y).(if(x < x0, y, if(x < x1, _curve_util(x0,x1,y0,y1,B,x), y1)));

    /*
    In the functions below, we use the equation `y=bias_curve(b,x)`.
    * `b`: bias in range (0,1]. Bias of 0.5 results in `y=x`. Bias above 0.5 pulls y upward. Bias below 0.5 pulls y downward.
    * `x`: input in range [0,1] that needs to be remapped/biased into `y`
    * `y`: `x` after it has been biased. Note that `(y==0 iff x==0) AND (y==1 iff x==1)`.
    */

    _bias_curve(b,x) = (x / ((((1.0/b) - 2.0)*(1.0 - x))+1.0));

    _curve_util(x0,x1,y0,y1,b,x) = y
    with {
        x_norm = (x-x0)/(x1-x0); // x_norm is in [0-1]
        y_norm = _bias_curve(b, x_norm); // y_norm is in [0-1]
        y = (1-y_norm)*y0 + y_norm*y1;  // linear interpolation between y0 and y1
    };
};


//-------------------`(ba.)listInterp`-------------------------
// Linearly interpolates between the elements of a list.
//
// #### Usage
//
// ```
// index = 1.69; // range is 0-4
// process = listInterp((800,400,350,450,325),index);
// ```
//
// Where:
//
// * `index`: the index (float) to interpolate between the different values.
// The range of `index` depends on the size of the list.
//------------------------------------------------------------
declare listInterp author "Romain Michon";

listInterp(v) =
    bpf.start(0,take(1,v)) :
    seq(i,count(v)-2,bpf.point(i+1,take(i+2,v))) :
    bpf.end(count(v)-1,take(count(v),v));


//-------------------`(ba.)bypass1`-------------------------
// Takes a mono input signal, route it to `e` and bypass it if `bpc = 1`.
// When bypassed, `e` is feed with zeros so that its state is cleanup up.
// `bypass1` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : bypass1(bpc,e) : _
// ```
//
// Where:
//
// * `bpc`: bypass switch (0/1)
// * `e`: a mono effect
//------------------------------------------------------------
declare bypass1 author "Julius Smith";
// License: STK-4.3

bypass1(bpc,e) = _ <: select2(bpc,(inswitch:e),_)
with {
    inswitch = select2(bpc,_,0);
};


//-------------------`(ba.)bypass2`-------------------------
// Takes a stereo input signal, route it to `e` and bypass it if `bpc = 1`.
// When bypassed, `e` is feed with zeros so that its state is cleanup up.
// `bypass2` is a standard Faust function.
//
// #### Usage
//
// ```
// _,_ : bypass2(bpc,e) : _,_
// ```
//
// Where:
//
// * `bpc`: bypass switch (0/1)
// * `e`: a stereo effect
//------------------------------------------------------------
declare bypass2 author "Julius Smith";
// License: STK-4.3

bypass2(bpc,e) = _,_ <: ((inswitch:e),_,_) : select2stereo(bpc)
with {
    inswitch = _,_ : (select2(bpc,_,0), select2(bpc,_,0)) : _,_;
};


//-------------------`(ba.)bypass1to2`-------------------------
// Bypass switch for effect `e` having mono input signal and stereo output.
// Effect `e` is bypassed if `bpc = 1`.When bypassed, `e` is feed with zeros
// so that its state is cleanup up.
// `bypass1to2` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : bypass1to2(bpc,e) : _,_
// ```
//
// Where:
//
// * `bpc`: bypass switch (0/1)
// * `e`: a mono-to-stereo effect
//------------------------------------------------------------
declare bypass1to2 author "Julius Smith";
// License: STK-4.3

bypass1to2(bpc,e) = _ <: ((inswitch:e),_,_) : select2stereo(bpc)
with {
    inswitch = select2(bpc,_,0);
};


//-------------------`(ba.)bypass_fade`-------------------------
// Bypass an arbitrary (N x N) circuit with 'n' samples crossfade.
// Inputs and outputs signals are faded out when 'e' is bypassed,
// so that 'e' state is cleanup up.
// Once bypassed the effect is replaced by `par(i,N,_)`.
// Bypassed circuits can be chained.
//
// #### Usage
//
// ```
// _ : bypass_fade(n,b,e) : _
// or
// _,_ : bypass_fade(n,b,e) : _,_
// ```
// * `n`: number of samples for the crossfade
// * `b`: bypass switch (0/1)
// * `e`: N x N circuit
//
// #### Example test program
//
// ```
// process = bypass_fade(ma.SR/10, checkbox("bypass echo"), echo);
// process = bypass_fade(ma.SR/10, checkbox("bypass reverb"), freeverb);
// ```
//---------------------------------------------------------------
bypass_fade(n, b, e) = par(i, ins, _)
            <: (par(i, ins, *(1-xb)) : e : par(i, outs, *(1-xb))), par(i, ins, *(xb))
            :> par(i, outs, _)
with {
    ins = inputs(e);
    outs = outputs(e);
    xb = ramp(n, b);
};


//----------------------------`(ba.)toggle`------------------------------------------
// Triggered by the change of 0 to 1, it toggles the output value
// between 0 and 1.
//
// #### Usage
//
// ```
// _ : toggle : _
// ```
// #### Example test program
//
// ```
// button("toggle") : toggle : vbargraph("output", 0, 1)
// (an.amp_follower(0.1) > 0.01) : toggle : vbargraph("output", 0, 1) // takes audio input
// ```
//
//------------------------------------------------------------------------------
declare toggle author "Vince";

toggle = trig : loop
with {
    trig(x) = (x-x') == 1;
    loop = != ~ _;
};


//----------------------------`(ba.)on_and_off`------------------------------------------
// The first channel set the output to 1, the second channel to 0.
//
// #### Usage
//
// ```
// _,_ : on_and_off : _
// ```
//
// #### Example test program
//
// ```
// button("on"), button("off") : on_and_off : vbargraph("output", 0, 1)
// ```
//
//------------------------------------------------------------------------------
declare on_and_off author "Vince";

on_and_off(a, b) = (a : trig) : loop(b)
with {
    trig(x) = (x-x') == 1;
    loop(b) = + ~ (_ >= 1) * ((b : trig) == 0);
};


//----------------------------`(ba.)bitcrusher`------------------------------------------
// Produce distortion by reduction of the signal resolution.
//
// #### Usage
//
// ```
// _ : bitcrusher(nbits) : _
// ```
//
// Where:
//
// * `nbits`: the number of bits of the wanted resolution
//
//------------------------------------------------------------------------------
declare bitcrusher author "Julius O. Smith III, revised by Stephane Letz";

bitcrusher(nbits,x) = round(x * scaler) / scaler
with {
    scaler = float(2^nbits - 1);
};

//----------------------------`(ba.)mulaw_bitcrusher`------------------------------------------
// Produce distortion by reducing the signal resolution using μ-law compression.
//
// #### Usage
//
// ```
// _ : mulaw_bitcrusher(mu,nbits) : _
// ```
//
// Where:
//
// * `mu`: controls the degree of μ-law compression, larger values result in stronger compression
// * `nbits`: the number of bits of the wanted resolution
//
// #### Description
//
// The `mulaw_bitcrusher` applies a combination of μ-law compression, quantization, and expansion 
// to create a non-linear bitcrushed effect. This method retains finer detail in lower-amplitude signals 
// compared to linear bitcrushing, making it suitable for creative sound design.
//
// #### Theory
//
// 1. **μ-law Compression**:
//    emphasizes lower-amplitude signals by applying a logarithmic curve to the signal.
//    The formula used is:
//    ```
//    F(x) = ma.signum(x) * log(1 + mu * abs(x)) / log(1 + mu);
//    ```
// 2. **Quantization**:
//    reduces the signal resolution to `nbits` by rounding values to the nearest step within the specified bit depth.
//
// 3. **μ-law Expansion**:
//    reverses the compression applied earlier to restore the signal to its original dynamic range:
//    ```
//    F⁻¹(y) = ma.signum(y) * (pow(1 + mu, abs(y)) - 1) / mu;
//    ```
//
// #### Example test program
//
// ```
// process = os.osc(440) : mulaw_bitcrusher(255, 8);
// ```
//
// In this example, a sine wave at 440 Hz is passed through the μ-law bitcrusher, with a compression
// parameter `mu` of 255 and 8-bit quantization. This creates a distorted, "lo-fi" effect.
//
// #### References
//
// * <https://en.wikipedia.org/wiki/Μ-law_algorithm>

//------------------------------------------------------------------------------
declare mulaw_bitcrusher author "Stephane Letz";

mulaw_bitcrusher(mu,nbits,x) =  x : muLawCompress(mu) : bitcrusher(nbits) : muLawExpand(mu)
with {
    // μ-law compression
    muLawCompress(mu, x) = ma.signum(x) * log(1 + mu * abs(x)) / log(1 + mu);

    // μ-law expansion (decompression)
    muLawExpand(mu, x) = ma.signum(x) * (pow(1 + mu, abs(x)) - 1) / mu;
};
   

//=================================Sliding Reduce=========================================
// Provides various operations on the last n samples using a high order
// `slidingReduce(op,n,maxN,disabledVal,x)` fold-like function:
//
// * `slidingSum(n)`: the sliding sum of the last n input samples, CPU-light
// * `slidingSump(n,maxN)`: the sliding sum of the last n input samples, numerically stable "forever"
// * `slidingMax(n,maxN)`: the sliding max of the last n input samples
// * `slidingMin(n,maxN)`: the sliding min of the last n input samples
// * `slidingMean(n)`: the sliding mean of the last n input samples, CPU-light
// * `slidingMeanp(n,maxN)`: the sliding mean of the last n input samples, numerically stable "forever"
// * `slidingRMS(n)`: the sliding RMS of the last n input samples, CPU-light
// * `slidingRMSp(n,maxN)`: the sliding RMS of the last n input samples, numerically stable "forever"
//
// #### Working Principle
//
// If we want the maximum of the last 8 values, we can do that as:
//
// ```
// simpleMax(x) =
//  (
//    (
//      max(x@0,x@1),
//      max(x@2,x@3)
//    ) :max
//  ),
//  (
//    (
//      max(x@4,x@5),
//      max(x@6,x@7)
//    ) :max
//  )
//  :max;
// ```
//
// `max(x@2,x@3)` is the same as `max(x@0,x@1)@2` but the latter re-uses a
// value we already computed,so is more efficient. Using the same trick for
// values 4 trough 7, we can write:
//
// ```
// efficientMax(x)=
//  (
//    (
//      max(x@0,x@1),
//      max(x@0,x@1)@2
//    ) :max
//  ),
//  (
//    (
//      max(x@0,x@1),
//      max(x@0,x@1)@2
//    ) :max@4
//  )
//  :max;
// ```
//
// We can rewrite it recursively, so it becomes possible to get the maximum at
// have any number of values, as long as it's a power of 2.
//
// ```
// recursiveMax =
//  case {
//    (1,x) => x;
//    (N,x) => max(recursiveMax(N/2,x), recursiveMax(N/2,x)@(N/2));
//  };
// ```
//
// What if we want to look at a number of values that's not a power of 2?
// For each value, we will have to decide whether to use it or not.
// If n is bigger than the index of the value, we use it, otherwise we replace
// it with (`0-(ma.MAX)`):
//
// ```
// variableMax(n,x) =
//  max(
//    max(
//      (
//        (x@0 : useVal(0)),
//        (x@1 : useVal(1))
//      ):max,
//      (
//        (x@2 : useVal(2)),
//        (x@3 : useVal(3))
//      ):max
//    ),
//    max(
//      (
//        (x@4 : useVal(4)),
//        (x@5 : useVal(5))
//      ):max,
//      (
//        (x@6 : useVal(6)),
//        (x@7 : useVal(7))
//      ):max
//    )
//  )
// with {
//  useVal(i) = select2((n>=i) , (0-(ma.MAX)),_);
// };
// ```
//
// Now it becomes impossible to re-use any values. To fix that let's first look
// at how we'd implement it using recursiveMax, but with a fixed n that is not
// a power of 2. For example, this is how you'd do it with `n=3`:
//
// ```
// binaryMaxThree(x) =
//  (
//    recursiveMax(1,x)@0, // the first x
//    recursiveMax(2,x)@1  // the second and third x
//  ):max;
// ```
//
// `n=6`
//
// ```
// binaryMaxSix(x) =
//  (
//    recursiveMax(2,x)@0, // first two
//    recursiveMax(4,x)@2  // third trough sixth
//  ):max;
// ```
//
// Note that `recursiveMax(2,x)` is used at a different delay then in
// `binaryMaxThree`, since it represents 1 and 2, not 2 and 3. Each block is
// delayed the combined size of the previous blocks.
//
// `n=7`
//
// ```
// binaryMaxSeven(x) =
//  (
//    (
//      recursiveMax(1,x)@0, // first x
//      recursiveMax(2,x)@1  // second and third
//    ):max,
//    (
//      recursiveMax(4,x)@3  // fourth trough seventh
//    )
//  ):max;
// ```
//
// To make a variable version, we need to know which powers of two are used,
// and at which delay time.
//
// Then it becomes a matter of:
//
// * lining up all the different block sizes in parallel: `sequentialOperatorParOut()`
// * delaying each the appropriate amount: `sumOfPrevBlockSizes()`
// * turning it on or off: `useVal()`
// * getting the maximum of all of them: `parallelOp()`
//
// In Faust, we can only do that for a fixed maximum number of values: `maxN`, known at compile time.

//========================================================================================
// Section contributed by Bart Brouns (bart@magnetophon.nl).
// SPDX-License-Identifier: GPL-3.0
// Copyright (C) 2018 Bart Brouns


//-----------------------------`(ba.)slidingReduce`-----------------------------
// Fold-like high order function. Apply a commutative binary operation `op` to
// the last `n` consecutive samples of a signal `x`. For example :
// `slidingReduce(max,128,128,0-(ma.MAX))` will compute the maximum of the last
// 128 samples. The output is updated each sample, unlike reduce, where the
// output is constant for the duration of a block.
//
// #### Usage
//
// ```
// _ : slidingReduce(op,n,maxN,disabledVal) : _
// ```
//
// Where:
//
// * `n`: the number of values to process
// * `maxN`: the maximum number of values to process (int, known at compile time, maxN > 0)
// * `op`: the operator. Needs to be a commutative one.
// * `disabledVal`: the value to use when we want to ignore a value.
//
// In other words, `op(x,disabledVal)` should equal to `x`. For example,
// `+(x,0)` equals `x` and `min(x,ma.MAX)` equals `x`. So if we want to
// calculate the sum, we need to give 0 as `disabledVal`, and if we want the
// minimum, we need to give `ma.MAX` as `disabledVal`.
//------------------------------------------------------------------------------
slidingReduce(op,n,0,disabledVal) = 0:!;
slidingReduce(op,n,1,disabledVal) = _;
slidingReduce(op,n,maxN,disabledVal) =
    sequentialOperatorParOut(maxNrBits(maxN)-1,op) : par(i, maxNrBits(maxN), _@sumOfPrevBlockSizes(i) : useVal(i)) : parallelOp(op, maxNrBits(maxN))
    with {
       sequentialOperatorParOut(N,op) = seq(i, N, operator(i));
        operator(i) = si.bus(i), (_<: _ , op(_,_@(pow2(i))));

        // The sum of all the sizes of the previous blocks
        sumOfPrevBlockSizes(0) = 0;
        sumOfPrevBlockSizes(i) = (ba.subseq((allBlockSizes),0,i):>_);

        allBlockSizes = par(i, maxNrBits(maxN-1), (pow2(i)) * isUsed(i));
        maxNrBits(n) = int2nrOfBits(n);

        // Decide wether or not to use a certain value, based on n
        useVal(i) = select2(isUsed(i), disabledVal, _);

        isUsed(i) = ba.take(i+1, (int2bin(n,(maxN-1)*2+1)));
        pow2(i) = 1<<i;
        // same as:
        // pow2(i) = int(pow(2,i));
        // but in the block diagram, it will be displayed as a number, instead of a formula

        // convert n into a list of ones and zeros
        int2bin(n,maxN) = par(j, maxNrBits(maxN-1), int(floor((n)/(pow2(j))))%2);
        // calculate how many ones and zeros are needed to represent maxN
        int2nrOfBits(n) = int(floor(log(n)/log(2))+1);
    };


//------------------------------`(ba.)slidingSum`------------------------------
// The sliding sum of the last n input samples.
//
// It will eventually run into numerical trouble when there is a persistent dc component.
// If that matters in your application, use the more CPU-intensive `ba.slidingSump`.
//
// #### Usage
//
// ```
// _ : slidingSum(n) : _
// ```
//
// Where:
//
// * `n`: the number of values to process
//------------------------------------------------------------------------------
slidingSum(n) = fi.integrator <: _, _@int(max(0,n)) :> -;


//------------------------------`(ba.)slidingSump`------------------------------
// The sliding sum of the last n input samples.
//
// It uses a lot more CPU than `ba.slidingSum`, but is numerically stable "forever" in return.
//
// #### Usage
//
// ```
// _ : slidingSump(n,maxN) : _
// ```
//
// Where:
//
// * `n`: the number of values to process
// * `maxN`: the maximum number of values to process (int, known at compile time, maxN > 0)
//------------------------------------------------------------------------------
slidingSump(n,maxN) = slidingReduce(+,n,maxN,0);


//----------------------------`(ba.)slidingMax`--------------------------------
// The sliding maximum of the last n input samples.
//
// #### Usage
//
// ```
// _ : slidingMax(n,maxN) : _
// ```
//
// Where:
//
// * `n`: the number of values to process
// * `maxN`: the maximum number of values to process (int, known at compile time, maxN > 0)
//------------------------------------------------------------------------------
slidingMax(n,maxN) = slidingReduce(max,n,maxN,0-(ma.MAX));

//----------------------------`(ba.)slidingMin`--------------------------------
// The sliding minimum of the last n input samples.
//
// #### Usage
//
// ```
// _ : slidingMin(n,maxN) : _
// ```
//
// Where:
//
// * `n`: the number of values to process
// * `maxN`: the maximum number of values to process (int, known at compile time, maxN > 0)
//------------------------------------------------------------------------------
slidingMin(n,maxN) = slidingReduce(min,n,maxN,ma.MAX);


//----------------------------`(ba.)slidingMean`-------------------------------
// The sliding mean of the last n input samples.
//
// It will eventually run into numerical trouble when there is a persistent dc component.
// If that matters in your application, use the more CPU-intensive `ba.slidingMeanp`.
//
// #### Usage
//
// ```
// _ : slidingMean(n) : _
// ```
//
// Where:
//
// * `n`: the number of values to process
//------------------------------------------------------------------------------
slidingMean(n) = slidingSum(n)/n;


//----------------------------`(ba.)slidingMeanp`-------------------------------
// The sliding mean of the last n input samples.
//
// It uses a lot more CPU than `ba.slidingMean`, but is numerically stable "forever" in return.
//
// #### Usage
//
// ```
// _ : slidingMeanp(n,maxN) : _
// ```
//
// Where:
//
// * `n`: the number of values to process
// * `maxN`: the maximum number of values to process (int, known at compile time, maxN > 0)
//------------------------------------------------------------------------------
slidingMeanp(n,maxN) = slidingSump(n,maxN)/n;


//---------------------------`(ba.)slidingRMS`---------------------------------
// The root mean square of the last n input samples.
//
// It will eventually run into numerical trouble when there is a persistent dc component.
// If that matters in your application, use the more CPU-intensive `ba.slidingRMSp`.

//
// #### Usage
//
// ```
// _ : slidingRMS(n) : _
// ```
//
// Where:
//
// * `n`: the number of values to process
//------------------------------------------------------------------------------
slidingRMS(n) = pow(2) : slidingMean(n) : sqrt;


//---------------------------`(ba.)slidingRMSp`---------------------------------
// The root mean square of the last n input samples.
//
// It uses a lot more CPU than `ba.slidingRMS`, but is numerically stable "forever" in return.
//
// #### Usage
//
// ```
// _ : slidingRMSp(n,maxN) : _
// ```
//
// Where:
//
// * `n`: the number of values to process
// * `maxN`: the maximum number of values to process (int, known at compile time, maxN > 0)
//------------------------------------------------------------------------------
slidingRMSp(n,maxn) = pow(2) : slidingMeanp(n,maxn) : sqrt;


//========================================================================================
// section contributed by Bart Brouns (bart@magnetophon.nl).
// spdx-license-identifier: gpl-3.0
// copyright (c) 2020 Bart Brouns

//=================================Parallel Operators=========================================
// Provides various operations on N parallel inputs using a high order
// `parallelOp(op,N,x)` function:
//
// * `parallelMax(N)`: the max of n parallel inputs
// * `parallelMin(N)`: the min of n parallel inputs
// * `parallelMean(N)`: the mean of n parallel inputs
// * `parallelRMS(N)`: the RMS of n parallel inputs

//-----------------------------`(ba.)parallelOp`-----------------------------
// Apply a commutative binary operation `op` to N parallel inputs.
//
// #### usage
//
// ```
// si.bus(N) : parallelOp(op,N) : _
// ```
//
// where:
//
// * `N`: the number of parallel inputs known at compile time
// * `op`: the operator which needs to be commutative
//
//------------------------------------------------------------------------------

parallelOp(op,1) = _;
parallelOp(op,2) = op;
parallelOp(op,n) = op(parallelOp(op,n-1));

declare parallelOp author "Bart Brouns";
declare parallelOp licence "GPL-3.0";
declare parallelOp copyright "Copyright (c) 2020 Bart Brouns <bart@magnetophon.nl>";


//---------------------------`(ba.)parallelMax`---------------------------------
// The maximum of N parallel inputs.
//
// #### Usage
//
// ```
// si.bus(N) : parallelMax(N) : _
// ```
//
// Where:
//
// * `N`: the number of parallel inputs known at compile time
//------------------------------------------------------------------------------
parallelMax(n) = parallelOp(max,n);

declare parallelMax author "Bart Brouns";
declare parallelMax licence "GPL-3.0";
declare parallelMax copyright "Copyright (c) 2020 Bart Brouns <bart@magnetophon.nl>";


//---------------------------`(ba.)parallelMin`---------------------------------
// The minimum of N parallel inputs.
//
// #### Usage
//
// ```
// si.bus(N) : parallelMin(N) : _
// ```
//
// Where:
//
// * `N`: the number of parallel inputs known at compile time
//------------------------------------------------------------------------------
parallelMin(n) = parallelOp(min,n);

declare parallelMin author "Bart Brouns";
declare parallelMin licence "GPL-3.0";
declare parallelMin copyright "Copyright (c) 2020 Bart Brouns <bart@magnetophon.nl>";


//---------------------------`(ba.)parallelMean`---------------------------------
// The mean of N parallel inputs.
//
// #### Usage
//
// ```
// si.bus(N) : parallelMean(N) : _
// ```
//
// Where:
//
// * `N`: the number of parallel inputs known at compile time
//------------------------------------------------------------------------------
parallelMean(n) = si.bus(n):>_/n;

declare parallelMean author "Bart Brouns";
declare parallelMean licence "GPL-3.0";
declare parallelMean copyright "Copyright (c) 2020 Bart Brouns <bart@magnetophon.nl>";

//---------------------------`(ba.)parallelRMS`---------------------------------
// The RMS of N parallel inputs.
//
// #### Usage
//
// ```
// si.bus(N) : parallelRMS(N) : _
// ```
//
// Where:
//
// * `N`: the number of parallel inputs known at compile time
//------------------------------------------------------------------------------
parallelRMS(n) = par(i, n, pow(2)) : parallelMean(n) : sqrt;

declare parallelRMS author "Bart Brouns";
declare parallelRMS licence "GPL-3.0";
declare parallelRMS copyright "Copyright (c) 2020 Bart Brouns <bart@magnetophon.nl>";

//////////////////////////////////Deprecated Functions////////////////////////////////////
// This section implements functions that used to be in music.lib but that are now
// considered as "deprecated".
//////////////////////////////////////////////////////////////////////////////////////////

millisec = ma.SR/1000.0;

time1s = hslider("time", 0, 0,  1000, 0.1)*millisec;
time2s = hslider("time", 0, 0,  2000, 0.1)*millisec;
time5s = hslider("time", 0, 0,  5000, 0.1)*millisec;
time10s = hslider("time", 0, 0, 10000, 0.1)*millisec;
time21s = hslider("time", 0, 0, 21000, 0.1)*millisec;
time43s = hslider("time", 0, 0, 43000, 0.1)*millisec;