//#################################### signals.lib ########################################
// A library of basic elements to handle signals in Faust. Its official prefix is `si`.
//
// #### References
// * <https://github.com/grame-cncm/faustlibraries/blob/master/signals.lib>
//########################################################################################

/************************************************************************
************************************************************************
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.
************************************************************************
************************************************************************/

ba = library("basics.lib");
ro = library("routes.lib");
ma = library("maths.lib");
si = library("signals.lib");

declare name "Faust Signal Routing Library";
declare version "1.6.0";

//=============================Functions Reference========================================
//========================================================================================

//--------------------------------`(si.)bus`-------------------------------------
// Put N cables in parallel.
// `bus` is a standard Faust function.
//
// #### Usage
//
// ```
// bus(N)
// bus(4) : _,_,_,_
// ```
//
// Where:
//
// * `N`: is an integer known at compile time that indicates the number of parallel cables
//-----------------------------------------------------------------------------
bus(0) = 0:!;
bus(1) = _;    // avoids a lot of "bus(1)" labels in block diagrams
bus(2) = _,_;  // avoids a lot of "bus(2)" labels in block diagrams
bus(N) = par(i, N, _);

//--------------`(si.)block`--------------
// Block - terminate N signals.
// `block` is a standard Faust function.
//
// #### Usage
//
// ```
// si.bus(N) : block(N)
// ```
//
// Where:
//
// * `N`: the number of signals to be blocked known at compile time 
//--------------------------------------
block(N) = par(i, N, !);

//-----------------------------`(si.)interpolate`-------------------------------
// Linear interpolation between two signals.
//
// #### Usage
//
// ```
// _,_ : interpolate(i) : _
// ```
//
// Where:
//
// * `i`: interpolation control between 0 and 1 (0: first input; 1: second input)
//-----------------------------------------------------------------------------
interpolate(i,x,y) = x + i*(y-x);  

//---------------`(si.)repeat`----------------------------------------------
// Repeat an effect N time(s) and take the parallel sum of all
// intermediate buses.
//
// #### References
// * <https://github.com/orlarey/presentation-compilateur-faust/blob/master/slides.pdf>
//
// #### Usage
//
// ```
// si.bus(inputs(FX)) : repeat(N, FX) : si.bus(outputs(FX))
// ```
//
// Where:
//
// * `N`: Number of repetitions, minimum of 1, a constant numerical expression 
// * `FX`: an arbitrary effect (N inputs and N outputs) that will be repeated
//
// Example 1:
// ```
// process = repeat(2, dm.zita_light) : _*.5,_*.5;
// ```
//
// Example 2:
// ```
// N = 4;
// C = 2;
// fx(i) = i+1, par(j, C, @(i*5000));
// process = 0, si.bus(C) : repeat(N, fx) : !, par(i, C, _*.2/N);
// ```
//---------------------------------------------------------------------------
declare repeat author "Yann Orlarey and Stéphane Letz, revised by David Braun";

repeat(1, FX) = FX;
repeat(n, FX) = FX <: si.bus(N), repeat(n-1, FX) :> si.bus(N)
with {
    N = outputs(FX);
};

//------------------------`(si.)smoo`---------------------------------------
// Smoothing function based on `smooth` ideal to smooth UI signals
// (sliders, etc.) down. Approximately, this is a 7 Hz one-pole
// low-pass considering the coefficient calculation:
//    exp(-2pi*CF/SR).
//
// `smoo` is a standard Faust function.
//
// #### Usage
//
// ```
// hslider(...) : smoo;
// ```
//---------------------------------------------------------------------
smoo = si.smooth(1 - 44.1/ma.SR);

//-----------------------`(si.)polySmooth`--------------------------------
// A smoothing function based on `smooth` that doesn't smooth when a
// trigger signal is given. This is very useful when making
// polyphonic synthesizer to make sure that the value of the parameter
// is the right one when the note is started.
//
// #### Usage
//
// ```
// hslider(...) : polySmooth(g,s,d) : _
// ```
//
// Where:
//
// * `g`: the gate/trigger signal used when making polyphonic synths
// * `s`: the smoothness (see `smooth`)
// * `d`: the number of samples to wait before the signal start being
//     smoothed after `g` switched to 1
//-------------------------------------------------------------------
polySmooth(g,s,d) = si.smooth(s*((g==(g@d)) | (g == 0)));

//-----------------------`(si.)smoothAndH`--------------------------------
// A smoothing function based on `smooth` that holds its output
// signal when a trigger is sent to it. This feature is convenient
// when implementing polyphonic instruments to prevent some
// smoothed parameter to change when a note-off event is sent.
//
// #### Usage
//
// ```
// hslider(...) : smoothAndH(g,s) : _
// ```
//
// Where:
//
// * `g`: the hold signal (0 for hold, 1 for bypass)
// * `s`: the smoothness (see `smooth`)
//-------------------------------------------------------------------
smoothAndH(t,s) = si.smooth(s*t) : ba.sAndH(t);

//-----------------------------`(si.)bsmooth`------------------------------
// Block smooth linear interpolation during a block of samples (given by the `ma.BS` value).
//
// #### Usage
//
// ```
// hslider(...) : bsmooth : _
// ```
//-----------------------------------------------------------------------
bsmooth(c) = +(i) ~ _
with {
    i = (c-c@n)/n;
    n = min(4096, max(1, ma.BS));
};

//-------------------------------`(si.)dot`--------------------------------------
// Dot product for two vectors of size N.
//
// #### Usage
//
// ```
// si.bus(N), si.bus(N) : dot(N) : _
// ```
//
// Where:
//
// * `N`: size of the vectors (int, must be known at compile time)
//-----------------------------------------------------------------------------
dot(N) = ro.interleave(N,2) : par(i,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.

All 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!)

************************************************************************/
//-------------------`(si.)smooth`-----------------------------------
// Exponential smoothing by a unity-dc-gain one-pole lowpass.
// `smooth` is a standard Faust function.
//
// #### Usage:
//
// ```
// _ : si.smooth(ba.tau2pole(tau)) : _
// ```
//
// Where:
//
// * `tau`: desired smoothing time constant in seconds, or
//
// ```
// hslider(...) : smooth(s) : _
// ```
//
// Where:
//
// * `s`: smoothness between 0 and 1. s=0 for no smoothing, s=0.999 is "very smooth",
// s>1 is unstable, and s=1 yields the zero signal for all inputs.
// The exponential time-constant is approximately 1/(1-s) samples, when s is close to
// (but less than) 1.
//
// #### References:
//
// * <https://ccrma.stanford.edu/~jos/mdft/Convolution_Example_2_ADSR.html>
// * <https://ccrma.stanford.edu/~jos/aspf/Appendix_B_Inspecting_Assembly.html>
//-------------------------------------------------------------

// See [grame-cncm/faustlibraries]: Minor improvement to si.smoo. (Discussion #106)

smooth_imp = case { 
  // y[n] = (1 - s) * x[n] + s * y[n - 1]
  (0,s) => \(x).(fb ~ _ with { fb(y) = (1.0 - s) * x + s * y; });
  
  // y[n] = s * (y[n - 1] - x[n]) + x[n]
  (1,s) => \(x).(fb ~ _ with { fb(y) = s * (y - x) + x; });
  
  // y[n] = y[n - 1] + (1 - s) * (x[n] - y[n - 1])
  (2,s) => \(x).(fb ~ _ with { fb(y) = y + (1.0 - s) * (x - y); });
};

// The best compromise on modern CPUs where two independent multiplications can be done in parallel.
// Other versions could possibly be faster on embedded devices. 
smooth = smooth_imp(0);


//--------------------------------`(si.)smoothq`-------------------------------------
// Smoothing with continuously variable curves from Exponential to Linear, with a constant time.
//
// #### Usage
//
// ```
// _ : smoothq(time, q) : _;
// ```
//
// Where:
//
// * `time`: seconds to reach target
// * `q`: curve shape (between 0..1, 0 is Exponential, 1 is Linear)
//-----------------------------------------------------------------------------
declare smoothq author "Andrew John March";
declare smoothq licence "STK-4.3";

smoothq(time, q, tar) = envelope 
    with {
        ratio = pow(q, 5) * 32 : max(0.001);
        coef = (ratio / (1 + ratio))^(1 / (ma.SR * time));
        fb(cBase, cTar, cSrc, cDistance, cDir, cY, cOut) = base, tar, src, distance, dir, y, out
            with {
                trig = cTar != tar;
                src = select2(trig, cSrc, cOut);
                dir = select2(trig, cDir, tar > src);
                base = select2(trig, cBase, select2(dir,
                    (0.0 - ratio) * (1.0 - coef),
                    (1.0 + ratio) * (1.0 - coef)
                ));
                distance = select2(trig, cDistance, abs(tar - src));
                y = select2(dir,
                    select2(trig, 
                        max( 0, base + cY * coef),
                        1
                    ),
                    select2(trig,
                        min( 1, base + cY * coef),
                        0
                    )
                );
                out = select2(dir,
                    tar + y * distance,
                    src + y * distance
                );
            };
        envelope = fb ~ (_,_,_,_,_,_,_) : !,!,!,!,!,!,_;
    };


//--------------------------------`(si.)cbus`-------------------------------------
// N parallel cables for complex signals.
// `cbus` is a standard Faust function.
//
// #### Usage
//
// ```
// cbus(N)
// cbus(4) : (r0,i0), (r1,i1), (r2,i2), (r3,i3)
// ```
//
// Where:
//
// * `N`: is an integer known at compile time that indicates the number of parallel cables.
// * each complex number is represented by two real signals as (real,imag)
//-----------------------------------------------------------------------------
cbus(1) = (_,_);
cbus(N) = par(i, N, (_,_));

//--------------------------------`(si.)cmul`-------------------------------------
// Multiply two complex signals pointwise.
// `cmul` is a standard Faust function.
//
// #### Usage
//
// ```
// (r1,i1) : cmul(r2,i2) : (_,_)
// ```
//
// Where:
//
// * Each complex number is represented by two real signals as (real,imag), so
// - `(r1,i1)` = real and imaginary parts of signal 1
// - `(r2,i2)` = real and imaginary parts of signal 2
//-----------------------------------------------------------------------------
cmul(r1,i1,r2,i2) = (r1*r2 - i1*i2), (r1*i2 + r2*i1);

//--------------------------------`(si.)cconj`-------------------------------------
// Complex conjugation of a (complex) signal.
// `cconj` is a standard Faust function.
//
// #### Usage
//
// ```
// (r1,i1) : cconj : (_,_)
// ```
//
// Where:
//
// * Each complex number is represented by two real signals as (real,imag), so
// - `(r1,i1)` = real and imaginary parts of the input signal
// - `(r1,-i1)` = real and imaginary parts of the output signal
//-----------------------------------------------------------------------------
cconj = _, *(-1);

// end jos section
/************************************************************************
FAUST library file, further contributions section
All contributions below should indicate both the contributor and terms
of license.  If no such indication is found, "git blame" will say who
last edited each line, and that person can be emailed to inquire about
license disposition, if their license choice is not already indicated
elsewhere among the libraries.  It is expected that all software will be
released under LGPL, STK-4.3, MIT, BSD, or a similar FOSS license.
************************************************************************/

//-------------`(si.)onePoleSwitching`---------------
// One pole filter with independent attack and release times.
//
// #### Usage
//
// ```
// _ : onePoleSwitching(att,rel) : _
// ```
// 
// Where:
//
// * `att`: the attack tau time constant in second
// * `rel`: the release tau time constant in second
//
//----------------------------------------------------
declare onePoleSwitching author "Jonatan Liljedahl, revised by Dario Sanfilippo";
declare onePoleSwitching licence "STK-4.3";

onePoleSwitching(att, rel, x) = loop ~ _
    with {
        loop(yState) = (1.0 - coeff) * x + coeff * yState
            with {
                coeff = ba.if(x > yState, ba.tau2pole(att), ba.tau2pole(rel));
            };
    };
    
// Kept for backward compatibility
lag_ud = onePoleSwitching;


// end further further contributions section

//-------------`(si.)rev`---------------
// Reverse the input signal by blocks of n>0 samples. `rev(1)` is the indentity
// function. `rev(n)` has a latency of `n-1` samples.
//
// #### Usage
//
// ```
// _ : rev(n) : _
// ```
//
// Where:
//
// * `n`: the block size in samples
//----------------------------------------------------
declare rev author "Yann Orlarey";

rev(n) = @(phase(n)*2) 
with {
    phase(n) = 1 : (+ : %(n)) ~ _ : max(0) : min(n-1);
};

//--------------------`(si.)vecOp`----------------------------------------------
//
// This function is a generalisation of Faust's iterators such as `prod` and
// `sum`, and it allows to perform operations on an arbitrary number of
// vectors, provided that they all have the same length. Unlike Faust's
// iterators `prod` and `sum` where the vector size is equal to one and the 
// vector space dimension must be specified by the user, this function will 
// infer the vector space dimension and vector size based on the vectors list 
// that we provide.
//
// The outputs of the function are equal to the vector size, whereas the
// number of inputs is dependent on whether the elements of the vectors
// provided expect an incoming signal themselves or not. We will see a
// clarifying example later; in general, the number of total inputs will
// be the sum of the inputs in each input vector.
//
// Note that we must provide a list of at least two vectors, each with a size 
// that is greater or equal to one.
//
// #### Usage
//
// ```
//      si.bus(inputs(vectorsList)) : vecOp((vectorsList), op) : si.bus(outputs(ba.take(1, vectorsList)));
// ```
//
// #### Where
//
// * `vectorsList`: is a list of vectors
// * `op`: is a two-input, one-output operator
//
// For example, consider the following vectors lists:
//
//      v0 = (0 , 1 , 2 , 3);
//      v1 = (4 , 5 , 6 , 7);
//      v2 = (8 , 9 , 10 , 11);
//      v3 = (12 , 13 , 14 , 15);
//      v4 = (+(16) , _ , 18 , *(19));
//      vv = (v0 , v1 , v2 , v3);
//
// Although Faust has limitations for list processing, these vectors can be
// combined or processed individually.
//
// If we do:
// 
//      process = vecOp(v0, +);
//
// the function will deduce a vector space of dimension equal to four and 
// a vector length equal to one. Note that this is equivalent to writing:
//
//      process = v0 : sum(i, 4, _);
//
// Similarly, we can write:
//
//      process = vecOp((v0 , v1), *) :> _;
//
// and we have a dimension-two space and length-four vectors. This is the dot 
// product between vectors v0 and v1, which is equivalent to writing:
//
//      process = v0 , v1 : dot(4);
//
// The examples above have no inputs, as none of the elements of the vectors
// expect inputs. On the other hand, we can write:
//
//      process = vecOp((v4 , v4), +);
//
// and the function will have six inputs and four outputs, as each vector
// has three of the four elements expecting an input, times two, as the two
// input vectors are identical.
//
// Finally, we can write:
//
//      process = vecOp(vv, &);
//
// to perform the bitwise AND on all the elements at the same position in 
// each vector, having dimension equal to the vector length equal to four.
// 
// Or even:
//
//      process = vecOp((vv , vv), &);
//
// which gives us a dimension equal to two, and a vector size equal to sixteen.
//
// For a more practical use-case, this is how we can implement a time-invariant
// feedback delay network with Hadamard matrix:
//
//      N = 4;
//      normalisation = 1.0 / sqrt(N);
//      coeffVec = par(i, N, .99 * normalisation);
//      delVec = par(i, N, (i + 1) * 3);
//      process = vecOp((si.bus(N) , si.bus(N)), +) ~ 
//          vecOp((vecOp((ro.hadamard(N) , coeffVec), *) , delVec), @);
//
//------------------------------------------------------------------------------
declare vecOp author "Dario Sanfilippo";
declare vecOp copyright "Copyright (C) 2022 Dario Sanfilippo
      <sanfilippo.dario@gmail.com>";
declare vecOp license "MIT License";
vecOp(vectorsList, op) =
    vectorsList : seq(i, vecDim - 1, vecOp2D , vecBus(vecDim - 2 - i))
    with {
        vecBus(0) = par(i, vecLen, 0 : !);
        vecBus(dim) = par(i, dim, si.bus(vecLen));
        vecOp2D = ro.interleave(vecLen, 2) : par(i, vecLen, op);
        vecDim = outputs(vectorsList) / vecLen;
        vecLen = outputs(ba.take(1, vectorsList));
    };

/*
Balanced par, sum, and prod, faster versions of built-in par, sum, prod
when the repeated expression doesn't depend on the variable.
*/

//-------------`(si.)bpar`---------------
// Balanced `par` where the repeated expression doesn't depend on a variable.
// The built-in `par` is implemented as an unbalanced tree, and also has
// to substitute the variable into the repeated expression, which is expensive
// even when the variable doesn't appear. This version is implemented as a
// balanced tree (which allows node reuse during tree traversal) and also
// doesn't search for the variable. This can be much faster than `par` to compile.
//
// #### Usage
//
// ```
// si.bus(N * inputs(f)) : bpar(N, f) : si.bus(N * outputs(f))
// ```
//
// Where:
//
// * `N`: number of repetitions, minimum 1, a constant numerical expression
// * `f`: an arbitrary expression
//
// Example:
// ```
// // square each of 4000 inputs
// process = si.bpar(4000, (_ <: _, _ : *));
// ```
//
//------------------------------------------
declare bpar author "Haggai Nuchi";
declare bpar copyright "Copyright (C) 2024 Haggai Nuchi <haggai@haggainuchi.com>";
declare bpar license "MIT License";
bpar(1, f) = f;
bpar(N, f) = bpar(left, f), bpar(right, f) with {
    left = ma.np2(N) >> 1; // biggest power of 2 less than N
    right = N - left;
};

//-------------`(si.)bsum`---------------
// Balanced `sum`, see `si.bpar`.
//
// #### Usage
//
// ```
// si.bus(N * inputs(f)) : bsum(N, f) : _
// ```
//
// Where:
//
// * `N`: number of repetitions, minimum 1, a constant numerical expression
// * `f`: an arbitrary expression with 1 output.
//
// Example:
// ```
// // square each of 1000 inputs and add the results
// process = si.bsum(1000, (_ <: _, _ : *));
// ```
//----------------------------------------------------
declare bsum author "Haggai Nuchi";
declare bsum copyright "Copyright (C) 2024 Haggai Nuchi <haggai@haggainuchi.com>";
declare bsum license "MIT License";
bsum(1, f) = f;
bsum(N, f) = bsum(left, f) + bsum(right, f) with {
    left = ma.np2(N) >> 1; // biggest power of 2 less than N
    right = N - left;
};

//-------------`(si.)bprod`---------------
// Balanced `prod`, see `si.bpar`.
//
// #### Usage
//
// ```
// si.bus(N * inputs(f)) : bprod(N, f) : _
// ```
//
// Where:
//
// * `N`: number of repetitions, minimum 1, a constant numerical expression
// * `f`: an arbitrary expression with 1 output.
//
// Example:
// ```
// // Add 8000 consecutive inputs (in pairs) and multiply the results
// process = si.bprod(4000, +);
// ```
//----------------------------------------------------
declare bprod author "Haggai Nuchi";
declare bprod copyright "Copyright (C) 2024 Haggai Nuchi <haggai@haggainuchi.com>";
declare bprod license "MIT License";
bprod(1, f) = f;
bprod(N, f) = bprod(left, f) * bprod(right, f) with {
    left = ma.np2(N) >> 1; // biggest power of 2 less than N
    right = N - left;
};