//################################## misceffects.lib ##########################################
// Collection of audio effects library. Its official prefix is `ef`.
//
// The library is organized into 7 sections:
//
// * [Dynamic](#Dynamic)
// * [Fibonacci](#fibonacci)
// * [Filtering](#filtering)
// * [Meshes](#meshes)
// * [Mixing](#mixing)
// * [Time Based](#time-based)
// * [Pitch Shifting](#pitch-shifting)
// * [Saturators](#saturators)
//
// #### References
// * <https://github.com/grame-cncm/faustlibraries/blob/master/misceffects.lib>
//########################################################################################

/************************************************************************
************************************************************************
FAUST library file
Copyright (C) 2003-2016 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.
************************************************************************
************************************************************************/

ma = library("maths.lib");
ba = library("basics.lib");
de = library("delays.lib");
si = library("signals.lib");
an = library("analyzers.lib");
fi = library("filters.lib");
ro = library("routes.lib");
aa = library("aanl.lib");
ef = library("misceffects.lib"); // for compatible copy/paste out of this file

declare name "Misc Effects Library";
declare version "2.5.1"; 

//======================================Dynamic===========================================
//========================================================================================

//---------------------`(ef.)cubicnl`-----------------------
// Cubic nonlinearity distortion.
// `cubicnl` is a standard Faust function.
//
// #### Usage:
//
// ```
// _ : cubicnl(drive,offset) : _
// _ : cubicnl_nodc(drive,offset) : _
// ```
//
// Where:
//
// * `drive`: distortion amount, between 0 and 1
// * `offset`: constant added before nonlinearity to give even harmonics. Note: offset
//    can introduce a nonzero mean - feed cubicnl output to dcblocker to remove this.
//
// #### References:
//
// * <https://ccrma.stanford.edu/~jos/pasp/Cubic_Soft_Clipper.html>
// * <https://ccrma.stanford.edu/~jos/pasp/Nonlinear_Distortion.html>
//------------------------------------------------------------
cubicnl(drive,offset) = *(pregain) : +(offset) : clip(-1,1) : cubic
with {
    pregain = pow(10.0,2*drive);
    clip(lo,hi) = min(hi) : max(lo);
    cubic(x) = x - x*x*x/3;
    postgain = max(1.0,1.0/pregain);
};

cubicnl_nodc(drive,offset) = cubicnl(drive,offset) : fi.dcblocker;

declare cubicnl author "Julius O. Smith III";
declare cubicnl license "STK-4.3";
declare cubicnl_nodc author "Julius O. Smith III";
declare cubicnl_nodc license "STK-4.3";

//-----------------`(ef.)gate_mono`-------------------
// Mono signal gate.
// `gate_mono` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : gate_mono(thresh,att,hold,rel) : _
// ```
//
// Where:
//
// * `thresh`: dB level threshold above which gate opens (e.g., -60 dB)
// * `att`: attack time = time constant (sec) for gate to open (e.g., 0.0001 s = 0.1 ms)
// * `hold`: hold time = time (sec) gate stays open after signal level < thresh (e.g., 0.1 s)
// * `rel`: release time = time constant (sec) for gate to close (e.g., 0.020 s = 20 ms)
//
// #### References
//
// * <http://en.wikipedia.org/wiki/Noise_gate>
// * <http://www.soundonsound.com/sos/apr01/articles/advanced.asp>
// * <http://en.wikipedia.org/wiki/Gating_(sound_engineering)>
//------------------------------------------------------------
gate_mono(thresh,att,hold,rel,x) = x * gate_gain_mono(thresh,att,hold,rel,x);

declare gate_mono author "Julius O. Smith III";
declare gate_mono license "STK-4.3";

//-----------------`(ef.)gate_stereo`-------------------
// Stereo signal gates.
// `gate_stereo` is a standard Faust function.
//
// #### Usage
//
// ```
//  _,_ : gate_stereo(thresh,att,hold,rel) : _,_
// ```
//
// Where:
//
// * `thresh`: dB level threshold above which gate opens (e.g., -60 dB)
// * `att`: attack time = time constant (sec) for gate to open (e.g., 0.0001 s = 0.1 ms)
// * `hold`: hold time = time (sec) gate stays open after signal level < thresh (e.g., 0.1 s)
// * `rel`: release time = time constant (sec) for gate to close (e.g., 0.020 s = 20 ms)
//
// #### References
//
// * <http://en.wikipedia.org/wiki/Noise_gate>
// * <http://www.soundonsound.com/sos/apr01/articles/advanced.asp>
// * <http://en.wikipedia.org/wiki/Gating_(sound_engineering)>
//------------------------------------------------------------
gate_stereo(thresh,att,hold,rel,x,y) = ggm*x, ggm*y with {
  ggm = gate_gain_mono(thresh,att,hold,rel,abs(x)+abs(y));
};

gate_gain_mono(thresh,att,hold,rel,x) = x : extendedrawgate : an.amp_follower_ar(att,rel) with {
  extendedrawgate(x) = max(float(rawgatesig(x)),holdsig(x));
  rawgatesig(x) = inlevel(x) > ba.db2linear(thresh);
  minrate = min(att,rel);
  inlevel = an.amp_follower_ar(minrate,minrate);
  holdcounter(x) = (max(holdreset(x) * holdsamps,_) ~-(1));
  holdsig(x) = holdcounter(x) > 0;
  holdreset(x) = rawgatesig(x) < rawgatesig(x)'; // reset hold when raw gate falls
  holdsamps = int(hold*ma.SR);
};

declare gate_stereo author "Julius O. Smith III";
declare gate_stereo license "STK-4.3";
declare gate_gain_mono author "Julius O. Smith III";
declare gate_gain_mono license "STK-4.3";

//=====================================Fibonacci==========================================
//========================================================================================

//---------------`(ef.)fibonacci`---------------------------
// Fibonacci system where the current output is the current
// input plus the sum of the previous N outputs.
//
// #### Usage
//
// ```
// _ : fibonacci(N) : _
// ```
//
// Where:
//
// * `N`: the Fibonacci system's order, where 2 is standard
//
// #### Example
// Generate the famous series: [1, 1, 2, 3, 5, 8, 13, ...]
//
// ```
// 1. : ba.impulsify : fibonacci(2)
// ```
//------------------------------------------------------------
declare fibonacci author "Dario Sanfilippo";

fibonacci(order) = +~(_<:sum(i, order, @(i)):>_);


//---------------`(ef.)fibonacciGeneral`----------------------
// Fibonacci system with customizable coefficients.
// The order of the system is inferred from the number of coefficients.
//
// #### Usage
//
// ```
// _ : fibonacciGeneral(wave) : _
// ```
//
// Where:
//
// * `wave`: a waveform such as `waveform{1, 1}`
//
// #### Example:
// Use the update equation `y = 2*y' + 3*y'' + 4*y'''`
//
// ```
// 1. : ba.impulsify : fibonacciGeneral(waveform{2, 3, 4})
// ```
//------------------------------------------------------------
declare fibonacciGeneral author "Dario Sanfilippo and David Braun";

fibonacciGeneral(wave) = +~(_<:sum(i, N, func(i)):>_)
with {
    N = wave : _, !;
    func(i) = @(i) : _ * (wave, i : rdtable);
};

//---------------`(ef.)fibonacciSeq`---------------------------
// First N numbers of the Fibonacci sequence [1, 1, 2, 3, 5, 8, ...]
// as parallel channels.
//
// #### Usage
//
// ```
// fibonacciSeq(N) : si.bus(N)
// ```
//
// Where:
//
// * `N`: The number of Fibonacci numbers to generate as channels.
// 
//------------------------------------------------------------

fibonacciSeq(N) = iterate(N, (1, 1))
with {
  iterate(1, (A0, A1)) = A0;
  iterate(N, (A0, A1)) = A0 , iterate(N - 1, (A1, A0 + A1));
};

declare fibonacciSeq author "Dario Sanfilippo";


//=====================================Filtering==========================================
//========================================================================================

//-------------------------`(ef.)speakerbp`-------------------------------
// Dirt-simple speaker simulator (overall bandpass eq with observed
// roll-offs above and below the passband). `speakerbp` is a standard Faust function.
//
// Low-frequency speaker model = +12 dB/octave slope breaking to
// flat near f1. Implemented using two dc blockers in series.
//
// High-frequency model = -24 dB/octave slope implemented using a
// fourth-order Butterworth lowpass.

//
//
// #### Usage
// ```
// _ : speakerbp(f1,f2) : _
// ```
// #### Example
//
// Based on measured Celestion G12 (12" speaker):
// ```
// speakerbp(130,5000)
// ```
//------------------------------------------------------------
// TODO: perhaps this should be moved to physmodels.lib
// [JOS: I don't think so because it's merely a bandpass filter tuned to speaker bandwidth]
speakerbp(f1,f2) = fi.dcblockerat(f1) : fi.dcblockerat(f1) : fi.lowpass(4,f2);

declare speakerbp author "Julius O. Smith III";
declare speakerbp license "STK-4.3";

//------------`(ef.)piano_dispersion_filter`---------------
// Piano dispersion allpass filter in closed form.
//
// #### Usage
//
// ```
// piano_dispersion_filter(M,B,f0)
// _ : piano_dispersion_filter(1,B,f0) : +(totalDelay),_ : fdelay(maxDelay) : _
// ```
//
// Where:
//
// * `M`: number of first-order allpass sections (compile-time only)
//    Keep below 20. 8 is typical for medium-sized piano strings.
// * `B`: string inharmonicity coefficient (0.0001 is typical)
// * `f0`: fundamental frequency in Hz
//
// #### Outputs
//
// * MINUS the estimated delay at `f0` of allpass chain in samples,
//     provided in negative form to facilitate subtraction
//     from delay-line length.
// * Output signal from allpass chain
//
// #### Reference
//
// * "Dispersion Modeling in Waveguide Piano Synthesis Using Tunable
// Allpass Filters", by Jukka Rauhala and Vesa Valimaki, DAFX-2006, pp. 71-76
// * <http://lib.tkk.fi/Diss/2007/isbn9789512290666/article2.pdf>
//   An erratum in Eq. (7) is corrected in Dr. Rauhala's encompassing
//   dissertation (and below).
// * <http://www.acoustics.hut.fi/research/asp/piano/>
//------------------------------------------------------------
// TODO: perhaps this should be moved to physmodels.lib?
// [JOS: I vote yes when there is a piano model in physmodels.lib.]

piano_dispersion_filter(M,B,f0) = -Df0*M,seq(i,M,fi.tf1(a1,1,a1))
with {
 a1 = (1-D)/(1+D); // By Eq. 3, have D >= 0, hence a1 >= 0 also
 D = exp(Cd - Ikey(f0)*kd);
 trt = pow(2.0,1.0/12.0); // 12th root of 2
 logb(b,x) = log(x) / log(b); // log-base-b of x
 Ikey(f0) = logb(trt,f0*trt/27.5);
 Bc = max(B,0.000001);
 kd = exp(k1*log(Bc)*log(Bc) + k2*log(Bc)+k3);
 Cd = exp((m1*log(M)+m2)*log(Bc)+m3*log(M)+m4);
 k1 = -0.00179;
 k2 = -0.0233;
 k3 = -2.93;
 m1 = 0.0126;
 m2 = 0.0606;
 m3 = -0.00825;
 m4 = 1.97;
 wT = 2*ma.PI*f0/ma.SR;
 polydel(a) = atan(sin(wT)/(a+cos(wT)))/wT;
 Df0 = polydel(a1) - polydel(1.0/a1);
};

declare piano_dispersion_filter author "Julius O. Smith III";
declare piano_dispersion_filter license "STK-4.3";

//-------------------------`(ef.)stereo_width`---------------------------
// Stereo Width effect using the Blumlein Shuffler technique.
// `stereo_width` is a standard Faust function.
//
// #### Usage
//
// ```
// _,_ : stereo_width(w) : _,_
// ```
//
// Where:
//
// * `w`: stereo width between 0 and 1
//
// At `w=0`, the output signal is mono ((left+right)/2 in both channels).
// At `w=1`, there is no effect (original stereo image).
// Thus, w between 0 and 1 varies stereo width from 0 to "original".
//
// #### Reference
//
// * "Applications of Blumlein Shuffling to Stereo Microphone Techniques"
// Michael A. Gerzon, JAES vol. 42, no. 6, June 1994
//------------------------------------------------------------
stereo_width(w) = shuffle : *(mgain),*(sgain) : shuffle
with {
     shuffle =  _,_ <: +,-; // normally scaled by 1/sqrt(2) for orthonormality,
     mgain = 1-w/2;  // but we pick up the needed normalization here.
     sgain = w/2;
};

declare stereo_width author "Julius O. Smith III";
declare stereo_width license "STK-4.3";

//===========================================Meshes=======================================
//========================================================================================


// TODO: the following should be in physmodels.lib when it will be operational
// [JOS: I think a new "Meshes" section would fit well after Modal Percussions.]
//----------------------------------`(ef.)mesh_square`------------------------------
// Square Rectangular Digital Waveguide Mesh.
//
// #### Usage
//
// ```
// bus(4*N) : mesh_square(N) : bus(4*N)
// ```
//
// Where:
//
// * `N`: number of nodes along each edge - a power of two (1,2,4,8,...)
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Mesh.html>
//
// #### Signal Order In and Out
//
// The mesh is constructed recursively using 2x2 embeddings. Thus,
// the top level of `mesh_square(M)` is a block 2x2 mesh, where each
// block is a `mesh(M/2)`. Let these blocks be numbered 1,2,3,4 in the
// geometry NW,NE,SW,SE, i.e., as:
//
//         1 2
//         3 4
// 
// Each block has four vector inputs and four vector outputs, where the
// length of each vector is `M/2`. Label the input vectors as Ni,Ei,Wi,Si,
// i.e., as the inputs from the North, East South, and West,
// and similarly for the outputs. Then, for example, the upper
// left input block of M/2 signals is labeled 1Ni. Most of the
// connections are internal, such as 1Eo -> 2Wi. The `8*(M/2)` input
// signals are grouped in the order:
// 
//        1Ni 2Ni
//        3Si 4Si
//        1Wi 3Wi
//        2Ei 4Ei
// 
// and the output signals are:
//
//        1No 1Wo
//        2No 2Eo
//        3So 3Wo
//        4So 4Eo
// or:
//
//        In: 1No 1Wo 2No 2Eo 3So 3Wo 4So 4Eo
//        Out: 1Ni 2Ni 3Si 4Si 1Wi 3Wi 2Ei 4Ei
//
// Thus, the inputs are grouped by direction N,S,W,E, while the
// outputs are grouped by block number 1,2,3,4, which can also be
// interpreted as directions NW, NE, SW, SE.  A simple program
// illustrating these orderings is `process = mesh_square(2);`.
//
// #### Example
//
// Reflectively terminated mesh impulsed at one corner:
//
// ```
// mesh_square_test(N,x) = mesh_square(N)~(busi(4*N,x)) // input to corner
// with { 
//     busi(N,x) = bus(N) : par(i,N,*(-1)) : par(i,N-1,_), +(x); 
// };
// process = 1-1' : mesh_square_test(4); // all modes excited forever
// ```
//
// In this simple example, the mesh edges are connected as follows:
//
//        1No -> 1Ni, 1Wo -> 2Ni, 2No -> 3Si, 2Eo -> 4Si,
//        3So -> 1Wi, 3Wo -> 3Wi, 4So -> 2Ei, 4Eo -> 4Ei
//
// A routing matrix can be used to obtain other connection geometries.
//------------------------------------------------------------
// four-port scattering junction:
mesh_square(1) =
          si.bus(4) <: par(i,4,*(-1)), (si.bus(4) :> (*(.5)) <: si.bus(4)) :> si.bus(4);

// rectangular NxN square waveguide mesh:
mesh_square(N) = si.bus(4*N) : (route_inputs(N/2) : par(i,4,mesh_square(N/2)))
          ~(prune_feedback(N/2))
          : prune_outputs(N/2) : route_outputs(N/2) : si.bus(4*N)
with {
  // select block i of N, block size = M:
  s(i,N,M) = par(j, M*N, Sv(i, j))
     with { Sv(i,i) = si.bus(N); Sv(i,j) = si.block(N); };

  // prune mesh outputs down to the signals which make it out:
  prune_outputs(N)
    = si.bus(16*N) :
      si.block(N), si.bus(N),   si.block(N), si.bus(N),
      si.block(N), si.bus(N),   si.bus(N),   si.block(N),
      si.bus(N),   si.block(N), si.block(N), si.bus(N),
      si.bus(N),   si.block(N), si.bus(N),   si.block(N)
      : si.bus(8*N);

  // collect mesh outputs into standard order (N,W,E,S):
  route_outputs(N)
    = si.bus(8*N)
      <: s(4,N,8),s(5,N,8), s(0,N,8),s(2,N,8),
         s(3,N,8),s(7,N,8), s(1,N,8),s(6,N,8)
      : si.bus(8*N);

  // collect signals used as feedback:
  prune_feedback(N) = si.bus(16*N) :
      si.bus(N),   si.block(N), si.bus(N),   si.block(N),
      si.bus(N),   si.block(N), si.block(N), si.bus(N),
      si.block(N), si.bus(N),   si.bus(N),   si.block(N),
      si.block(N), si.bus(N),   si.block(N), si.bus(N) :
      si.bus(8*N);

  // route mesh inputs (feedback, external inputs):
  route_inputs(N) = si.bus(8*N), si.bus(8*N)
  <:s(8,N,16),s(4,N,16), s(12,N,16),s(3,N,16),
    s(9,N,16),s(6,N,16), s(1,N,16),s(14,N,16),
    s(0,N,16),s(10,N,16), s(13,N,16),s(7,N,16),
    s(2,N,16),s(11,N,16), s(5,N,16),s(15,N,16)
    : si.bus(16*N);
};

declare mesh_square author "Julius O. Smith III";
declare mesh_square license "STK-4.3";


//=====================================Mixing=============================================
//========================================================================================

// Implementation to share common code
dwmEnv(wetAmount, FX) = environment
{
    N = inputs(FX);
    wet(wg) = FX : par(i, N, *(wg));
    dry(dg) = par(i, N, *(dg));
    out(wg, dg) = si.bus(N) <: wet(wg), dry(dg) :> si.bus(N);
    
    dryWetMixer = out(wetGain, dryGain) 
    with { 
        wetGain = wetAmount; 
        dryGain = 1. - wetGain; 
    };
    
    dryWetMixerConstantPower = out(wetGain, dryGain)  
    with {
        theta = ma.PI*wetAmount/2.;
        dryGain = cos(theta)/sqrt(2.);
        wetGain = sin(theta)/sqrt(2.);
    };
};


//---------------`(ef.)dryWetMixer`-------------
// Linear dry-wet mixer for a N inputs and N outputs effect.
//
// #### Usage
//
// ```
// si.bus(inputs(FX)) : dryWetMixer(wetAmount, FX) : si.bus(inputs(FX))
// ```
//
// Where:
//
// * `wetAmount`: the wet amount (0-1). 0 produces only the dry signal and 1 produces only the wet signal
// * `FX`: an arbitrary effect (N inputs and N outputs) to apply to the input bus
//------------------------------------------------------------
declare dryWetMixer author "David Braun, revised by Stéphane Letz";

dryWetMixer(wetAmount, FX) = dwmEnv(wetAmount, FX).dryWetMixer;


//---------------`(ef.)dryWetMixerConstantPower`-------------
// Constant-power dry-wet mixer for a N inputs and N outputs effect.
//
// #### Usage
//
// ```
// si.bus(inputs(FX)) : dryWetMixerConstantPower(wetAmount, FX) :si.bus(inputs(FX))
// ```
//
// Where:
//
// * `wetAmount`: the wet amount (0-1). 0 produces only the dry signal and 1 produces only the wet signal
// * `FX`: an arbitrary effect (N inputs and N outputs) to apply to the input bus
//------------------------------------------------------------
declare dryWetMixerConstantPower author "David Braun, revised by Stéphane Letz";

dryWetMixerConstantPower(wetAmount, FX) = dwmEnv(wetAmount, FX).dryWetMixerConstantPower;


mixingEnv = environment
{
  // Note that i goes from 0 to N-1.
  // m goes from 0 to N-1 typically, but the output should be periodic with size N.
  // In other words the output with m=-4*N is the same as -2*N, -1*N, 0, 1*N, 2*N etc.
  phaseLoop(N, m, i) = select2(abs(phase1)<abs(phase2), phase2, phase1)
  with {
    phase1 = fmod(i-m,N);
    phase2 = phase1+ba.if(phase1<0,N,-N);
  };

  phaseClamp(N, m, i) = i-aa.clip(0,N-1,m);

  // We divide by sqrt(2) at the end so that for m=0.5,1.5,2.5 etc,
  // the total gain is 1.0, matching phase2LinearWeight. However,
  // this means for m=0,1,2,3, etc, the gain is (1./sqrt(2)~=0.7071).
  phase2PowerWeight = aa.clip(-1, 1) : cos(_*ma.PI*.5) / sqrt(2.);

  phase2LinearWeight = aa.clip(-1, 1) : 1-abs(_);

  //------------------------`weightsPowerLoop`---------------------------
  // "Fan out" an index into N weights between 0 and 1. At any given
  // moment, two weights may be non-zero. Suppose they are N_m and N_{m+1}.
  // Then `cos(N_m)^2+sin(N_{m+1})^2==0.5`. 
  //
  // #### Usage
  //
  // ```
  // _ : weightsPowerLoop(N) : si.bus(N)
  // ```
  // Where:
  //
  // * `N`: number of output weights
  // * `m`: [0;N-1] (float) blend index. If m is outside [0;N-1], the behavior will loop.
  //.       So m=-N, m=0, and m=N should give the same output.
  weightsPowerLoop(N, m) = par(i, N, gain(i))
  with {
    gain(i) = phaseLoop(N, m, i) : phase2PowerWeight;
  };

  // Same as above, but the two weights being blended at any moment SUM to 1.
  weightsLinearLoop(N, m) = par(i, N, gain(i))
  with {
    gain(i) = phaseLoop(N, m, i) : phase2LinearWeight;
  };

  // Same as weightsPowerLoop, but m is clamped to [0;N-1]
  weightsPowerClamp(N, m) = par(i, N, gain(i))
  with {
    gain(i) = phaseClamp(N, m, i) : phase2PowerWeight;
  };

  // Same as weightsLinearLoop, but m is clamped to [0;N-1]
  weightsLinearClamp(N, m) = par(i, N, gain(i))
  with {
    gain(i) = phaseClamp(N, m, i) : phase2LinearWeight;
  };

  dryWetMixer(wetAmount, FX) = si.vecOp((weights, sounds), *) :> si.bus(C)
  with {
    N = 2; // We know in advance that there are 2 sounds (the dry and wet).
    C = inputs(FX);
    weights = weightsLinearClamp(N, wetAmount) <: ro.interleave(N, C);
    sounds = si.bus(C) <: si.bus(C), FX;
  };

  dryWetMixerConstantPower(wetAmount, FX) = si.vecOp((weights, sounds), *) :> si.bus(C)
  with {
    N = 2; // We know in advance that there are 2 sounds (the dry and wet).
    C = inputs(FX);
    weights = weightsPowerClamp(N, wetAmount) <: ro.interleave(N, C);
    sounds = si.bus(C) <: si.bus(C), FX;
  };

  // Suppose `sounds` is N buses, each of C channels.
  // We want to linearly mix the buses using index `m` [0;N-1]
  mixLinearClamp(N, C, m) = si.vecOp((weights, si.bus(N)), *) :> si.bus(C)
  with {
    weights = weightsLinearClamp(N, m) <: ro.interleave(N, C);
  };

  mixLinearLoop(N, C, m) = si.vecOp((weights, si.bus(N)), *) :> si.bus(C)
  with {
    weights = weightsLinearLoop(N, m) <: ro.interleave(N, C);
  };

  mixPowerClamp(N, C, m) = si.vecOp((weights, si.bus(N)), *) :> si.bus(C)
  with {
    weights = weightsPowerClamp(N, m) <: ro.interleave(N, C);
  };

  mixPowerLoop(N, C, m) = si.vecOp((weights, si.bus(N)), *) :> si.bus(C)
  with {
    weights = weightsPowerLoop(N, m) <: ro.interleave(N, C);
  };
};


//---------------`(ef.)mixLinearClamp`-------------------------------------------------
// Linear mixer for `N` buses, each with `C` channels. The output will be a sum of 2 buses
// determined by the mixing index `mix`. 0 produces the first bus, 1 produces the
// second, and so on. `mix` is clamped automatically. For example, `mixLinearClamp(4, 1, 1)`
// will weight its 4 inputs by `(0, 1, 0, 0)`. Similarly, `mixLinearClamp(4, 1, 1.1)`
// will weight its 4 inputs by `(0,.9,.1,0)`.
//
// #### Usage
//
// ```
// si.bus(N*C) : mixLinearClamp(N, C, mix) : si.bus(C)
// ```
//
// Where:
//
// * `N`: the number of input buses
// * `C`: the number of channels in each bus
// * `mix`: the mixing index, continuous in [0;N-1].
//---------------------------------------------------------------------------------------
declare mixLinearClamp author "David Braun";

mixLinearClamp = mixingEnv.mixLinearClamp;


//---------------`(ef.)mixLinearLoop`-------------------------------------------------
// Linear mixer for `N` buses, each with `C` channels. Refer to `mixLinearClamp`. `mix`
// will loop for multiples of `N`. For example, `mixLinearLoop(4, 1, 0)` has the same
// effect as `mixLinearLoop(4, 1, -4)` and `mixLinearLoop(4, 1, 4)`.
//
// #### Usage
//
// ```
// si.bus(N*C) : mixLinearLoop(N, C, mix) : si.bus(C)
// ```
//
// Where:
//
// * `N`: the number of input buses
// * `C`: the number of channels in each bus
// * `mix`: the mixing index (N-1) selects the last bus, and 0 or N selects the 0th bus.
//---------------------------------------------------------------------------------------
declare mixLinearLoop author "David Braun";

mixLinearLoop = mixingEnv.mixLinearLoop;


//---------------`(ef.)mixPowerClamp`-------------------------------------------------
// Constant-power mixer for `N` buses, each with `C` channels. The output will be a sum of 2 buses
// determined by the mixing index `mix`. 0 produces the first bus, 1 produces the
// second, and so on. `mix` is clamped automatically. `mixPowerClamp(4, 1, 1)`
// will weight its 4 inputs by `(0, 1./sqrt(2), 0, 0)`. Similarly, `mixPowerClamp(4, 1, 1.5)`
// will weight its 4 inputs by `(0,.5,.5,0)`.
//
// #### Usage
//
// ```
// si.bus(N*C) : mixPowerClamp(N, C, mix) : si.bus(C)
// ```
//
// Where:
//
// * `N`: the number of input buses
// * `C`: the number of channels in each bus
// * `mix`: the mixing index, continuous in [0;N-1].
//---------------------------------------------------------------------------------------
declare mixPowerClamp author "David Braun";

mixPowerClamp = mixingEnv.mixPowerClamp;


//---------------`(ef.)mixPowerLoop`-----------------------------------------------------
// Constant-power mixer for `N` buses, each with `C` channels. Refer to `mixPowerClamp`. `mix`
// will loop for multiples of `N`. For example, `mixPowerLoop(4, 1, 0)` has the same effect
// as `mixPowerLoop(4, 1, -4)` and `mixPowerLoop(4, 1, 4)`.
//
// #### Usage
//
// ```
// si.bus(N*C) : mixPowerLoop(N, C, mix) : si.bus(C)
// ```
//
// Where:
//
// * `N`: the number of input buses
// * `C`: the number of channels in each bus
// * `mix`: the mixing index (N-1) selects the last bus, and 0 or N selects the 0th bus.
//---------------------------------------------------------------------------------------
declare mixPowerLoop author "David Braun";

mixPowerLoop = mixingEnv.mixPowerLoop;


//========================================Time Based======================================
//========================================================================================

//----------`(ef.)echo`----------
// A simple echo effect.
// `echo` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : echo(maxDuration,duration,feedback) : _
// ```
//
// Where:
//
// * `maxDuration`: the max echo duration in seconds
// * `duration`: the echo duration in seconds
// * `feedback`: the feedback coefficient
//----------------------------------------------------
declare echo author "Romain Michon";

echo(maxDuration,duration,feedback) = +~de.delay(maxDel,del)*feedback
with{
  maxDel = ma.SR*maxDuration;
  del = ma.SR*duration;
};

// TODO demo function for echo

//--------------------`(ef.)reverseEchoN`-------------------
// Reverse echo effect.
//
// #### Usage
//
// ```
// _ : ef.reverseEchoN(N,delay) : si.bus(N)
// ```
//
// Where:
//
// * `N`: Number of output channels desired (1 or more), a constant numerical expression
// * `delay`: echo delay (integer power of 2)
//
// #### Demo
//
// ```
// _ : dm.reverseEchoN(N) : _,_
// ```
//
// #### Description
//
// The effect uses N instances of `reverseDelayRamped` at different phases.
//
//------------------------------------------------------------
reverseEchoN(N,delMax) = _<: par(i,N,ef.reverseDelayRamped(delMax,i/N));

declare reverseEchoN author "Julius O. Smith III";
declare reverseEchoN license "STK-4.3";


//-------------------`(ef.)reverseDelayRamped`------------------
// Reverse delay with amplitude ramp.
//
// #### Usage
//
// ```
// _ : ef.reverseDelayRamped(delay,phase) : _
// ```
//
// Where:
//
// * `delay`: echo delay (integer power of 2)
// * `phase`: float between 0 and 1 giving ramp delay phase*delay
//
// #### Demo
//
// ```
// _ : ef.reverseDelayRamped(32,0.6) : _,_
// ```
//
//------------------------------------------------------------
reverseDelayRamped(delMax,phs) = rampGain * de.delay(delMax,del) with {
  rampGain = 4 * (del/delMax) * (1 - del/delMax); // suppress click when delay-line wraps around
  delOffset = int(floor(0.5 + delMax * max(0,min(0.999999,phs)))); // starting point in delay line
  startPulse = (1-1') * delOffset;
  del = int(startPulse : + ~ +(2)) & (delMax-1);
};

declare reverseDelayRamped author "Julius O. Smith III";
declare reverseDelayRamped license "STK-4.3";

//-------------------`(ef.)uniformPanToStereo`------------------
// Pan nChans channels to the stereo field, spread uniformly left to right.
//
// #### Usage
//
// ```
// si.bus(N) : ef.uniformPanToStereo(N) : _,_
// ```
//
// Where:
//
// * `N`: Number of input channels to pan down to stereo, a constant numerical expression
//
// #### Demo
//
// ```
// _,_,_ : ef.uniformPanToStereo(3) : _,_
// ```
//
//------------------------------------------------------------
uniformPanToStereo(N) = si.bus(N) <: par(i,2*N,_) :
              (par(i,N,*(i/(N-1))) :> _),
              (par(i,N,*(1-i/(N-1))) :> _);

declare uniformPanToStereo author "Julius O. Smith III";
declare uniformPanToStereo license "STK-4.3";

//---------------------`(ef.)tapeStop`-----------------------------------------
// A tape-stop effect, like putting a finger on a vinyl record player.
//
// #### Usage:
//
// ```
// _,_ : tapeStop(2, LAGRANGE_ORDER, MAX_TIME_SAMP, 
//               crossfade, gainAlpha, stopAlpha, stopTime, stop) : _,_
// ```
//
// ```
// _ : tapeStop(1, LAGRANGE_ORDER, MAX_TIME_SAMP, 
//              crossfade, gainAlpha, stopAlpha, stopTime, stop) : _
// ```
//
// Where:
//
// * `C`: The number of input and output channels.
// * `LAGRANGE_ORDER`: The order of the Lagrange interpolation on the delay line. [2-3] recommended.
// * `MAX_TIME_SAMP`: Maximum stop time in samples
// * `crossfade`: A crossfade in samples to apply when resuming normal playback. Crossfade is not applied during the enabling of the tape-stop.
// * `gainAlpha`: During the tape-stop, lower alpha stays louder longer. Safe values are in the range [.01,2].
// * `stopAlpha`: `stopAlpha==1` represents a linear deceleration (constant force). `stopAlpha<1` represents an initially weaker, then stronger force. `stopAlpha>1` represents an initially stronger, then weaker force. Safe values are in the range [.01,2].
// * `stopTime`: Desired duration of the stop time, in samples.
// * `stop`: When `stop` becomes positive, the tape-stop effect will start. When `stop` becomes zero, normal audio will resume via crossfade.
//-----------------------------------------------------------------------------
tapeStop(C, LAGRANGE_ORDER, MAX_TIME_SAMP, crossfade, gainAlpha, stopAlpha, stopTime, stop) = 
(tapeStopTick(C) ~ _) : !,si.bus(C)
with {
    tapeStopTick(C, _delaySamples) = delaySamples, circuitFinal
    with {
        // Where `stopCounter` goes from 0 to stopTime (or higher)
        // When `stopCounter` is 0, curve's output is 1.
        // When `stopCounter` is stopTime, curve's output is 0.
        curve(alpha) = 1-stopCounter/stopTime : max(0) : pow(_, alpha)
        with {
            // when stop is pulsed, count samples starting at 0
            stopCounter = *(ba.if(stop&(1-stop'),0,1))+1~_ : -(1);
        };
        minDelay = (LAGRANGE_ORDER-1)/2;
        delayFunc(curDel) = par(i, C, de.fdelayltv(LAGRANGE_ORDER, MAX_TIME_SAMP, max(curDel, minDelay)));

        delaySamples = ba.if(stop&(1-stop'), minDelay, _delaySamples) + delayDelta
        with {
            /*
            Velocity describes the velocity of the read-index in units of samples per sample.
            If the velocity is 1, then the read-index is moving as fast as the write-index
            is moving, and there is no delay. If the velocity is 0, then the read-index is "stuck"
            on a particular location. During a tape-stop, our technique is to animate velocity
            from 1 to 0 according to a curve based on stopAlpha. We discretize the accumulated
            delay with delayDelta. Note that when velocity is zero, then delayDelta is 1. At this
            moment the delay line wrote 1 new sample (as always), but our delayDelta INCREASED by one.
            This means it's playing same sample twice in a row, and the record player is motionless.
            When `stop` triggers by becoming 1, then delaySamples is reset to `minDelay`. At this moment
            we should have already been listening to the circuitNormal which is using `minDelay`.
            Therefore, there isn't a click.
            */
            velocity = curve(stopAlpha);
            delayDelta = 1-velocity;
        };

        circuitNormal = delayFunc(0); // Don't use si.bus(C) because of minDelay
        circuitSlowdown = delayFunc(delaySamples) : par(i, C, _*g)
        with {
            g = curve(gainAlpha);
        };
        circuitFinal = ba.selectmulti(actualCrossfade, (circuitNormal, circuitSlowdown), stop)
        with {
            actualCrossfade = ba.if(stop,0,crossfade); // only crossfade when resuming normal playback
        };
    };
};

declare tapeStop author "David Braun";
declare tapeStop copyright "Copyright (C) 2024 by David Braun <braun@ccrma.stanford.edu>";
declare tapeStop license "MIT-style STK-4.3 license";

//=======================================Pitch Shifting===================================
//========================================================================================

//--------------`(ef.)transpose`----------------
// A simple pitch shifter based on 2 delay lines.
// `transpose` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : transpose(w, x, s) : _
// ```
//
// Where:
//
// * `w`: the window length (samples)
// * `x`: crossfade duration duration (samples)
// * `s`: shift (semitones)
//-----------------------------------------
transpose(w, x, s, sig) = de.fdelay(maxDelay,d,sig)*ma.fmin(d/x,1) +
    de.fdelay(maxDelay,d+w,sig)*(1-ma.fmin(d/x,1))
with {
    maxDelay = 65536;
    i = 1 - pow(2, s/12);
    d = i : (+ : +(w) : fmod(_,w)) ~ _;
};

//=======================================Saturators=======================================
//========================================================================================

nonlinearityEnv = environment
{
  // Turn a function f(x) into a odd function g(x) such that
  // g(x)==f(x) for x>=0, AND g(x)==-g(-x).
  makeOdd(f, x) = x : optionalNegate : f : optionalNegate
  with {
    optionalNegate(y) = ba.if(x<0, -y, y);
  };

  // Wavefolder/Wavefolding
  // coded by David Braun
  // Reference: Chapter 10: Nonlinear Processing. Figure 10.7
  // U. Zölzer: Digital Audio Signal Processing. John Wiley & Sons Ltd, 2022.
  wavefold(width, x) = makeOdd(f, x)
  with {
    f(x) = ba.if(x>(1-2*a), tri, x) : *(g)
    with {
      // a is the adjusted width.
      a = width : aa.clip(0, 1) : *(.4);
      g = 1/(1-2*a); // peak-level normalization

      // triangle-like function
      tri = 1 - 2.5*a + a*abs(ma.frac((x-(1-2*a))/(2*a))-.5);
    };
  };
};

//---------------`(ef.)softclipQuadratic`-------------------------------------------------
// Quadratic softclip nonlinearity.
//
// #### Usage
//
// ```
// _ : softclipQuadratic : _
// ```
//
// #### References
//
// * U. Zölzer: Digital Audio Signal Processing. John Wiley & Sons Ltd, 2022.
//---------------------------------------------------------------------------------------
declare softclipQuadratic author "David Braun";
declare softclipQuadratic copyright "Copyright (C) 2024 David Braun";
declare softclipQuadratic license "MIT license";

softclipQuadratic(x) = ba.ifNcNo(2, 1, 
    absX < 1/3, 2*x,
    absX <= 2/3, ma.signum(x)*(3-(2-absX*3)^2)/3,
    ma.signum(x)
    )
with {
    absX = abs(x);
};


//---------------`(ef.)wavefold`-------------------------------------------------
// Wavefolding nonlinearity.
//
// #### Usage
//
// ```
// _ : wavefold(width) : _
// ```
//
// Where:
//
// * `width`: The width of the folded section [0..1] (float).
//---------------------------------------------------------------------------------------
declare wavefold author "David Braun";

wavefold = nonlinearityEnv.wavefold;

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

echo1s  = vgroup("echo 1000", +~(de.delay(65536, int(hslider("millisecond", 0, 0, 1000, 0.10)*ba.millisec)-1) * (hslider("feedback", 0, 0, 100, 0.1)/100.0)));
echo2s  = vgroup("echo 2000", +~(de.delay(131072, int(hslider("millisecond", 0, 0, 2000, 0.25)*ba.millisec)-1) * (hslider("feedback", 0, 0, 100, 0.1)/100.0)));
echo5s  = vgroup("echo 5000", +~(de.delay(262144, int(hslider("millisecond", 0, 0, 5000, 0.50)*ba.millisec)-1) * (hslider("feedback", 0, 0, 100, 0.1)/100.0)));
echo10s = vgroup("echo 10000", +~(de.delay(524288, int(hslider("millisecond", 0, 0, 10000, 1.00)*ba.millisec)-1) * (hslider("feedback", 0, 0, 100, 0.1)/100.0)));
echo21s = vgroup("echo 21000", +~(de.delay(1048576, int(hslider("millisecond", 0, 0, 21000, 1.00)*ba.millisec)-1) * (hslider("feedback", 0, 0, 100, 0.1)/100.0)));
echo43s = vgroup("echo 43000", +~(de.delay(2097152, int(hslider("millisecond", 0, 0, 43000, 1.00)*ba.millisec)-1) * (hslider("feedback", 0, 0, 100, 0.1)/100.0)));