faust-common 0.9.95~repack1-2 / usr / share / faust / delay.lib

//#################################### delay.lib #########################################
// This library contains a collection of delay functions.
//
// It should be used using the `de` environment:
//
// ```
// de = library("delay.lib");
// process = de.functionCall;
// ```
//
// Another option is to import `stdfaust.lib` which already contains the `de`
// environment:
//
// ```
// import("stdfaust.lib");
// process = de.functionCall;
// ```
//########################################################################################

/************************************************************************
************************************************************************
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("math.lib");
ba = library("basic.lib");
si = library("signal.lib");
fi = library("filter.lib");


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

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

//-------`delay`----------
// Simple `d` samples delay where `n` is the maximum delay length as a number of 
// samples (it needs to be a power of 2). Unlike the `@` delay operator, this
// function allows to preallocate memory which means that `d` can be changed dynamically
// at run time as long as it remains smaller than `n`.
// `delay` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : delay(n,d) : _
// ``` 
//
// Where:
//
// * `n`: the max delay length as a power of 2
// * `d`: the delay length as a number of samples (integer)
//-----------------------------
// TODO: add MBH np2
delay(n,d,x) = x@(int(d)&(n-1));


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


//--------------------------`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,dt) : _
// ```
//
// Where :
//
// * `N`: maximal delay in samples (must be a constant power of 2, for example 65536)
// * `it`: interpolation time (in samples) for example 1024
// * `dt`: delay time (in samples)
//--------------------------------------------------------------------------
sdelay(N, it, dt) = ctrl(it,dt),_ : 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 (must be a power of 2, for example 65536)
	// <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,dt)
	// Control logic for a Double Delay with Interpolation according to two 
	//
	// USAGE : ctrl(it,dt)
	// where : 
	// <it> an interpolation time (in samples, for example 256)
	// <dt> 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 dt, 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, dt) = \(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) & (dt!=d0),  1.0/it,		// if true xfade from dl0 to dl1
			ba.if ((ip==1.0) & (dt!=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!=dt), dt, d0);

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


// ----------`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>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
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===================================
//========================================================================================

//----------------------`fdelaylti` and `fdelayltv`-------------------------
// Fractional delay line using Lagrange interpolation.
//
// #### Usage
// 
// ```
// _ : fdelaylt[i|v](order, maxdelay, delay, inputsignal) : _
// ```
//
// Where `order=1,2,3,...` is the order of the Lagrange interpolation polynomial.
//
// `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>
// * 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.
//------------------------------------------------------------
// TODO: author JOS, revised by RM
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);
};

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 >=0 [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
};

//------------------`fdelay[n]`----------------------------
// For convenience, `fdelay1`, `fdelay2`, `fdelay3`, `fdelay4`, `fdelay5`
// are also available where n is the order of the interpolation.
//------------------------------------------------------------
// TODO: author JOS, revised by RM
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>
//========================================================================================

//----------------`fdelay[n]a`-------------
// Delay lines interpolated using Thiran allpass interpolation.
//
// #### Usage 
//
// ```
// _ : fdelay[N]a(maxdelay, delay, inputsignal) : _
// ```
//
// (exactly like `fdelay`)
//
// Where: 
//
// * `N`=1,2,3, or 4 is the order of the Thiran interpolation filter,
// and the delay argument is at least N - 1/2.
//
// #### 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.
// The delay range `[N-1/2`,`N+1/2]` is not optimal. What is?
//
// 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.)
//
// First-order allpass interpolation, delay d in [0.5,1.5]
//------------------------------------------------------------
// TODO: author JOS, revised by RM
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
};

// second-order allpass delay in [1.5,2.5]
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);
};

// third-order allpass delay in [2.5,3.5]
// delay d should be at least 2.5
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);
};

// fourth-order allpass delay in [3.5,4.5]
// delay d should be at least 3.5
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);
};

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