//#################################### interpolators.lib ########################################
// A library to handle interpolation. Its official prefix is `it`.
//
// This library provides several basic interpolation functions as well as interpolators 
// taking a `gen` circuit of N outputs producing values to be interpolated, triggered
// by a `idv` read index signal. Two points and four points interpolations are implemented.
//
// The `idv` parameter is to be used as a read index. In float (= singleprecision) mode, 
// a technique based on 2 signals with the pure integer index and a fractional part in the [0,1] 
// range is used to avoid accumulating errors. In -double (= doubleprecision)/-quad (= quadprecision) modes, 
// a standard implementation with a single fractional index signal is used. 
// Three functions `int_part`, `frac_part` and `mak_idv` are available to manipule the read index signal.
//
// Use-case with `waveform`. Here the signal given to `interpolator_XXX` uses the `idv` model.
// 
// ```
// waveform_interpolator(wf, step, interp) = interp(gen, idv)
// with {
//    gen(idx) = wf, (idx:max(0):min(size-1)) : rdtable with { size = wf:(_,!); };   /* waveform size */
//    index = (+(step)~_)-step;  /* starting from 0 */
//    idv = it.make_idv(index);  /* build the signal for interpolation in a generic way */
// };
// 
// waveform_linear(wf, step) = waveform_interpolator(wf, step, it.interpolator_linear);
// waveform_cosine(wf, step) = waveform_interpolator(wf, step, it.interpolator_cosine);
// waveform_cubic(wf, step) = waveform_interpolator(wf, step, it.interpolator_cubic);
// 
// waveform_interp(wf, step, selector) = waveform_interpolator(wf, step, interp_select(selector))
// with {
//    /* adapts the argument order */
//    interp_select(sel, gen, idv) = it.interpolator_select(gen, idv, sel);
// };
// 
// waveform and index 
// waveform_interpolator1(wf, idv, interp) = interp(gen, idv)
// with {
//    gen(idx) = wf, (idx:max(0):min(size-1)) : rdtable with { size = wf:(_,!); };   /* waveform size */
// };
// 
// waveform_linear1(wf, idv) = waveform_interpolator1(wf, idv, it.interpolator_linear);
// waveform_cosine1(wf, idv) = waveform_interpolator1(wf, idv, it.interpolator_cosine);
// waveform_cubic1(wf, idv) = waveform_interpolator1(wf, idv, it.interpolator_cubic);
// 
// waveform_interp1(wf, idv, selector) = waveform_interpolator1(wf, idv, interp_select(selector))
// with {
//    /* adapts the argument order */
//    interp_select(sel, gen, idv) = it.interpolator_select(gen, idv, sel);
// };
// ```
// 
// Some tests here:
// 
// ```
// wf = waveform {0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 50.0, 40.0, 30.0, 20.0, 10.0, 0.0};
// 
// process = waveform_linear(wf, step), waveform_cosine(wf, step), waveform_cubic(wf, step) with { step = 0.25; };
// 
// process = waveform_interp(wf, 0.25, nentry("algo", 0, 0, 3, 1));
// 
// process = waveform_interp1(wf, idv, nentry("algo", 0, 0, 3, 1))
// with {
//    step = 0.1;
//    idv_aux = (+(step)~_)-step;  /* starting from 0 */
//    idv = it.make_idv(idv_aux);  /* build the signal for interpolation in a generic way */
// };
// 
// /* Test linear interpolation between 2 samples with a `(idx,dv)` signal built using a waveform */
// linear_test = (idx,dv), it.interpolator_linear(gen, (idx,dv))
// with {
//    /* signal to interpolate (only 2 points here) */
//    gen(id) = waveform {3.0, -1.0}, (id:max(0)) : rdtable;
//    dv = waveform {0.0, 0.25, 0.50, 0.75, 1.0}, index : rdtable;
//    idx = 0; 
//    /* test index signal */
//    index = (+(1)~_)-1;   /* starting from 0 */
// };
// 
// /* Test cosine interpolation between 2 samples with a `(idx,dv)` signal built using a waveform */
// cosine_test = (idx,dv), it.interpolator_cosine(gen, (idx,dv))
// with {
//    /* signal to interpolate (only 2 points here) */
//    gen(id) = waveform {3.0, -1.0}, (id:max(0)) : rdtable;
//    dv = waveform {0.0, 0.25, 0.50, 0.75, 1.0}, index : rdtable;
//    idx = 0;
//    /* test index signal */
//    index = (+(1)~_)-1;   /* starting from 0 */
// };
// 
// /* Test cubic interpolation between 4 samples with a `(idx,dv)` signal built using a waveform */
// cubic_test = (idx,dv), it.interpolator_cubic(gen, (idx,dv))
// with {
//    /* signal to interpolate (only 4 points here) */
//    gen(id) = waveform {-1.0, 2.0, 1.0, 4.0}, (id:max(0)) : rdtable;
//    dv = waveform {0.0, 0.25, 0.50, 0.75, 1.0}, index : rdtable;
//    idx = 0;
//    /* test index signal */
//    index = (+(1)~_)-1;   /* starting from 0 */
// };
// ```
//########################################################################################

/************************************************************************
************************************************************************
FAUST library file
Copyright (C) 2019-2020 GRAME, Centre National de Creation Musicale
----------------------------------------------------------------------
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.
************************************************************************
************************************************************************/

declare name "Faust Interpolator Library";
declare version "0.3";

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

reset(trig) = (trig-trig') <= 0;
 
// The following 3 functions allow to adapt a 'single signal' fractional idv in this [idx, dv] model.

singleprecision int_part(idv) = idv : _,!;
singleprecision frac_part(idv) = idv : !,_;
singleprecision make_idv(id) = int(id), ma.frac(id);

// Infinite raising index
singleprecision raise(trig, step, length) = id, dv
letrec {
    'id = (id + int(step) + int(dv + ma.frac(step))) * reset(trig);
    'dv = ma.frac(dv + ma.frac(step)) * reset(trig);
};

// Modulo based raising index
singleprecision raise_modulo(trig, step, length) = id, dv
letrec {
    'id = fmod(id + int(step) + int(dv + ma.frac(step)), length) * reset(trig);
    'dv = ma.frac(dv + ma.frac(step)) * reset(trig);
};

// Decreasing index starting at 'length'
singleprecision decrease(trig, step, length) = raise(trig, -step, length) : (+(length), _);

// Modulo decreasing index starting at 'length'
singleprecision decrease_modulo(trig, step, length) = raise_modulo(trig, -step, length) : (+(length), _);

// The following 3 functions allow to adapt a 'single signal' fractional idv in this [idx, dv] model.

doubleprecision quadprecision int_part(idv) = int(idv);
doubleprecision quadprecision frac_part(idv) = ma.frac(idv);
doubleprecision quadprecision make_idv(id) = id;

// Infinite raising index
doubleprecision quadprecision raise(trig, step, length) = idv
letrec {
    'idv = (idv + step) * reset(trig);
};

// Modulo based raising index
doubleprecision quadprecision raise_modulo(trig, step, length) = idv
letrec {
    'idv = fmod(idv + step, length) * reset(trig);
};

// Decreasing index starting at 'length'
doubleprecision quadprecision decrease(trig, step, length) = raise(trig, -step, length) + length;

// Modulo decreasing index starting at 'length'
doubleprecision quadprecision decrease_modulo(trig, step, length) = raise_modulo(trig, -step, length) + length;
 
//=========================Two points interpolation functions=============================
//========================================================================================

//-------`(it.)interpolate_linear`----------
// Linear interpolation between 2 values.
//
// #### Usage
//
// ```
// interpolate_linear(dv,v0,v1) : _
// ```
//
// Where:
//
// * `dv`: in the fractional value in [0..1] range
// * `v0`: is the first value
// * `v1`: is the second value
//
//
// #### Reference:
//
// <https://github.com/jamoma/JamomaCore/blob/master/Foundation/library/includes/TTInterpolate.h>
//
//--------------------------------------------

interpolate_linear(dv,v0,v1) = v0 + dv*(v1-v0);  // (faster than v0*(1-dv)+v1*dv which is currently not optimized...)


//-------`(it.)interpolate_cosine`----------
// Cosine interpolation between 2 values.
//
// #### Usage
//
// ```
// interpolate_cosine(dv,v0,v1) : _
// ```
//
// Where:
//
// * `dv`: in the fractional value in [0..1] range
// * `v0`: is the first value
// * `v1`: is the second value
//
//
// #### Reference:
//
// <https://github.com/jamoma/JamomaCore/blob/master/Foundation/library/includes/TTInterpolate.h>
//
//--------------------------------------------

interpolate_cosine(dv,v0,v1) = v0 + a2*(v1-v0) with { a2 = 0.5 * (1.0 - cos(dv*ma.PI)); };


//=========================Four points interpolation functions============================
//========================================================================================

//-------`(it.)interpolate_cubic`----------
// Cubic interpolation between 4 values.
//
// #### Usage
//
// ```
// interpolate_cubic(dv,v0,v1,v2,v3) : _
// ```
//
// Where:
//
// * `dv`: in the fractional value in [0..1] range
// * `v0`: is the first value
// * `v1`: is the second value
// * `v2`: is the third value
// * `v3`: is the fourth value
//
//
// #### Reference:
//
// <https://www.paulinternet.nl/?page=bicubic>
//
//--------------------------------------------

interpolate_cubic(dv,v0,v1,v2,v3) 
	= v1 + 0.5*dv*(v2 - v0 + dv*(2.0*v0 - 5.0*v1 + 4.0*v2 - v3 + dv*(3.0*(v1 - v2) + v3 - v0)));


//=========================Two points interpolators=======================================
//========================================================================================

//-------`(it.)interpolator_two_points`----------
// Generic interpolator on two points (current and next index), assuming an increasing index.
//
// #### Usage
//
// ```
// interpolator_two_points(gen, idv, interpolate_two_points) : _,_... (equal to N = outputs(gen))
// ```
//
// Where:
//
// * `gen`: a circuit with an 'idv' reader input that produces N outputs
// * `idv`: a fractional read index expressed as a float value, or a (int,frac) pair
// * `interpolate_two_points`: a two points interpolation function
//
//--------------------------------------------

interpolator_two_points(gen, idv, interpolate_two_points) = (gen(id0), gen(id1))
                                                            : ro.interleave(outputs(gen), 2)
                                                            : par(i, outputs(gen), interpolate_two_points(dv))
with {
    id0 = int_part(idv); // index integer part
    id1 = id0 + 1;       // next index
    dv = frac_part(idv); // index fractional part in [0..1]
};


//-------`(it.)interpolator_linear`----------
// Linear interpolator for a 'gen' circuit triggered by an 'idv' input to generate values.
//
// #### Usage
//
// ```
// interpolator_linear(gen, idv) : _,_... (equal to N = outputs(gen))
// ```
//
// Where:
//
// * `gen`: a circuit with an 'idv' reader input that produces N outputs
// * `idv`: a fractional read index expressed as a float value, or a (int,frac) pair
//
//--------------------------------------------

interpolator_linear(gen, idv) = interpolator_two_points(gen, idv, interpolate_linear);


//-------`(it.)interpolator_cosine`----------
// Cosine interpolator for a 'gen' circuit triggered by an 'idv' input to generate values.
//
// #### Usage
//
// ```
// interpolator_cosine(gen, idv) : _,_... (equal to N = outputs(gen))
// ```
//
// Where:
//
// * `gen`: a circuit with an 'idv' reader input that produces N outputs
// * `idv`: a fractional read index expressed as a float value, or a (int,frac) pair
//
//--------------------------------------------

interpolator_cosine(gen, idv) = interpolator_two_points(gen, idv, interpolate_cosine);


// To be used in 'interpolator_select'
interpolator_null(gen, idv) = interpolator_two_points(gen, idv, \(dv,v0,v1).(v0));


//=========================Four points interpolators======================================
//========================================================================================

//-------`(it.)interpolator_two_points`----------
// Generic interpolator on interpolator_four_points points (previous, current and two next indexes), assuming an increasing index.
//
// #### Usage
//
// ```
// interpolator_four_points(gen, idv, interpolate_four_points) : _,_... (equal to N = outputs(gen))
// ```
//
// Where:
//
// * `gen`: a circuit with an 'idv' reader input that produces N outputs
// * `idv`: a fractional read index expressed as a float value, or a (int,frac) pair
// * `interpolate_four_points`: a four points interpolation function
//
//--------------------------------------------

interpolator_four_points(gen, idv, interpolate_four_points) = (gen(id0), gen(id1), gen(id2), gen(id3))
                                                            : ro.interleave(outputs(gen), 4)
                                                            : par(i, outputs(gen), interpolate_four_points(dv))
with {
    id0 = id1 - 1;       // previous index
    id1 = int_part(idv); // index integer part
    id2 = id1 + 1;       // next index
    id3 = id2 + 1;       // next index
    dv = frac_part(idv); // index fractional part in [0..1]
};


//-------`(it.)interpolator_cubic`----------
// Cubic interpolator for a 'gen' circuit triggered by an 'idv' input to generate values
//
// #### Usage
//
// ```
// interpolator_cubic(gen, idv) : _,_... (equal to N = outputs(gen))
// ```
//
// Where:
//
// * `gen`: a circuit with an 'idv' reader input that produces N outputs
// * `idv`: a fractional read index expressed as a float value, or a (int,frac) pair
//
//--------------------------------------------

interpolator_cubic(gen, idv) = interpolator_four_points(gen, idv, interpolate_cubic);

// Enum of interpolation algorithms
MAX_INTER = 4;

linear = 0;
cosine = 1;
cubic  = 2;
nointerp = MAX_INTER-1;


//-------`(it.)interpolator_select`----------
// Generic configurable interpolator (with selector between in [0..3]). The value 3 is used for no interpolation.
//
// #### Usage
//
// ```
// interpolator_select(gen, idv, sel) : _,_... (equal to N = outputs(gen))
// ```
//
// Where:
//
// * `gen`: a circuit with an 'idv' reader input that produces N outputs
// * `idv`: a fractional read index expressed as a float value, or a (int,frac) pair
// * `sel`: an interpolation algorithm selector in [0..3] (0 = linear, 1 = cosine, 2 = cubic, 3 = nointerp)
//
//--------------------------------------------

interpolator_select(gen, idv, sel) = ba.selectmulti(ma.SR/10, interpolators, sel)
with {
	interpolators = (interpolator_linear(gen, idv),
                    interpolator_cosine(gen, idv),
                    interpolator_cubic(gen, idv),
                    interpolator_null(gen, idv));
};