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

//#################################### reverb.lib ########################################
// A library of reverb effects.
//
// It should be used using the `re` environment:
//
// ```
// re = library("reverb.lib");
// process = re.functionCall;
// ```
//
// Another option is to import `stdfaust.lib` which already contains the `re`
// environment:
//
// ```
// import("stdfaust.lib");
// process = re.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");
de = library("delay.lib");
ro = library("route.lib");
si = library("signal.lib");
fi = library("filter.lib");

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

//=============================Functions Reference========================================
//========================================================================================

//------------------------------`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 : _,_,_,_
// ```
//------------------------------------------------------------
// TODO: author JOS, revised by RM
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),_ :> _;
  };
};


//------------------------------`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 : _,_
// ```
//------------------------------------------------------------
// TODO: author JOS, revised by RM
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);
};


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


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


//--------------------------------`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>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
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)); // filter.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)
  };
};


//-------------------------------`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>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
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)))); // filter.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); // filter.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);
};


//----------------------------`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)
//------------------------------------------------------------
// TODO: author JOS, revised by RM
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,!);
};


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