//#################################### 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 "0.9";

//=============================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`----------
// Converts a loudness in dB to a linear gain (0-1).
// `db2linear` is a standard Faust function.
//
// #### Usage
//
// ```
// db2linear(l) : _
// ```
//
// Where:
//
// * `l`: loudness in dB
//-----------------------------
db2linear(l) = pow(10.0, l/20.0);


//-------`(ba.)linear2db`----------
// Converts a linear gain (0-1) to a loudness 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)~_) - 1;


//-------`(ba.)ramp`----------
// A linear ramp with a slope of '(+/-)1/n' samples to reach the next 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 linear ramp to reach a next value in 'n' samples.
// Note that the interpolation process is restarted every time
// the desired output value changes, the interpolation time is sampled only then.
//
// #### Usage
//
// ```
// _ : line(n) : _
// ```
// Where:
//
// * `n`: number of samples to reach the next value
//------------------------
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.)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`.
// `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
//
// * `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(ba.time)*(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), 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";


//============================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.)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.)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.)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.)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.)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
//
// 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}`
//
// `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)));
};


//-------------------`(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);
    round(x) = floor(x + 0.5);
};


//=================================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;