delays.lib

//#################################### delays.lib #########################################
// This library contains a collection of delay functions. Its official prefix is `de`.
//
// #### References
// * <https://github.com/grame-cncm/faustlibraries/blob/master/delays.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");
si = library("signals.lib");
fi = library("filters.lib");

declare name "Faust Delay Library";
declare version "1.2.0";

//==================================Basic Delay Functions=================================
//========================================================================================

//-------`(de.)delay`----------
// Simple `d` samples delay where `n` is the maximum delay length as a number of
// samples. Unlike the `@` delay operator, here the delay signal `d` is explicitly
// bounded to the interval [0..n]. The consequence is that delay will compile even
// if the interval of d can't be computed by the compiler.
// `delay` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : delay(n,d) : _
// ```
//
// Where:
//
// * `n`: the max delay length in samples
// * `d`: the delay length in samples (integer)
//-----------------------------
// TODO: add MBH np2
delay(n,d,x) = x @ min(n, max(0,d));


//-------`(de.)fdelay`----------
// Simple `d` samples fractional delay based on 2 interpolated delay lines where `n` is
// the maximum delay length as a number of samples.
// 
// `fdelay` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : fdelay(n,d) : _
// ```
//
// Where:
//
// * `n`: the max delay length in samples
// * `d`: the delay length in samples (float)
//-----------------------------
fdelay(n,d,x) = delay(n+1,int(d),x)*(1 - ma.frac(d)) + delay(n+1,int(d)+1,x)*ma.frac(d);


//--------------------------`(de.)sdelay`----------------------------
// s(mooth)delay: a mono delay that doesn't click and doesn't
// transpose when the delay time is changed.
//
// #### Usage
//
// ```
// _ : sdelay(n,it,d) : _
// ```
//
// Where:
//
// * `n`: the max delay length in samples
// * `it`: interpolation time (in samples), for example 1024
// * `d`: the delay length in samples (float)
//--------------------------------------------------------------------------
sdelay(n, it, d) = ctrl(it,d),_ : ddi(n)
with {
    // ddi(n,i,d0,d1)
    // DDI (Double Delay with Interpolation) : the input signal is sent to two
    // delay lines. The outputs of these delay lines are crossfaded with
    // an interpolation stage. By acting on this interpolation value one
    // can move smoothly from one delay to another. When <i> is 0 we can
    // freely change the delay time <d1> of line 1, when it is 1 we can freely change
    // the delay time <d0> of line 0.
    //
    // <n>  = maximal delay in samples
    // <i>  = interpolation value between 0 and 1 used to crossfade the outputs of the
    //        two delay lines (0.0: first delay line, 1.0: second delay line)
    // <d0> = delay time of delay line 0 in samples between 0 and <N>-1
    // <d1> = delay time of delay line 1 in samples between 0 and <N>-1
    // <  > = the input signal we want to delay
    ddi(n, i, d0, d1) = _ <: delay(n,d0), delay(n,d1) : si.interpolate(i);

    // ctrl(it,d)
    // Control logic for a Double Delay with Interpolation according to two
    //
    // USAGE : ctrl(it,d)
    // where :
    // <it> an interpolation time (in samples, for example 256)
    // <d> a delay time (in samples)
    //
    // ctrl produces 3 outputs : an interpolation value <i> and two delay
    // times <d0> and <d1>. These signals are used to control a ddi (Double Delay with Interpolation).
    // The principle is to detect changes in the input delay time d, then to
    // change the delay time of the delay line currently unused and then by a
    // smooth crossfade to remove the first delay line and activate the second one.
    //
    // The control logic has an internal state controlled by 4 elements
    // <v> : the interpolation variation (0, 1/it, -1/it)
    // <i> : the interpolation value (between 0 and 1)
    // <d0>: the delay time of line 0
    // <d1>: the delay time of line 1
    //
    // Please note that the last stage (!,_,_,_) cut <v> because it is only
    // used internally.
    ctrl(it, d) = \(v,ip,d0,d1).((nv, nip, nd0, nd1)
    with {
        // interpolation variation
        nv = ba.if (v!=0.0,                        // if variation we are interpolating
            ba.if ((ip>0.0) & (ip<1.0), v, 0),     // should we continue or not ?
            ba.if ((ip==0.0) & (d!=d0), 1.0/it,    // if true xfade from dl0 to dl1
            ba.if ((ip==1.0) & (d!=d1), -1.0/it,   // if true xfade from dl1 to dl0
            0)));                                  // nothing to change
        // interpolation value
        nip = ip+nv : min(1.0) : max(0.0);

        // update delay time of line 0 if needed
        nd0 = ba.if ((ip >= 1.0) & (d1!=d), d, d0);

        // update delay time of line 0 if needed
        nd1 = ba.if ((ip <= 0.0) & (d0!=d), d, d1);
    }) ~ (_,_,_,_) : (!,_,_,_);
};


//------------------`(de.)prime_power_delays`----------------
// Prime Power Delay Line Lengths.
//
// #### Usage
//
// ```
// si.bus(N) : prime_power_delays(N,pathmin,pathmax) : si.bus(N);
// ```
//
// Where:
//
// * `N`: positive integer up to 16 (for higher powers of 2, extend 'primes' array below)
// * `pathmin`: minimum acoustic ray length in the reverberator (in meters)
// * `pathmax`: maximum acoustic ray length (meters) - think "room size"
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Prime_Power_Delay_Line.html>
//------------------------------------------------------------
declare prime_power_delays author "Julius O. Smith III";

prime_power_delays(N,pathmin,pathmax) = par(i,N,delayvals(i)) with {
  Np = 16;
  primes = 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53;
  prime(n) = primes : ba.selector(n,Np); // math.lib

  // Prime Power Bounds [matlab: floor(log(maxdel)./log(primes(53)))]
  maxdel = 8192; // more than 63 meters at 44100 samples/sec & 343 m/s
  ppbs = 13,8,5,4, 3,3,3,3, 2,2,2,2, 2,2,2,2; // 8192 is enough for all
  ppb(i) = ba.take(i+1,ppbs);

  // Approximate desired delay-line lengths using powers of distinct primes:
  c = 343; // soundspeed in m/s at 20 degrees C for dry air
  dmin = ma.SR*pathmin/c;
  dmax = ma.SR*pathmax/c;
  dl(i) = dmin * (dmax/dmin)^(i/float(N-1)); // desired delay in samples
  ppwr(i) = floor(0.5+log(dl(i))/log(prime(i))); // best prime power
  delayvals(i) = prime(i)^ppwr(i); // each delay a power of a distinct prime
};


//===============================Lagrange Interpolation===================================
//========================================================================================

//----------------------`(de.)fdelaylti` and `(de.)fdelayltv`-------------------------
// Fractional delay line using Lagrange interpolation.
//
// #### Usage
//
// ```
// _ : fdelaylt[i|v](N, n, d) : _
// ```
//
// Where:
//
// * `N=1,2,3,...` is the order of the Lagrange interpolation polynomial (constant numerical expression)
// * `n`: the max delay length in samples
// * `d`: the delay length in samples
//
// `fdelaylti` is most efficient, but designed for constant/slowly-varying delay.
// `fdelayltv` is more expensive and more robust when the delay varies rapidly.
//
// Note: the requested delay should not be less than `(N-1)/2`.
//
// #### References
//
// * <https://ccrma.stanford.edu/~jos/pasp/Lagrange_Interpolation.html>
//   - [fixed-delay case](https://ccrma.stanford.edu/~jos/Interpolation/Efficient_Time_Invariant_Lagrange_Interpolation.html)
//   - [variable-delay case](https://ccrma.stanford.edu/~jos/Interpolation/Time_Varying_Lagrange_Interpolation.html)
// * Timo I. Laakso et al., "Splitting the Unit Delay - Tools for Fractional
//         Delay Filter Design", IEEE Signal Processing Magazine,
//         vol. 13, no. 1, pp. 30-60, Jan 1996.
// * Philippe Depalle and Stephan Tassart, "Fractional Delay Lines using
//         Lagrange Interpolators", ICMC Proceedings, pp. 341-343, 1996.
//------------------------------------------------------------
declare fdelaylti author "Julius O. Smith III";

fdelaylti(N,n,d,x) = delay(n,id,x) <: seq(i,N,section(i)) : !,_
with {
  o = (N-1.00001)/2; // offset to ~center FIR interpolator
  dmo = d - o; // assumed nonnegative [d > (N-1)/2]
  id = int(dmo);
  fd = o + ma.frac(dmo);
  section(i,x,y) = (x-x') * c(i) <: _,+(y);
  c(i) = (i - fd)/(i+1);
};

declare fdelayltv author "Julius O. Smith III";

fdelayltv(N,n,d,x) = sum(i, N+1, delay(n,id+i,x) * h(N,fd,i))
with {
  o = (N-1.00001)/2;  // ~center FIR interpolator
  dmo = d - o; // assumed nonnegative [d > (N-1)/2]
  id = int(dmo);
  fd = o + ma.frac(dmo);
  h(N,d,n) = facs1(N,d,n) * facs2(N,d,n);
  facs1(N,d,n) = select2(n,1,prod(k,max(1,n),select2(k<n,1,fac(d,n,k))));
  facs2(N,d,n) = select2(n<N,1,prod(l,max(1,N-n),fac(d,n,l+n+1)));
  fac(d,n,k) = (d-k)/((n-k)+(n==k));
  // Explicit formula for Lagrange interpolation coefficients:
  // h_d(n) = \prod_{\stackrel{k=0}{k\neq n}}^N \frac{d-k}{n-k}, n=0:N
};

//------------------`(de.)fdelay[N]`----------------------------
// For convenience, `fdelay1`, `fdelay2`, `fdelay3`, `fdelay4`, `fdelay5`
// are also available where `N` is the order of the interpolation, built using `fdelayltv`.
//------------------------------------------------------------
declare fdelay1 author "Julius O. Smith III";
declare fdelay2 author "Julius O. Smith III";
declare fdelay3 author "Julius O. Smith III";
declare fdelay4 author "Julius O. Smith III";
declare fdelay5 author "Julius O. Smith III";

fdelay1 = fdelayltv(1);
fdelay2 = fdelayltv(2);
fdelay3 = fdelayltv(3);
fdelay4 = fdelayltv(4);
fdelay5 = fdelayltv(5);

//==============================Thiran Allpass Interpolation==============================
// Thiran Allpass Interpolation.
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Thiran_Allpass_Interpolators.html>
//========================================================================================

//----------------`(de.)fdelay[N]a`-------------
// Delay lines interpolated using Thiran allpass interpolation.
//
// #### Usage
//
// ```
// _ : fdelay[N]a(n, d) : _
// ```
//
// (exactly like `fdelay`)
//
// Where:
//
// * `N=1,2,3, or 4` is the order of the Thiran interpolation filter (constant numerical expression),
//    and the delay argument is at least `N-1/2`. First-order: `d` at least 0.5, second-order: `d` at least 1.5,
//    third-order: `d` at least 2.5, fourth-order: `d` at least 3.5.
// * `n`: the max delay length in samples
// * `d`: the delay length in samples
//
// #### Note
//
// The interpolated delay should not be less than `N-1/2`.
// (The allpass delay ranges from `N-1/2` to `N+1/2`).
// This constraint can be alleviated by altering the code,
// but be aware that allpass filters approach zero delay
// by means of pole-zero cancellations.
//
// Delay arguments too small will produce an UNSTABLE allpass!
//
// Because allpass interpolation is recursive, it is not as robust
// as Lagrange interpolation under time-varying conditions
// (you may hear clicks when changing the delay rapidly).
//
//------------------------------------------------------------
declare fdelay1a author "Julius O. Smith III";

fdelay1a(n,d,x) = delay(n,id,x) : fi.tf1(eta,1,eta)
with {
  o = 0.49999; // offset to make life easy for allpass
  dmo = d - o; // assumed nonnegative
  id = int(dmo);
  fd = o + ma.frac(dmo);
  eta = (1-fd)/(1+fd); // allpass coefficient
};

declare fdelay2a author "Julius O. Smith III";
fdelay2a(n,d,x) = delay(n,id,x) : fi.tf2(a2,a1,1,a1,a2)
with {
  o = 1.49999;
  dmo = d - o; // delay range is [order-1/2, order+1/2]
  id = int(dmo);
  fd = o + ma.frac(dmo);
  a1o2 = (2-fd)/(1+fd); // share some terms (the compiler does this anyway)
  a1 = 2*a1o2;
  a2 = a1o2*(1-fd)/(2+fd);
};

declare fdelay3a author "Julius O. Smith III";
fdelay3a(n,d,x) = delay(n,id,x) : fi.iir((a3,a2,a1,1),(a1,a2,a3))
with {
  o = 2.49999;
  dmo = d - o;
  id = int(dmo);
  fd = o + ma.frac(dmo);
  a1o3 = (3-fd)/(1+fd);
  a2o3 = a1o3*(2-fd)/(2+fd);
  a1 = 3*a1o3;
  a2 = 3*a2o3;
  a3 = a2o3*(1-fd)/(3+fd);
};

declare fdelay4a author "Julius O. Smith III";
fdelay4a(n,d,x) = delay(n,id,x) : fi.iir((a4,a3,a2,a1,1),(a1,a2,a3,a4))
with {
  o = 3.49999;
  dmo = d - o;
  id = int(dmo);
  fd = o + ma.frac(dmo);
  a1o4 = (4-fd)/(1+fd);
  a2o6 = a1o4*(3-fd)/(2+fd);
  a3o4 = a2o6*(2-fd)/(3+fd);
  a1 = 4*a1o4;
  a2 = 6*a2o6;
  a3 = 4*a3o4;
  a4 = a3o4*(1-fd)/(4+fd);
};

//==================================Others================================================
//========================================================================================

//----------------`(de.)multiTapSincDelay`-------------
//
// Variable delay line using multi-tap sinc interpolation.
//
// This function implements a continuously variable delay line by superposing (2K+2) auxiliary delayed signals
// whose positions and gains are determined by a sinc-based interpolation method. It extends the traditional
// crossfade delay technique to significantly reduce spectral coloration artifacts, which are problematic in
// applications like Wave Field Synthesis (WFS) and auralization.
//
// Operation:
//
//   - If tau1 and tau2 are very close (|tau2 - tau1| ≈ 0), a simple fixed fractional delay is applied
//   - Otherwise, a variable delay is synthesized by:
//
//      - Computing (2K+2) taps symmetrically distributed around tau1 and tau2
//      - Applying sinc-based weighting to each tap, based on its offset from the target interpolated delay tau
//      - Summing all the weighted taps to produce the output
//
// Features:
//
//   - Smooth delay variation without introducing Doppler pitch shifts
//   - Significant reduction of comb-filter coloration compared to classical crossfading
//   - Switching between fixed and variable delay modes to ensure stability
//
// #### Usage
// ```
//  _ : multiTapSincDelay(K, MaxDelay, tau1, tau2, alpha) : _
// 
// ```
//
// Where:
//
// * `K (integer)`: number of auxiliary tap pairs (a constant numerical expression). Total number of taps = 2*K + 2
// * `MaxDelay`: maximum allowable delay in samples (buffer size)
// * `tau1`: initial delay in samples (can be fractional)
// * `tau2`: target delay in samples (can be fractional)
// * `alpha`: interpolation factor between tau1 and tau2 (in [0,1] with 0 = tau1, 1 = tau2)
//
// #### Reference
//   - T. Carpentier, "Implementation of a continuously variable delay line by crossfading between several tap delays", 2024: <https://hal.science/hal-04646939>
//------------------------------------------------------------

multiTapSincDelay(K, MaxDelay, tau1, tau2, alpha, input) = output with {

    // Total number of taps used in the sum (2K + 2).
    numTaps = 2 * K + 2;

    // Normalized Sinc function: sinc(x) = sin(pi*x)/(pi*x)
    // Handles the case where x ≈ 0 to avoid division by zero and return 1.
    sinc(x) = ba.if(abs(x) < ma.EPSILON, 1, sin(ma.PI * x) / (ma.PI * x));

    // Computes the difference between two delays: delta = tau2 - tau1
    delta(t1, t2) = t2 - t1;

    // Computes the target interpolated delay: tau = (1-alpha)*tau1 + alpha*tau2
    // This is the center of the sinc function for the hk gains.
    tau(t1, t2, a) = (1 - a) * t1 + a * t2;

    // Computes the position (delay in samples) of the i-th tap (tk),
    // according to equation (17) of the article.
    tk(i, t1, t2) = ba.if(i <= K,
                        t1 - (K - i) * delta(t1, t2),
                        t2 + (i - K - 1) * delta(t1, t2));

    // Computes the gain (hk) of the i-th tap, according to equation (19) of the paper.
    // hk = sinc((tk - tau) / delta)
    hk(i, t1, t2, a) = sinc(offset) with {
        d      = delta(t1, t2);
        denom  = ba.if(abs(d) < ma.EPSILON, 1, d);
        offset = (tk(i, t1, t2) - tau(t1, t2, a)) / denom;
    };

    // Computes the sum of delayed signals weighted by their gains.
    // This is the core of the variable multi-tap effect.
    tapSum(sig, t1, t2, a) =
        // Iterates from i = 0 to numTaps-1
        sum(i, numTaps,
            // For each tap 'i':
            // 1. Delays the signal 'sig' by tk(i) samples (fractional)
            //    using an internal delay line (de.fdelay).
            // 2. Multiplies the delayed signal by the gain hk(i).
            sig : de.fdelay(MaxDelay, tk(i, t1, t2)) * hk(i, t1, t2, a)
        );

    // Simple function to apply a fixed delay (potentially fractional).
    // Used when tau1 and tau2 are very close.
    fixedDelay(x, t) = x : de.fdelay(MaxDelay, t);

    // Determines if delays tau1 and tau2 are close enough to be
    // considered fixed (to avoid unstable calculations when delta ≈ 0).
    // isFixed equals 1 if delta ≈ 0, otherwise 0.
    isFixed = abs(delta(tau1, tau2)) < ma.EPSILON;

    // Computes the output corresponding to the fixed delay (used if isFixed=1).
    fd = fixedDelay(input, tau1);            

    // Computes the output corresponding to the variable multi-tap delay (used if isFixed=0).
    mt = tapSum(input, tau1, tau2, alpha);    

    // Selects the final output based on isFixed.
    output = ba.if(isFixed, fd, mt); 
};


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

delay1s(d) = delay(65536,d);
delay2s(d) = delay(131072,d);
delay5s(d) = delay(262144,d);
delay10s(d) = delay(524288,d);
delay21s(d) = delay(1048576,d);
delay43s(d) = delay(2097152,d);

fdelay1s(d) = fdelay(65536,d);
fdelay2s(d) = fdelay(131072,d);
fdelay5s(d) = fdelay(262144,d);
fdelay10s(d) = fdelay(524288,d);
fdelay21s(d) = fdelay(1048576,d);
fdelay43s(d) = fdelay(2097152,d);