reverbs.lib

//#################################### reverbs.lib ########################################
// A library of reverb effects. Its official prefix is `re`.
//########################################################################################

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

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

//########################################################################################
/************************************************************************
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!)

************************************************************************/

//=============================Schroeder Reverberators======================================
//==========================================================================================

//------------------------------`(re.)jcrev`------------------------------
// This artificial reverberator take a mono signal and output stereo
// (`satrev`) and quad (`jcrev`). They were implemented by John Chowning
// in the MUS10 computer-music language (descended from Music V by Max
// Mathews).  They are Schroeder Reverberators, well tuned for their size.
// Nowadays, the more expensive freeverb is more commonly used (see the
// Faust examples directory).
//
// `jcrev` reverb below was made from a listing of "RV", dated April 14, 1972,
// which was recovered from an old SAIL DART backup tape.
// John Chowning thinks this might be the one that became the
// well known and often copied JCREV.
//
// `jcrev` is a standard Faust function
//
// #### Usage
//
// ```
// _ : jcrev : _,_,_,_
// ```
//------------------------------------------------------------
jcrev = *(0.06) : allpass_chain <: comb_bank : mix_mtx with {
  rev1N = fi.rev1;
  rev12(len,g) = rev1N(2048,len,g);
  rev14(len,g) = rev1N(4096,len,g);
  allpass_chain =
    fi.rev2(512,347,0.7) :
    fi.rev2(128,113,0.7) :
    fi.rev2( 64, 37,0.7);
  comb_bank =
    rev12(1601,.802),
    rev12(1867,.773),
    rev14(2053,.753),
    rev14(2251,.733);
    mix_mtx = _,_,_,_ <: psum, -psum, asum, -asum : _,_,_,_ with {
    psum = _,_,_,_ :> _;
    asum = *(-1),_,*(-1),_ :> _;
  };
};


//------------------------------`(re.)satrev`------------------------------
// This artificial reverberator take a mono signal and output stereo
// (`satrev`) and quad (`jcrev`).  They were implemented by John Chowning
// in the MUS10 computer-music language (descended from Music V by Max
// Mathews).  They are Schroeder Reverberators, well tuned for their size.
// Nowadays, the more expensive freeverb is more commonly used (see the
// Faust examples directory).
//
// `satrev` was made from a listing of "SATREV", dated May 15, 1971,
// which was recovered from an old SAIL DART backup tape.
// John Chowning thinks this might be the one used on his
// often-heard brass canon sound examples, one of which can be found at
// <https://ccrma.stanford.edu/~jos/wav/FM_BrassCanon2.wav>
//
// #### Usage
//
// ```
// _ : satrev : _,_
// ```
//------------------------------------------------------------
satrev = *(0.2) <: comb_bank :> allpass_chain <: _,*(-1) with {
  rev1N = fi.rev1;
  rev11(len,g) = rev1N(1024,len,g);
  rev12(len,g) = rev1N(2048,len,g);
  comb_bank =
    rev11( 778,.827),
    rev11( 901,.805),
    rev11(1011,.783),
    rev12(1123,.764);
  rev2N = fi.rev2;
  allpass_chain =
    rev2N(128,125,0.7) :
    rev2N( 64, 42,0.7) :
    rev2N( 16, 12,0.7);
};

//======================Feedback Delay Network (FDN) Reverberators========================
//========================================================================================

//--------------------------------`(re.)fdnrev0`---------------------------------
// Pure Feedback Delay Network Reverberator (generalized for easy scaling).
// `fdnrev0` is a standard Faust function.
//
// #### Usage
//
// ```
// <1,2,4,...,N signals> <:
// fdnrev0(MAXDELAY,delays,BBSO,freqs,durs,loopgainmax,nonl) :>
// <1,2,4,...,N signals>
// ```
//
// Where:
//
// * `N`: 2, 4, 8, ...  (power of 2)
// * `MAXDELAY`: power of 2 at least as large as longest delay-line length
// * `delays`: N delay lines, N a power of 2, lengths perferably coprime
// * `BBSO`: odd positive integer = order of bandsplit desired at freqs
// * `freqs`: NB-1 crossover frequencies separating desired frequency bands
// * `durs`: NB decay times (t60) desired for the various bands
// * `loopgainmax`: scalar gain between 0 and 1 used to "squelch" the reverb
// * `nonl`: nonlinearity (0 to 0.999..., 0 being linear)
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/FDN_Reverberation.html>
//------------------------------------------------------------
fdnrev0(MAXDELAY, delays, BBSO, freqs, durs, loopgainmax, nonl)
  = (si.bus(2*N) :> si.bus(N) : delaylines(N)) ~
    (delayfilters(N,freqs,durs) : feedbackmatrix(N))
with {
  N = ba.count(delays);
  NB = ba.count(durs);
//assert(count(freqs)+1==NB);
  delayval(i) = ba.take(i+1,delays);
  dlmax(i) = MAXDELAY; // must hardwire this from argument for now
//dlmax(i) = 2^max(1,nextpow2(delayval(i))) // try when slider min/max is known
//           with { nextpow2(x) = ceil(log(x)/log(2.0)); };
// -1 is for feedback delay:
  delaylines(N) = par(i,N,(de.delay(dlmax(i),(delayval(i)-1))));
  delayfilters(N,freqs,durs) = par(i,N,filter(i,freqs,durs));
  feedbackmatrix(N) = bhadamard(N);
  vbutterfly(n) = si.bus(n) <: (si.bus(n):>bus(n/2)) , ((si.bus(n/2),(si.bus(n/2):par(i,n/2,*(-1)))) :> si.bus(n/2));
  bhadamard(2) = si.bus(2) <: +,-;
  bhadamard(n) = si.bus(n) <: (si.bus(n):>si.bus(n/2)) , ((si.bus(n/2),(si.bus(n/2):par(i,n/2,*(-1)))) :> si.bus(n/2))
                 : (bhadamard(n/2) , bhadamard(n/2));

  // Experimental nonlinearities:
  // nonlinallpass = apnl(nonl,-nonl);
  // s = nonl*PI;
  // nonlinallpass(x) = allpassnn(3,(s*x,s*x*x,s*x*x*x)); // filters.lib
     nonlinallpass = _; // disabled by default (rather expensive)

  filter(i,freqs,durs) = fi.filterbank(BBSO,freqs) : par(j,NB,*(g(j,i)))
                         :> *(loopgainmax) / sqrt(N) : nonlinallpass
  with {
    dur(j) = ba.take(j+1,durs);
    n60(j) = dur(j)*ma.SR; // decay time in samples
    g(j,i) = exp(-3.0*log(10.0)*delayval(i)/n60(j));
        // ~ 1.0 - 6.91*delayval(i)/(SR*dur(j)); // valid for large dur(j)
  };
};


//-------------------------------`(re.)zita_rev_fdn`-------------------------------
// Internal 8x8 late-reverberation FDN used in the FOSS Linux reverb zita-rev1
// by Fons Adriaensen <fons@linuxaudio.org>.  This is an FDN reverb with
// allpass comb filters in each feedback delay in addition to the
// damping filters.
//
// #### Usage
//
// ```
// bus(8) : zita_rev_fdn(f1,f2,t60dc,t60m,fsmax) : bus(8)
// ```
//
// Where:
//
// * `f1`: crossover frequency (Hz) separating dc and midrange frequencies
// * `f2`: frequency (Hz) above f1 where T60 = t60m/2 (see below)
// * `t60dc`: desired decay time (t60) at frequency 0 (sec)
// * `t60m`: desired decay time (t60) at midrange frequencies (sec)
// * `fsmax`: maximum sampling rate to be used (Hz)
//
// #### Reference
//
// * <http://www.kokkinizita.net/linuxaudio/zita-rev1-doc/quickguide.html>
// * <https://ccrma.stanford.edu/~jos/pasp/Zita_Rev1.html>
//------------------------------------------------------------
zita_rev_fdn(f1,f2,t60dc,t60m,fsmax) =
  ((si.bus(2*N) :> allpass_combs(N) : feedbackmatrix(N)) ~
   (delayfilters(N,freqs,durs) : fbdelaylines(N)))
with {
  N = 8;

  // Delay-line lengths in seconds:
  apdelays = (0.020346, 0.024421, 0.031604, 0.027333, 0.022904,
              0.029291, 0.013458, 0.019123); // feedforward delays in seconds
  tdelays = ( 0.153129, 0.210389, 0.127837, 0.256891, 0.174713,
              0.192303, 0.125000, 0.219991); // total delays in seconds
  tdelay(i) = floor(0.5 + ma.SR*ba.take(i+1,tdelays)); // samples
  apdelay(i) = floor(0.5 + ma.SR*ba.take(i+1,apdelays));
  fbdelay(i) = tdelay(i) - apdelay(i);
  // NOTE: Since SR is not bounded at compile time, we can't use it to
  // allocate delay lines; hence, the fsmax parameter:
  tdelaymaxfs(i) = floor(0.5 + fsmax*ba.take(i+1,tdelays));
  apdelaymaxfs(i) = floor(0.5 + fsmax*ba.take(i+1,apdelays));
  fbdelaymaxfs(i) = tdelaymaxfs(i) - apdelaymaxfs(i);
  nextpow2(x) = ceil(log(x)/log(2.0));
  maxapdelay(i) = int(2.0^max(1.0,nextpow2(apdelaymaxfs(i))));
  maxfbdelay(i) = int(2.0^max(1.0,nextpow2(fbdelaymaxfs(i))));

  apcoeff(i) = select2(i&1,0.6,-0.6);  // allpass comb-filter coefficient
  allpass_combs(N) =
    par(i,N,(fi.allpass_comb(maxapdelay(i),apdelay(i),apcoeff(i)))); // filters.lib
  fbdelaylines(N) = par(i,N,(de.delay(maxfbdelay(i),(fbdelay(i)))));
  freqs = (f1,f2); durs = (t60dc,t60m);
  delayfilters(N,freqs,durs) = par(i,N,filter(i,freqs,durs));
  feedbackmatrix(N) = ro.hadamard(N);

  staynormal = 10.0^(-20); // let signals decay well below LSB, but not to zero

  special_lowpass(g,f) = si.smooth(p) with {
    // unity-dc-gain lowpass needs gain g at frequency f => quadratic formula:
    p = mbo2 - sqrt(max(0,mbo2*mbo2 - 1.0)); // other solution is unstable
    mbo2 = (1.0 - gs*c)/(1.0 - gs); // NOTE: must ensure |g|<1 (t60m finite)
    gs = g*g;
    c = cos(2.0*ma.PI*f/float(ma.SR));
  };

  filter(i,freqs,durs) = lowshelf_lowpass(i)/sqrt(float(N))+staynormal
  with {
    lowshelf_lowpass(i) = gM*low_shelf1_l(g0/gM,f(1)):special_lowpass(gM,f(2));
    low_shelf1_l(G0,fx,x) = x + (G0-1)*fi.lowpass(1,fx,x); // filters.lib
    g0 = g(0,i);
    gM = g(1,i);
    f(k) = ba.take(k,freqs);
    dur(j) = ba.take(j+1,durs);
    n60(j) = dur(j)*ma.SR; // decay time in samples
    g(j,i) = exp(-3.0*log(10.0)*tdelay(i)/n60(j));
  };
};

// Stereo input delay used by zita_rev1 in both stereo and ambisonics mode:
zita_in_delay(rdel) = zita_delay_mono(rdel), zita_delay_mono(rdel) with {
  zita_delay_mono(rdel) = de.delay(8192,ma.SR*rdel*0.001) * 0.3;
};

// Stereo input mapping used by zita_rev1 in both stereo and ambisonics mode:
zita_distrib2(N) = _,_ <: fanflip(N) with {
   fanflip(4) = _,_,*(-1),*(-1);
   fanflip(N) = fanflip(N/2),fanflip(N/2);
};


//----------------------------`(re.)zita_rev1_stereo`---------------------------
// Extend `zita_rev_fdn` to include `zita_rev1` input/output mapping in stereo mode.
// `zita_rev1_stereo` is a standard Faust function.
//
// #### Usage
//
// ```
// _,_ : zita_rev1_stereo(rdel,f1,f2,t60dc,t60m,fsmax) : _,_
// ```
//
// Where:
//
// `rdel`  = delay (in ms) before reverberation begins (e.g., 0 to ~100 ms)
// (remaining args and refs as for `zita_rev_fdn` above)
//------------------------------------------------------------
zita_rev1_stereo(rdel,f1,f2,t60dc,t60m,fsmax) =
   zita_in_delay(rdel)
 : zita_distrib2(N)
 : zita_rev_fdn(f1,f2,t60dc,t60m,fsmax)
 : output2(N)
with {
 N = 8;
 output2(N) = outmix(N) : *(t1),*(t1);
 t1 = 0.37; // zita-rev1 linearly ramps from 0 to t1 over one buffer
 outmix(4) = !,ro.butterfly(2),!; // probably the result of some experimenting!
 outmix(N) = outmix(N/2),par(i,N/2,!);
};


//-----------------------------`(re.)zita_rev1_ambi`---------------------------
// Extend zita_rev_fdn to include zita_rev1 input/output mapping in
// "ambisonics mode", as provided in the Linux C++ version.
//
// #### Usage
//
// ```
// _,_ : zita_rev1_ambi(rgxyz,rdel,f1,f2,t60dc,t60m,fsmax) : _,_,_,_
// ```
//
// Where:
//
// `rgxyz` = relative gain of lanes 1,4,2 to lane 0 in output (e.g., -9 to 9)
//   (remaining args and references as for zita_rev1_stereo above)
//------------------------------------------------------------
zita_rev1_ambi(rgxyz,rdel,f1,f2,t60dc,t60m,fsmax) =
   zita_in_delay(rdel)
 : zita_distrib2(N)
 : zita_rev_fdn(f1,f2,t60dc,t60m,fsmax)
 : output4(N) // ambisonics mode
with {
  N=8;
  output4(N) = select4 : *(t0),*(t1),*(t1),*(t1);
  select4 = _,_,_,!,_,!,!,! : _,_,cross with { cross(x,y) = y,x; };
  t0 = 1.0/sqrt(2.0);
  t1 = t0 * 10.0^(0.05 * rgxyz);
};

// end jos section
/************************************************************************
************************************************************************
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.
************************************************************************
************************************************************************/

//===============================Freeverb===================================
//==========================================================================

//----------------------------`(re.)mono_freeverb`-------------------------
// A simple Schroeder reverberator primarily developed by "Jezar at Dreampoint" that
// is extensively used in the free-software world. It uses four Schroeder allpasses in
// series and eight parallel Schroeder-Moorer filtered-feedback comb-filters for each
// audio channel, and is said to be especially well tuned.
//
// `mono_freeverb` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : mono_freeverb(fb1, fb2, damp, spread) : _;
// ```
//
// Where:
//
// * `fb1`: coefficient of the lowpass comb filters (0-1)
// * `fb2`: coefficient of the allpass comb filters (0-1)
// * `damp`: damping of the lowpass comb filter (0-1)
// * `spread`: spatial spread in number of samples (for stereo)
//
// #### License
// While this version is licensed LGPL (with exception) along with other GRAME
// library functions, the file freeverb.dsp in the examples directory of older
// Faust distributions, such as faust-0.9.85, was released under the BSD license,
// which is less restrictive.
//------------------------------------------------------------
// TODO: author RM
mono_freeverb(fb1, fb2, damp, spread) = _ <: par(i,8,lbcf(combtuningL(i)+spread,fb1,damp))
	:> seq(i,4,fi.allpass_comb(1024, allpasstuningL(i)+spread, -fb2))
with{
    origSR = 44100;

    // Filters parameters
    combtuningL(0) = 1116*ma.SR/origSR : int;
    combtuningL(1) = 1188*ma.SR/origSR : int;
    combtuningL(2) = 1277*ma.SR/origSR : int;
    combtuningL(3) = 1356*ma.SR/origSR : int;
    combtuningL(4) = 1422*ma.SR/origSR : int;
    combtuningL(5) = 1491*ma.SR/origSR : int;
    combtuningL(6) = 1557*ma.SR/origSR : int;
    combtuningL(7) = 1617*ma.SR/origSR : int;

    allpasstuningL(0) = 556*ma.SR/origSR : int;
    allpasstuningL(1) = 441*ma.SR/origSR : int;
    allpasstuningL(2) = 341*ma.SR/origSR : int;
    allpasstuningL(3) = 225*ma.SR/origSR : int;
    // Lowpass Feedback Combfiler:
    // https://ccrma.stanford.edu/~jos/pasp/Lowpass_Feedback_Comb_Filter.html
    lbcf(dt, fb, damp) = (+:@(dt)) ~ (*(1-damp) : (+ ~ *(damp)) : *(fb));
};


//----------------------------`(re.)stereo_freeverb`-------------------------
// A simple Schroeder reverberator primarily developed by "Jezar at Dreampoint" that
// is extensively used in the free-software world. It uses four Schroeder allpasses in
// series and eight parallel Schroeder-Moorer filtered-feedback comb-filters for each
// audio channel, and is said to be especially well tuned.
//
// #### Usage
//
// ```
// _,_ : stereo_freeverb(fb1, fb2, damp, spread) : _,_;
// ```
//
// Where:
//
// * `fb1`: coefficient of the lowpass comb filters (0-1)
// * `fb2`: coefficient of the allpass comb filters (0-1)
// * `damp`: damping of the lowpass comb filter (0-1)
// * `spread`: spatial spread in number of samples (for stereo)
//------------------------------------------------------------
// TODO: author RM
stereo_freeverb(fb1, fb2, damp, spread) =  + <: mono_freeverb(fb1, fb2, damp,0), mono_freeverb(fb1, fb2, damp, spread);

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

// end further further contributions section