//################################ filter.lib ########################################
// A library of filters and of more advanced filter-based sound processor organized 
// in 18 sections:
//
// * Basic Filters
// * Comb Filters
// * Direct-Form Digital Filter Sections
// * Direct-Form Second-Order Biquad Sections
// * Ladder/Lattice Digital Filters
// * Useful Special Cases
// * Ladder/Lattice Allpass Filters
// * Digital Filter Sections Specified as Analog Filter Sections
// * Simple Resonator Filters
// * Butterworth Lowpass/Highpass Filters
// * Special Filter-Bank Delay-Equalizing Allpass Filters
// * Elliptic (Cauer) Lowpass Filters
// * Elliptic Highpass Filters
// * Butterworth Bandpass/Bandstop Filters
// * Elliptic Bandpass Filters
// * Parametric Equalizers (Shelf, Peaking)
// * Mth-Octave Filter-Banks
// * Arbritary-Crossover Filter-Banks and Spectrum Analyzers
//
// It should be used using the `fi` environment:
//
// ```
// fi = library("filter.lib");
// process = fi.functionCall;
// ```
//
// Another option is to import `stdfaust.lib` which already contains the `fi`
// environment:
//
// ```
// import("stdfaust.lib");
// process = fi.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");
an = library("analyzer.lib");

declare name "Faust Filter Library";
declare version "2.0";

//===============================Basic Filters============================================
//========================================================================================

//----------------------`zero`--------------------------
// One zero filter. Difference equation: y(n) = x(n) - z * x(n-1).
//
// #### Usage
//
// ```
// _ : zero(z) : _
// ```
//
// Where:
//
// * `z`: location of zero along real axis in z-plane
//
// #### Reference 
// <https://ccrma.stanford.edu/~jos/filters/One_Zero.html>
//----------------------------------------------------------
// TODO: author JOS, revised by RM 
zero(z) =  _ <: _,mem : _,*(z) : -;

//------------------------`pole`---------------------------
// One pole filter. Could also be called a "leaky integrator". 
// Difference equation: y(n) = x(n) + p * y(n-1).
//
// #### Usage
//
// ```
// _ : pole(z) : _
// ```
//
// Where:
//
// * `p`: pole location = feedback coefficient
//
// #### Reference 
// <https://ccrma.stanford.edu/~jos/filters/One_Pole.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
pole(p) = + ~ *(p);

//----------------------`integrator`--------------------------
// Same as `pole(1)` [implemented separately for block-diagram clarity].
//------------------------------------------------------------
// TODO: author JOS, revised by RM
integrator = + ~ _ ;

//-------------------`dcblockerat`-----------------------
// DC blocker with configurable break frequency. 
// The amplitude response is substantially flat above fb,
// and sloped at about +6 dB/octave below fb.
// Derived from the analog transfer function 
// H(s) = s / (s + 2*PI*fb) 
// by the low-frequency-matching bilinear transform method
// (i.e., the standard frequency-scaling constant 2*SR).
//
// #### Usage
//
// ```
// _ : dcblockerat(fb) : _
// ```
//
// Where:
//
// * `fb`: "break frequency" in Hz, i.e., -3 dB gain frequency.
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM 
dcblockerat(fb) = *(b0) : zero(1) : pole(p)
with {
  wn = ma.PI*fb/ma.SR;
  b0 = 1.0 / (1 + wn);
  p = (1 - wn) * b0;
};

//----------------------`dcblocker`--------------------------
// DC blocker. Default dc blocker has -3dB point near 35 Hz (at 44.1 kHz)
// and high-frequency gain near 1.0025 (due to no scaling).
// `dcblocker` is as standard Faust function.
//
// #### Usage
//
// ```
// _ : dcblocker : _
// ```
//------------------------------------------------------------
// TODO: author JOS, revised by RM 
dcblocker = zero(1) : pole(0.995);

//=======================================Comb Filters=====================================
//========================================================================================

//------`ff_comb`--------
// Feed-Forward Comb Filter. Note that `ff_comb` requires integer delays  
// (uses `delay`  internally).
// `ff_comb` is a standard Faust function.
// 
// #### Usage
//
// ```
// _ : ff_comb(maxdel,intdel,b0,bM) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `intdel`: current (integer) comb-filter delay between 0 and maxdel
// * `del`: current (float) comb-filter delay between 0 and maxdel
// * `b0`: gain applied to delay-line input
// * `bM`: gain applied to delay-line output and then summed with input
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Feedforward_Comb_Filters.html>
//------------------------------------------------------------
// TODO: author JOS
ff_comb (maxdel,M,b0,bM) = _ <: *(b0), bM *  de.delay(maxdel,M) : + ;


//------`ff_fcomb`--------
// Feed-Forward Comb Filter. Note that `ff_fcomb` takes floating-point delays 
// (uses `fdelay` internally).
// `ff_fcomb` is a standard Faust function.
// 
// #### Usage
//
// ```
// _ : ff_fcomb(maxdel,del,b0,bM) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `intdel`: current (integer) comb-filter delay between 0 and maxdel
// * `del`: current (float) comb-filter delay between 0 and maxdel
// * `b0`: gain applied to delay-line input
// * `bM`: gain applied to delay-line output and then summed with input
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Feedforward_Comb_Filters.html>
//------------------------------------------------------------
// TODO: author JOS
ff_fcomb(maxdel,M,b0,bM) = _ <: *(b0), bM * de.fdelay(maxdel,M) : + ;

//-----------`ffcombfilter`-------------------
// Typical special case of `ff_comb()` where: `b0 = 1`.
//------------------------------------------------------------
// TODO: author JOS, revised by RM
ffcombfilter(maxdel,del,g) = ff_comb(maxdel,del,1,g);


//-----------------------`fb_comb`-----------------------
// Feed-Back Comb Filter (integer delay).
//
// #### Usage 
//
// ```
// _ : fb_comb(maxdel,intdel,b0,aN) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `intdel`: current (integer) comb-filter delay between 0 and maxdel
// * `del`: current (float) comb-filter delay between 0 and maxdel
// * `b0`: gain applied to delay-line input and forwarded to output
// * `aN`: minus the gain applied to delay-line output before summing with the input 
// 	and feeding to the delay line
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Feedback_Comb_Filters.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
fb_comb (maxdel,N,b0,aN) = (+ <:  de.delay(maxdel,N-1),_) ~ *(-aN) : !,*(b0):mem;


//-----------------------`fb_fcomb`-----------------------
// Feed-Back Comb Filter (floating point delay).
//
// #### Usage 
//
// ```
// _ : fb_fcomb(maxdel,del,b0,aN) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `intdel`: current (integer) comb-filter delay between 0 and maxdel
// * `del`: current (float) comb-filter delay between 0 and maxdel
// * `b0`: gain applied to delay-line input and forwarded to output
// * `aN`: minus the gain applied to delay-line output before summing with the input 
// 	and feeding to the delay line
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Feedback_Comb_Filters.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
fb_fcomb(maxdel,N,b0,aN) = (+ <: de.fdelay(maxdel,float(N)-1.0),_) ~ *(-aN) : !,*(b0):mem;

//-----------------------`rev1`-----------------------
// Special case of `fb_comb` (`rev1(maxdel,N,g)`).
// The "rev1 section" dates back to the 1960s in computer-music reverberation.
// See the `jcrev` and `brassrev` in `reverb.lib` for usage examples.
//------------------------------------------------------------
// TODO: author JOS
rev1(maxdel,N,g) = fb_comb (maxdel,N,1,-g);

//-----`fbcombfilter` and `ffbcombfilter`------------
// Other special cases of Feed-Back Comb Filter.
//
// #### Usage 
//
// ```
// _ : fbcombfilter(maxdel,intdel,g) : _
// _ : ffbcombfilter(maxdel,del,g) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `intdel`: current (integer) comb-filter delay between 0 and maxdel
// * `del`: current (float) comb-filter delay between 0 and maxdel
// * `g`: feedback gain
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Feedback_Comb_Filters.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
fbcombfilter(maxdel,intdel,g) = (+ : de.delay(maxdel,intdel)) ~ *(g);
ffbcombfilter(maxdel,del,g) = (+ : de.fdelay(maxdel,del)) ~ *(g);


//-------------------`allpass_comb`-----------------
// Schroeder Allpass Comb Filter. Note that 
//
// ```
// allpass_comb(maxlen,len,aN) = ff_comb(maxlen,len,aN,1) : fb_comb(maxlen,len-1,1,aN);
// ```
//
// which is a direct-form-1 implementation, requiring two delay lines.
// The implementation here is direct-form-2 requiring only one delay line.
//
// #### Usage
//
// ``` 
// _ : allpass_comb (maxdel,intdel,aN) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `intdel`: current (integer) comb-filter delay between 0 and maxdel
// * `del`: current (float) comb-filter delay between 0 and maxdel
// * `aN`: minus the feedback gain
// 
// #### References
//
// * <https://ccrma.stanford.edu/~jos/pasp/Allpass_Two_Combs.html>
// * <https://ccrma.stanford.edu/~jos/pasp/Schroeder_Allpass_Sections.html>
// * <https://ccrma.stanford.edu/~jos/filters/Four_Direct_Forms.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
allpass_comb(maxdel,N,aN) = (+ <: de.delay(maxdel,N-1),*(aN)) ~ *(-aN) : mem,_ : + ;


//-------------------`allpass_fcomb`-----------------
// Schroeder Allpass Comb Filter. Note that 
//
// ```
// allpass_comb(maxlen,len,aN) = ff_comb(maxlen,len,aN,1) : fb_comb(maxlen,len-1,1,aN);
// ```
//
// which is a direct-form-1 implementation, requiring two delay lines.
// The implementation here is direct-form-2 requiring only one delay line.
//
// `allpass_fcomb` is a standard Faust library.
//
// #### Usage
//
// ``` 
// _ : allpass_comb (maxdel,intdel,aN) : _
// _ : allpass_fcomb(maxdel,del,aN) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `intdel`: current (float) comb-filter delay between 0 and maxdel
// * `del`: current (float) comb-filter delay between 0 and maxdel
// * `aN`: minus the feedback gain
// 
// #### References
//
// * <https://ccrma.stanford.edu/~jos/pasp/Allpass_Two_Combs.html>
// * <https://ccrma.stanford.edu/~jos/pasp/Schroeder_Allpass_Sections.html>
// * <https://ccrma.stanford.edu/~jos/filters/Four_Direct_Forms.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
allpass_fcomb(maxdel,N,aN) = (+ <: de.fdelay(maxdel,N-1),*(aN)) ~ *(-aN) : mem,_ : + ;


//-----------------------`rev2`-----------------------
// Special case of `allpass_comb` (`rev2(maxlen,len,g)`).
// The "rev2 section" dates back to the 1960s in computer-music reverberation.
// See the `jcrev` and `brassrev` in `reverb.lib` for usage examples.
//------------------------------------------------------------
rev2(maxlen,len,g) = allpass_comb(maxlen,len,-g);

//-------------------`allpass_fcomb5` and `allpass_fcomb1a`-----------------
// Same as `allpass_fcomb` but use `fdelay5` and `fdelay1a` internally 
// (Interpolation helps - look at an fft of faust2octave on
//
// ```
// `1-1' <: allpass_fcomb(1024,10.5,0.95), allpass_fcomb5(1024,10.5,0.95);`).
// ```
//------------------------------------------------------------
// TODO: author JOS, revised by RM
allpass_fcomb5(maxdel,N,aN) = (+ <: de.fdelay5(maxdel,N-1),*(aN)) ~ *(-aN) : mem,_ : + ;
allpass_fcomb1a(maxdel,N,aN) = (+ <: de.fdelay1a(maxdel,N-1),*(aN)) ~ *(-aN) : mem,_ : + ;

//========================Direct-Form Digital Filter Sections=============================
//========================================================================================

// Specified by transfer-function polynomials B(z)/A(z) as in matlab

//----------------------------`iir`-------------------------------
// Nth-order Infinite-Impulse-Response (IIR) digital filter,
// implemented in terms of the Transfer-Function (TF) coefficients.
// Such filter structures are termed "direct form".
//
// `iir` is a standard Faust function.
//
// #### Usage
//
// ```
//   _ : iir(bcoeffs,acoeffs) : _
// ```
//
// Where:
//
// * `order`: filter order (int) = max(#poles,#zeros)
// * `bcoeffs`: (b0,b1,...,b_order) = TF numerator coefficients
// * `acoeffs`: (a1,...,a_order) = TF denominator coeffs (a0=1)
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/filters/Four_Direct_Forms.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
iir(bv,av) = ma.sub ~ fir(av) : fir(bv);

//----------------------------- `fir` ---------------------------------
// FIR filter (convolution of FIR filter coefficients with a signal)
//
// #### Usage
//
// ``` 
// _ : fir(bv) : _
// ```
//
// `fir` is standard Faust function.
//
// Where: 
//
// * `bv` = b0,b1,...,bn is a parallel bank of coefficient signals.
//
// #### Note 
//
// `bv` is processed using pattern-matching at compile time,
//       so it must have this normal form (parallel signals).
//
// #### Example
//
// Smoothing white noise with a five-point moving average:
//
// ```
// bv = .2,.2,.2,.2,.2;
// process = noise : fir(bv);
// ```
//
// Equivalent (note double parens):
//
// ```
// process = noise : fir((.2,.2,.2,.2,.2));
// ```
//------------------------------------------------------------
// TODO: author JOS, revised by RM
//fir(bv) = conv(bv);
fir((b0,bv)) = _ <: *(b0), R(1,bv) :> _ with {
	R(n,(bn,bv)) = (@(n):*(bn)), R(n+1,bv);
	R(n, bn)     = (@(n):*(bn));              };
fir(b0) = *(b0);

//---------------`conv` and `convN`-------------------------------
// Convolution of input signal with given coefficients.
//
// #### Usage
//
// ```
// _ : conv((k1,k2,k3,...,kN)) : _; // Argument = one signal bank
// _ : convN(N,(k1,k2,k3,...)) : _; // Useful when N < count((k1,...))
// ```
//------------------------------------------------------------
// TODO: author JOS, revised by RM
//convN(N,kv,x) = sum(i,N,take(i+1,kv) * x@i); // take() defined in math.lib
convN(N,kv)     = sum(i,N, @(i)*take(i+1,kv)); // take() defined in math.lib
//conv(kv,x) = sum(i,count(kv),take(i+1,kv) * x@i); // count() from math.lib
conv(kv) = fir(kv);

//----------------`tf1`, `tf2` and `tf3`----------------------
// tfN = N'th-order direct-form digital filter.
//
// #### Usage
//
// ```
// _ : tf1(b0,b1,a1) : _
// _ : tf2(b0,b1,b2,a1,a2) : _
// _ : tf3(b0,b1,b2,b3,a1,a2,a3) : _
// ```
//
// Where:
// 
// * `a`: the poles
// * `b`: the zeros
//
// #### Reference
// <https://ccrma.stanford.edu/~jos/fp/Direct_Form_I.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
tf1(b0,b1,a1) = _ <: *(b0), (mem : *(b1)) :> + ~ *(0-a1);
tf2(b0,b1,b2,a1,a2) = iir((b0,b1,b2),(a1,a2)); // cf. TF2 in music.lib)
// tf2 is a variant of tf22 below with duplicated mems
tf3(b0,b1,b2,b3,a1,a2,a3) = iir((b0,b1,b2,b3),(a1,a2,a3));

// "Original" version for music.lib. This is here for comparison but people should
// use tf2 instead
TF2(b0,b1,b2,a1,a2) = sub ~ conv2(a1,a2) : conv3(b0,b1,b2)
with {
	conv3(k0,k1,k2,x) 	= k0*x + k1*x' + k2*x'';
	conv2(k0,k1,x) 		= k0*x + k1*x';
	sub(x,y)			= y-x;
};

//------------`notchw`--------------
// Simple notch filter based on a biquad (`tf2`).
// `notchw` is a standard Faust function. 
//
// #### Usage:
//
// ```
// _ : notchw(width,freq) : _
// ```
//
// Where:
//
// * `width`: "notch width" in Hz (approximate)
// * `freq`: "notch frequency" in Hz
//
// #### Reference
// <https://ccrma.stanford.edu/~jos/pasp/Phasing_2nd_Order_Allpass_Filters.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM 
notchw(width,freq) = tf2(b0,b1,b2,a1,a2)
with {
  fb = 0.5*width; // First design a dcblockerat(width/2)
  wn = PI*fb/SR;
  b0db = 1.0 / (1 + wn);
  p = (1 - wn) * b0db; // This is our pole radius.
  // Now place unit-circle zeros at desired angles:
  tn = 2*PI*freq/SR;
  a2 = p * p;
  a2p1 = 1+a2;
  a1 = -a2p1*cos(tn);
  b1 = a1;
  b0 = 0.5*a2p1;
  b2 = b0;
};

//======================Direct-Form Second-Order Biquad Sections==========================
// Direct-Form Second-Order Biquad Sections
// 
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/filters/Four_Direct_Forms.html>
//========================================================================================

//----------------`tf21`, `tf22`, `tf22t` and `tf21t`----------------------
// tfN = N'th-order direct-form digital filter where: 
//
// * `tf21` is tf2, direct-form 1
// * `tf22` is tf2, direct-form 2
// * `tf22t` is tf2, direct-form 2 transposed
// * `tf21t` is tf2, direct-form 1 transposed
//
// #### Usage
//
// ```
// _ : tf21(b0,b1,b2,a1,a2) : _
// _ : tf22(b0,b1,b2,a1,a2) : _
// _ : tf22t(b0,b1,b2,a1,a2) : _
// _ : tf21t(b0,b1,b2,a1,a2) : _
// ```
//
// Where:
// 
// * `a`: the poles
// * `b`: the zeros
//
// #### Reference
// <https://ccrma.stanford.edu/~jos/fp/Direct_Form_I.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
tf21(b0,b1,b2,a1,a2) = // tf2, direct-form 1:
    _ <:(mem<:((mem:*(b2)),*(b1))),*(b0) :>_
    : ((_,_,_:>_) ~(_<:*(-a1),(mem:*(-a2))));
tf22(b0,b1,b2,a1,a2) = // tf2, direct-form 2:
    _ : (((_,_,_:>_)~*(-a1)<:mem,*(b0))~*(-a2))
      : (_<:mem,*(b1)),_ : *(b2),_,_ :> _;
tf22t(b0,b1,b2,a1,a2) = // tf2, direct-form 2 transposed:
    _ : (_,_,(_ <: *(b2)',*(b1)',*(b0))
      : _,+',_,_ :> _)~*(-a1)~*(-a2) : _;
tf21t(b0,b1,b2,a1,a2) = // tf2, direct-form 1 transposed:
    tf22t(1,0,0,a1,a2) : tf22t(b0,b1,b2,0,0); // or write it out if you want

//=========================== Ladder/Lattice Digital Filters =============================
// Ladder and lattice digital filters generally have superior numerical
// properties relative to direct-form digital filters.  They can be derived
// from digital waveguide filters, which gives them a physical interpretation.

// #### Reference
// 
// * F. Itakura and S. Saito: "Digital Filtering Techniques for Speech Analysis and Synthesis",
//     7th Int. Cong. Acoustics, Budapest, 25 C 1, 1971.
// * J. D. Markel and A. H. Gray: Linear Prediction of Speech, New York: Springer Verlag, 1976.
// * <https://ccrma.stanford.edu/~jos/pasp/Conventional_Ladder_Filters.html>
//========================================================================================

//-------------------------------`av2sv`-----------------------------------
// Compute reflection coefficients sv from transfer-function denominator av.
//
// #### Usage
//
// ```
// sv = av2sv(av)
// ```
//
// Where:
//
// * `av`: parallel signal bank `a1,...,aN`
// * `sv`: parallel signal bank `s1,...,sN`
// 
// where `ro = ith` reflection coefficient, and
//       `ai` = coefficient of `z^(-i)` in the filter
//          transfer-function denominator `A(z)`.
//
// #### Reference
//
//   <https://ccrma.stanford.edu/~jos/filters/Step_Down_Procedure.html> 
//   (where reflection coefficients are denoted by k rather than s).
//------------------------------------------------------------
// TODO: author JOS, revised by RM
av2sv(av) = par(i,M,s(i+1)) with {
  M = ba.count(av);
  s(m) = sr(M-m+1); // m=1..M
  sr(m) = Ari(m,M-m+1); // s_{M-1-m}
  Ari(m,i) = ba.take(i+1,Ar(m-1));
  //step-down recursion for lattice/ladder digital filters:
  Ar(0) = (1,av); // Ar(m) is order M-m (i.e. "reverse-indexed")
  Ar(m) = 1,par(i,M-m, (Ari(m,i+1) - sr(m)*Ari(m,M-m-i))/(1-sr(m)*sr(m)));
};

//----------------------------`bvav2nuv`--------------------------------
// Compute lattice tap coefficients from transfer-function coefficients.
//
// #### Usage
//
// ```
// nuv = bvav2nuv(bv,av)
// ```
//
// Where:
//
// * `av`: parallel signal bank `a1,...,aN`
// * `bv`: parallel signal bank `b0,b1,...,aN`
// * `nuv`: parallel signal bank  `nu1,...,nuN`
//
// where `nui` is the i'th tap coefficient,
//       `bi` is the coefficient of `z^(-i)` in the filter numerator,
//       `ai` is the coefficient of `z^(-i)` in the filter denominator
//------------------------------------------------------------
// TODO: author JOS, revised by RM
bvav2nuv(bv,av) = par(m,M+1,nu(m)) with {
  M = ba.count(av);
  nu(m) = ba.take(m+1,Pr(M-m)); // m=0..M
  // lattice/ladder tap parameters:
  Pr(0) = bv; // Pr(m) is order M-m, 'r' means "reversed"
  Pr(m) = par(i,M-m+1, (Pri(m,i) - nu(M-m+1)*Ari(m,M-m-i+1)));
  Pri(m,i) = ba.take(i+1,Pr(m-1));
  Ari(m,i) = ba.take(i+1,Ar(m-1));
  //step-down recursion for lattice/ladder digital filters:
  Ar(0) = (1,av); // Ar(m) is order M-m (recursion index must start at constant)
  Ar(m) = 1,par(i,M-m, (Ari(m,i+1) - sr(m)*Ari(m,M-m-i))/(1-sr(m)*sr(m)));
  sr(m) = Ari(m,M-m+1); // s_{M-1-m}
};

//--------------------`iir_lat2`-----------------------
// Two-multiply latice IIR filter or arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_lat2(bv,av) : _
// ```
//
// Where:
//
// * bv: zeros as a bank of parallel signals
// * av: poles as a bank of parallel signals
//------------------------------------------------------------
// TODO: author JOS, revised by RM
iir_lat2(bv,av) = allpassnt(M,sv) : sum(i,M+1,*(ba.take(M-i+1,tg)))
with {
  M = ba.count(av);
  sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
  tg = bvav2nuv(bv,av); // tg = vector of tap gains
};

//-----------------------`allpassnt`--------------------------
// Two-multiply lattice allpass (nested order-1 direct-form-ii allpasses).
//
// #### Usage
//
// ```
// _ : allpassnt(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflexion coefficients (-1 1)
//------------------------------------------------------------
// TODO: author JOS, revised by RM
allpassnt(0,sv) = _;
allpassnt(n,sv) =
//0:   x <: ((+ <: (allpassnt(n-1,sv)),*(s))~(*(-s))) : _',_ :+
       _ : ((+ <: (allpassnt(n-1,sv),*(s)))~*(-s)) : fsec(n)
with {
  fsec(1) = ro.crossnn(1) : _, (_<:mem,_) : +,_;
  fsec(n) = ro.crossn1(n) : _, (_<:mem,_),par(i,n-1,_) : +, par(i,n,_);
  innertaps(n) = par(i,n,_);
  s = ba.take(n,sv); // reflection coefficient s = sin(theta)
};

//--------------------`iir_kl`-----------------------
// Kelly-Lochbaum ladder IIR filter or arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_kl(bv,av) : _
// ```
//
// Where:
//
// * bv: zeros as a bank of parallel signals
// * av: poles as a bank of parallel signals
//------------------------------------------------------------
// TODO: author JOS, revised by RM
iir_kl(bv,av) = allpassnklt(M,sv) : sum(i,M+1,*(tghr(i)))
with {
  M = ba.count(av);
  sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
  tg = bvav2nuv(bv,av); // tg = vector of tap gains for 2mul case
  tgr(i) = ba.take(M+1-i,tg);
  tghr(n) = tgr(n)/pi(n);
  pi(0) = 1;
  pi(n) = pi(n-1)*(1+ba.take(M-n+1,sv)); // all sign parameters '+'
};

//-----------------------`allpassnklt`--------------------------
// Kelly-Lochbaum ladder allpass.
//
// #### Usage:
//
// ```
// _ : allpassklt(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflexion coefficients (-1 1)
//------------------------------------------------------------
// TODO: author JOS, revised by RM
allpassnklt(0,sv) = _;
allpassnklt(n,sv) = _ <: *(s),(*(1+s) : (+
                   : allpassnklt(n-1,sv))~(*(-s))) : fsec(n)
with {
  fsec(1) = _, (_<:mem*(1-s),_) : sumandtaps(n);
  fsec(n) = _, (_<:mem*(1-s),_), par(i,n-1,_) : sumandtaps(n);
  s = ba.take(n,sv);
  sumandtaps(n) = +,par(i,n,_);
};

//--------------------`iir_lat1`-----------------------
// One-multiply latice IIR filter or arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_lat1(bv,av) : _
// ```
//
// Where:
//
// * bv: zeros as a bank of parallel signals
// * av: poles as a bank of parallel signals
//------------------------------------------------------------
// TODO: author JOS, revised by RM
iir_lat1(bv,av) = allpassn1mt(M,sv) : sum(i,M+1,*(tghr(i+1)))
with {
  M = ba.count(av);
  sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
  tg = bvav2nuv(bv,av); // tg = vector of tap gains
  tgr(i) = ba.take(M+2-i,tg); // i=1..M+1 (for "takability")
  tghr(n) = tgr(n)/pi(n);
  pi(1) = 1;
  pi(n) = pi(n-1)*(1+ba.take(M-n+2,sv)); // all sign parameters '+'
};

//-----------------------`allpassn1mt`--------------------------
// One-multiply lattice allpass with tap lines.
//
// #### Usage
//
// ```
// _ : allpassn1mt(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflexion coefficients (-1 1)
//------------------------------------------------------------
// TODO: author JOS, revised by RM
allpassn1mt(0,sv) = _;
allpassn1mt(n,sv)= _ <: _,_ : ((+:*(s) <: _,_),_ : _,+ : ro.crossnn(1)
		  : allpassn1mt(n-1,sv),_)~(*(-1)) : fsec(n)
with {
//0:  fsec(n) = _',_ : +
  fsec(1) = ro.crossnn(1) : _, (_<:mem,_) : +,_;
  fsec(n) = ro.crossn1(n) : _, (_<:mem,_),par(i,n-1,_) : +, par(i,n,_);
  innertaps(n) = par(i,n,_);
  s = ba.take(n,sv); // reflection coefficient s = sin(theta)
};

//-------------------------------`iir_nl`-------------------------
// Normalized ladder filter of arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_nl(bv,av) : _
// ```
//
// Where:
//
// * bv: zeros as a bank of parallel signals
// * av: poles as a bank of parallel signals
// 
// #### References
//
// * J. D. Markel and A. H. Gray, Linear Prediction of Speech, New York: Springer Verlag, 1976.
// * <https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
iir_nl(bv,av) = allpassnnlt(M,sv) : sum(i,M+1,*(tghr(i)))
with {
  M = ba.count(av);
  sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
  tg = bvav2nuv(bv,av); // tg = vector of tap gains for 2mul case
  tgr(i) = ba.take(M+1-i,tg);
  tghr(n) = tgr(n)/pi(n);
  pi(0) = 1;
  s(n) = ba.take(M-n+1,sv); // reflection coefficient = sin(theta) 
  c(n) = sqrt(max(0,1-s(n)*s(n))); // compiler crashes on sqrt(-)
  pi(n) = pi(n-1)*c(n);
};

//-------------------------------`allpassnnlt`-------------------------
// Normalized ladder allpass filter of arbitrary order.
//
// #### Usage:
//
// ```
// _ : allpassnnlt(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflexion coefficients (-1,1)
// 
// #### References
//
// * J. D. Markel and A. H. Gray, Linear Prediction of Speech, New York: Springer Verlag, 1976.
// * <https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
allpassnnlt(0,sv) = _;
allpassnnlt(n,scl*(sv)) = allpassnnlt(n,par(i,count(sv),scl*(sv(i))));
allpassnnlt(n,sv) = _ <: *(s),(*(c) : (+
                   : allpassnnlt(n-1,sv))~(*(-s))) : fsec(n)
with {
  fsec(1) = _, (_<:mem*(c),_) : sumandtaps(n);
  fsec(n) = _, (_<:mem*(c),_), par(i,n-1,_) : sumandtaps(n);
  s = ba.take(n,sv);
  c = sqrt(max(0,1-s*s));
  sumandtaps(n) = +,par(i,n,_);
};

//=============================Useful Special Cases=======================================
//========================================================================================

//--------------------------------`tf2np`------------------------------------
// Biquad based on a stable second-order Normalized Ladder Filter
// (more robust to modulation than `tf2` and protected against instability).
//
// #### Usage
//
// ```
// _ : tf2np(b0,b1,b2,a1,a2) : _
// ```
//
// Where:
// 
// * `a`: the poles
// * `b`: the zeros
//------------------------------------------------------------
// TODO: author JOS, revised by RM
tf2np(b0,b1,b2,a1,a2) = allpassnnlt(M,sv) : sum(i,M+1,*(tghr(i)))
with {
  smax = 0.9999; // maximum reflection-coefficient magnitude allowed
  s2 = max(-smax, min(smax,a2)); // Project both reflection-coefficients
  s1 = max(-smax, min(smax,a1/(1+a2))); // into the defined stability-region.
  sv = (s1,s2); // vector of sin(theta) reflection coefficients
  M = 2;
  nu(2) = b2;
  nu(1) = b1 - b2*a1;
  nu(0) = (b0-b2*a2) - nu(1)*s1;
  tg = (nu(0),nu(1),nu(2));
  tgr(i) = ba.take(M+1-i,tg); // vector of tap gains for 2mul case
  tghr(n) = tgr(n)/pi(n);  // apply pi parameters for NLF case
  pi(0) = 1;
  s(n) = ba.take(M-n+1,sv);
  c(n) = sqrt(1-s(n)*s(n));
  pi(n) = pi(n-1)*c(n);
};

//-----------------------------`wgr`---------------------------------
// Second-order transformer-normalized digital waveguide resonator.
//
// #### Usage
//
// ``` 
// _ : wgr(f,r) : _
// ```
//
// Where:
//
// * `f`: resonance frequency (Hz)
// * `r`: loss factor for exponential decay (set to 1 to make a numerically stable oscillator)
//
// #### References
//
// * <https://ccrma.stanford.edu/~jos/pasp/Power_Normalized_Waveguide_Filters.html>
// * <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
wgr(f,r,x) = (*(G),_<:_,((+:*(C))<:_,_),_:+,_,_:+(x),-) ~ cross : _,*(0-gi)
with { 
  C = cos(2*ma.PI*f/ma.SR);
  gi = sqrt(max(0,(1+C)/(1-C))); // compensate amplitude (only needed when
  G = r*(1-1' + gi')/gi;         //   frequency changes substantially)
  cross = _,_ <: !,_,_,!;
};

//-----------------------------`nlf2`--------------------------------
// Second order normalized digital waveguide resonator.
//
// #### Usage 
//
// ```
// _ : nlf2(f,r) : _
// ```
//
// Where:
//
// * `f`: resonance frequency (Hz)
// * `r`: loss factor for exponential decay (set to 1 to make a sinusoidal oscillator)
//
// #### Reference
//
// <https://ccrma.stanford.edu/~jos/pasp/Power_Normalized_Waveguide_Filters.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
nlf2(f,r,x) = ((_<:_,_),(_<:_,_) : (*(s),*(c),*(c),*(0-s)) :> 
              (*(r),+(x))) ~ cross
with { 
  th = 2*ma.PI*f/ma.SR;
  c = cos(th);
  s = sin(th);
  cross = _,_ <: !,_,_,!;
};


//------------`apnl`---------------
// Passive Nonlinear Allpass based on Pierce switching springs idea.
// Switch between allpass coefficient `a1` and `a2` at signal zero crossings.
//
// #### Usage
//
// ```
// _ : apnl(a1,a2) : _
// ```
//
// Where:
//
// * `a1` and `a2`: allpass coefficients
//
// #### Reference
// 
// * "A Passive Nonlinear Digital Filter Design ..." by John R. Pierce and Scott 
// A. Van Duyne, JASA, vol. 101, no. 2, pp. 1120-1126, 1997
//------------------------------------------------------------
// TODO: author JOS and RM
apnl(a1,a2,x) = nonLinFilter
with{
   condition = _>0;
   nonLinFilter = (x - _ <: _*(condition*a1 + (1-condition)*a2),_')~_ :> +;
};


//============================Ladder/Lattice Allpass Filters==============================
// An allpass filter has gain 1 at every frequency, but variable phase.
// Ladder/lattice allpass filters are specified by reflection coefficients.
// They are defined here as nested allpass filters, hence the names allpassn*.
//
// #### References
// 
// * <https://ccrma.stanford.edu/~jos/pasp/Conventional_Ladder_Filters.html>
// * <https://ccrma.stanford.edu/~jos/pasp/Nested_Allpass_Filters.html>
// * Linear Prediction of Speech, Markel and Gray, Springer Verlag, 1976
//========================================================================================

//---------------`allpassn`-----------------
// Two-multiply lattice - each section is two multiply-adds.
//
// #### Usage:
//
// ```
// _ : allpassn(n,sv) : _
// ```
// #### Where:
//
// * `n`: the order of the filter
// * `sv`: the reflexion coefficients (-1 1)
//
// #### References
//
// * J. O. Smith and R. Michon, "Nonlinear Allpass Ladder Filters in FAUST", in 
// Proceedings of the 14th International Conference on Digital Audio Effects 
// (DAFx-11), Paris, France, September 19-23, 2011.
//----------------------------------------------
// TODO: author JOS, revised by RM
allpassn(0,sv) = _;
allpassn(n,sv) = _ <: ((+ <: (allpassn(n-1,sv)),*(s))~(*(-s))) : _',_ :+
with { s = ba.take(n,sv); };

//---------------`allpassnn`-----------------
// Normalized form - four multiplies and two adds per section,
// but coefficients can be time varying and nonlinear without
// "parametric amplification" (modulation of signal energy).
//
// #### Usage:
//
// ```
// _ : allpassnn(n,tv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `tv`: the reflexion coefficients (-PI PI)
//----------------------------------------------
// TODO: author JOS, revised by RM
// power-normalized (reflection coefficients s = sin(t)):
allpassnn(0,tv) = _;
allpassnn(n,tv) = _ <: *(s), (*(c) : (+ 
        : allpassnn(n-1,tv))~(*(-s))) : _, mem*c : +
with { c=cos(ba.take(n,tv));  s=sin(ba.take(n,tv)); };

//---------------`allpasskl`-----------------
// Kelly-Lochbaum form - four multiplies and two adds per
// section, but all signals have an immediate physical
// interpretation as traveling pressure waves, etc.
//
// #### Usage:
//
// ```
// _ : allpassnkl(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflexion coefficients (-1 1)
//----------------------------------------------
// TODO: author JOS, revised by RM
// Kelly-Lochbaum:
allpassnkl(0,sv) = _;
allpassnkl(n,sv) = _ <: *(s),(*(1+s) : (+
                   : allpassnkl(n-1,sv))~(*(-s))) : _, mem*(1-s) : +
with { s = ba.take(n,sv); };

//---------------`allpass1m`-----------------
// One-multiply form - one multiply and three adds per section.
// Normally the most efficient in special-purpose hardware.
//
// #### Usage:
//
// ```
// _ : allpassn1m(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflexion coefficients (-1 1)
//----------------------------------------------
// TODO: author JOS, revised by RM
// one-multiply:
allpassn1m(0,sv) = _;
allpassn1m(n,sv)= _ <: _,_ : ((+:*(s) <: _,_),_ : _,+ : cross
		  : allpassn1m(n-1,sv),_)~(*(-1)) : _',_ : +
with {s = ba.take(n,sv); cross = _,_ <: !,_,_,!; };

//===========Digital Filter Sections Specified as Analog Filter Sections==================
//========================================================================================

//-------------------------`tf2s` and `tf2snp`--------------------------------
// Second-order direct-form digital filter,
// specified by ANALOG transfer-function polynomials B(s)/A(s),
// and a frequency-scaling parameter.  Digitization via the
// bilinear transform is built in.
//
// #### Usage 
//
// ```
// _ : tf2s(b2,b1,b0,a1,a0,w1) : _
// ```
// Where: 
//
// ```
//         b2 s^2 + b1 s + b0
// H(s) = --------------------
//            s^2 + a1 s + a0
// ```
//
// and `w1` is the desired digital frequency (in radians/second)
// corresponding to analog frequency 1 rad/sec (i.e., `s = j`).
//
// #### Example
//
// A second-order ANALOG Butterworth lowpass filter,
//          normalized to have cutoff frequency at 1 rad/sec,
//          has transfer function
//
// ```
//              1
// H(s) = -----------------
//         s^2 + a1 s + 1
// ```
//
// where `a1 = sqrt(2)`.  Therefore, a DIGITAL Butterworth lowpass 
// cutting off at `SR/4` is specified as `tf2s(0,0,1,sqrt(2),1,PI*SR/2);`
//
// #### Method 
//
// Bilinear transform scaled for exact mapping of w1.
//
// #### Reference 
//
// <https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html>
//----------------------------------------------
// TODO: author JOS, revised by RM
tf2s(b2,b1,b0,a1,a0,w1) = tf2(b0d,b1d,b2d,a1d,a2d)
with {
  c   = 1/tan(w1*0.5/ma.SR); // bilinear-transform scale-factor
  csq = c*c;
  d   = a0 + a1 * c + csq;
  b0d = (b0 + b1 * c + b2 * csq)/d;
  b1d = 2 * (b0 - b2 * csq)/d;
  b2d = (b0 - b1 * c + b2 * csq)/d;
  a1d = 2 * (a0 - csq)/d;
  a2d = (a0 - a1*c + csq)/d;
};

// tf2snp = tf2s but using a protected normalized ladder filter for tf2:
tf2snp(b2,b1,b0,a1,a0,w1) = tf2np(b0d,b1d,b2d,a1d,a2d)
with {
  c   = 1/tan(w1*0.5/ma.SR); // bilinear-transform scale-factor
  csq = c*c;
  d   = a0 + a1 * c + csq;
  b0d = (b0 + b1 * c + b2 * csq)/d;
  b1d = 2 * (b0 - b2 * csq)/d;
  b2d = (b0 - b1 * c + b2 * csq)/d;
  a1d = 2 * (a0 - csq)/d;
  a2d = (a0 - a1*c + csq)/d;
};

//-----------------------------`tf3slf`-------------------------------
// Analogous to tf2s above, but third order, and using the typical
// low-frequency-matching bilinear-transform constant 2/T ("lf" series)
// instead of the specific-frequency-matching value used in tf2s and tf1s.
// Note the lack of a "w1" argument.
//
// #### Usage
//
// ```
// _ : tf3slf(b3,b2,b1,b0,a3,a2,a1,a0) : _
// ```
//----------------------------------------------
// TODO: author JOS, revised by RM
tf3slf(b3,b2,b1,b0,a3,a2,a1,a0) = tf3(b0d,b1d,b2d,b3d,a1d,a2d,a3d) with {
  c   = 2.0 * ma.SR; // bilinear-transform scale-factor ("lf" case)
  csq = c*c;
  cc = csq*c;
  // Thank you maxima:
  b3d = (b3*c^3-b2*c^2+b1*c-b0)/d;
  b2d = (-3*b3*c^3+b2*c^2+b1*c-3*b0)/d;
  b1d = (3*b3*c^3+b2*c^2-b1*c-3*b0)/d;
  b0d = (-b3*c^3-b2*c^2-b1*c-b0)/d;
  a3d = (a3*c^3-a2*c^2+a1*c-a0)/d;
  a2d = (-3*a3*c^3+a2*c^2+a1*c-3*a0)/d;
  a1d = (3*a3*c^3+a2*c^2-a1*c-3*a0)/d;
  d = (-a3*c^3-a2*c^2-a1*c-a0);
};

//-----------------------------`tf1s`--------------------------------
// First-order direct-form digital filter,
// specified by ANALOG transfer-function polynomials B(s)/A(s),
// and a frequency-scaling parameter.
//
// #### Usage
//
// ```
// tf1s(b1,b0,a0,w1)
// ```
// Where: 
//
//        b1 s + b0
// H(s) = ----------
//           s + a0
//
// and `w1` is the desired digital frequency (in radians/second)
// corresponding to analog frequency 1 rad/sec (i.e., `s = j`).
//
// #### Example
// A first-order ANALOG Butterworth lowpass filter,
//          normalized to have cutoff frequency at 1 rad/sec,
//          has transfer function
//
//           1
// H(s) = -------
//         s + 1
//
// so `b0 = a0 = 1` and `b1 = 0`.  Therefore, a DIGITAL first-order 
// Butterworth lowpass with gain -3dB at `SR/4` is specified as 
//
// ```
// tf1s(0,1,1,PI*SR/2); // digital half-band order 1 Butterworth
// ```
//
// #### Method
//
// Bilinear transform scaled for exact mapping of w1.
//
// #### Reference 
// <https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html>
//----------------------------------------------
// TODO: author JOS, revised by RM
tf1s(b1,b0,a0,w1) = tf1(b0d,b1d,a1d)
with {
  c   = 1/tan(w1*0.5/ma.SR); // bilinear-transform scale-factor
  d   = a0 + c;
  b1d = (b0 - b1*c) / d;
  b0d = (b0 + b1*c) / d;
  a1d = (a0 - c) / d;
};

//-----------------------------`tf2sb`--------------------------------
// Bandpass mapping of `tf2s`: In addition to a frequency-scaling parameter
// `w1` (set to HALF the desired passband width in rad/sec),
// there is a desired center-frequency parameter wc (also in rad/s).
// Thus, `tf2sb` implements a fourth-order digital bandpass filter section
// specified by the coefficients of a second-order analog lowpass prototpe
// section.  Such sections can be combined in series for higher orders.
// The order of mappings is (1) frequency scaling (to set lowpass cutoff w1),
// (2) bandpass mapping to wc, then (3) the bilinear transform, with the 
// usual scale parameter `2*SR`.  Algebra carried out in maxima and pasted here.
//
// #### Usage
//
// ```
// _ : tf2sb(b2,b1,b0,a1,a0,w1,wc) : _
// ```
//----------------------------------------------
// TODO: author JOS, revised by RM
tf2sb(b2,b1,b0,a1,a0,w1,wc) = 
  iir((b0d/a0d,b1d/a0d,b2d/a0d,b3d/a0d,b4d/a0d),(a1d/a0d,a2d/a0d,a3d/a0d,a4d/a0d)) with {
  T = 1.0/float(ma.SR);
  b0d = (4*b0*w1^2+8*b2*wc^2)*T^2+8*b1*w1*T+16*b2;
  b1d = 4*b2*wc^4*T^4+4*b1*wc^2*w1*T^3-16*b1*w1*T-64*b2;
  b2d = 6*b2*wc^4*T^4+(-8*b0*w1^2-16*b2*wc^2)*T^2+96*b2;
  b3d = 4*b2*wc^4*T^4-4*b1*wc^2*w1*T^3+16*b1*w1*T-64*b2;
  b4d = (b2*wc^4*T^4-2*b1*wc^2*w1*T^3+(4*b0*w1^2+8*b2*wc^2)*T^2-8*b1*w1*T +16*b2)
        + b2*wc^4*T^4+2*b1*wc^2*w1*T^3;
  a0d = wc^4*T^4+2*a1*wc^2*w1*T^3+(4*a0*w1^2+8*wc^2)*T^2+8*a1*w1*T+16;
  a1d = 4*wc^4*T^4+4*a1*wc^2*w1*T^3-16*a1*w1*T-64;
  a2d = 6*wc^4*T^4+(-8*a0*w1^2-16*wc^2)*T^2+96;
  a3d = 4*wc^4*T^4-4*a1*wc^2*w1*T^3+16*a1*w1*T-64;
  a4d = wc^4*T^4-2*a1*wc^2*w1*T^3+(4*a0*w1^2+8*wc^2)*T^2-8*a1*w1*T+16;
};

//-----------------------------`tf1sb`--------------------------------
// First-to-second-order lowpass-to-bandpass section mapping,
// analogous to tf2sb above.
// 
// #### Usage
//
// ```
// _ : tf1sb(b1,b0,a0,w1,wc) : _
// ```
//----------------------------------------------
// TODO: author JOS, revised by RM
tf1sb(b1,b0,a0,w1,wc) = tf2(b0d/a0d,b1d/a0d,b2d/a0d,a1d/a0d,a2d/a0d) with {
  T = 1.0/float(ma.SR);
  a0d = wc^2*T^2+2*a0*w1*T+4;
  b0d = b1*wc^2*T^2 +2*b0*w1*T+4*b1;
  b1d = 2*b1*wc^2*T^2-8*b1;
  b2d = b1*wc^2*T^2-2*b0*w1*T+4*b1;
  a1d = 2*wc^2*T^2-8;
  a2d = wc^2*T^2-2*a0*w1*T+4;
};

//==============================Simple Resonator Filters==================================
//========================================================================================

//------------------`resonlp`-----------------
// Simple resonant lowpass filter based on `tf2s` (virtual analog).
// `resonlp` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : resonlp(fc,Q,gain) : _
// _ : resonhp(fc,Q,gain) : _
// _ : resonbp(fc,Q,gain) : _
// 
// ```
//
// Where:
//
// * `fc`: center frequency (Hz)
// * `Q`: q
// * `gain`: gain (0-1)
//---------------------------------------------------------------------
// TODO: author JOS, revised by RM
// resonlp = 2nd-order lowpass with corner resonance:
resonlp(fc,Q,gain) = tf2s(b2,b1,b0,a1,a0,wc)
with {
     wc = 2*ma.PI*fc;
     a1 = 1/Q;
     a0 = 1;
     b2 = 0;
     b1 = 0;
     b0 = gain;
};


//------------------`resonhp`-----------------
// Simple resonant highpass filters based on `tf2s` (virtual analog).
// `resonhp` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : resonlp(fc,Q,gain) : _
// _ : resonhp(fc,Q,gain) : _
// _ : resonbp(fc,Q,gain) : _
// 
// ```
//
// Where:
//
// * `fc`: center frequency (Hz)
// * `Q`: q
// * `gain`: gain (0-1)
//---------------------------------------------------------------------
// TODO: author JOS, revised by RM
// resonhp = 2nd-order highpass with corner resonance:
resonhp(fc,Q,gain,x) = gain*x-resonlp(fc,Q,gain,x);


//------------------`resonbp`-----------------
// Simple resonant bandpass filters based on `tf2s` (virtual analog).
// `resonbp` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : resonlp(fc,Q,gain) : _
// _ : resonhp(fc,Q,gain) : _
// _ : resonbp(fc,Q,gain) : _
// 
// ```
//
// Where:
//
// * `fc`: center frequency (Hz)
// * `Q`: q
// * `gain`: gain (0-1)
//---------------------------------------------------------------------
// TODO: author JOS, revised by RM
// resonbp = 2nd-order bandpass
resonbp(fc,Q,gain) = tf2s(b2,b1,b0,a1,a0,wc) 
with {
     wc = 2*ma.PI*fc;
     a1 = 1/Q;
     a0 = 1;
     b2 = 0;
     b1 = gain;
     b0 = 0;
};

//======================Butterworth Lowpass/Highpass Filters==============================
//========================================================================================

//----------------`lowpass`--------------------
// Nth-order Butterworth lowpass filter.
// `lowpass` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : lowpass(N,fc) : _
// ```
//
// Where:
//
// * `N`: filter order (number of poles) [nonnegative constant integer]
// * `fc`: desired cut-off frequency (-3dB frequency) in Hz
//
// #### References
// 
// * <https://ccrma.stanford.edu/~jos/filters/Butterworth_Lowpass_Design.html>
// * `butter` function in Octave `("[z,p,g] = butter(N,1,'s');")`
//------------------------------
// TODO: author JOS and Orlarey, revised by RM
lowpass(N,fc) = lowpass0_highpass1(0,N,fc);


//----------------`highpass`--------------------
// Nth-order Butterworth highpass filters.
// `highpass` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : highpass(N,fc) : _
// ```
//
// Where:
//
// * `N`: filter order (number of poles) [nonnegative constant integer]
// * `fc`: desired cut-off frequency (-3dB frequency) in Hz
//
// #### References
// 
// * <https://ccrma.stanford.edu/~jos/filters/Butterworth_Lowpass_Design.html>
// * `butter` function in Octave `("[z,p,g] = butter(N,1,'s');")`
//------------------------------
// TODO: author JOS and Orlarey, revised by RM
highpass(N,fc) = lowpass0_highpass1(1,N,fc);


//-------------`lowpass0_highpass1`--------------
// TODO
//------------------------------
// TODO: author JOS and Orlarey
lowpass0_highpass1(s,N,fc) = lphpr(s,N,N,fc)
with {
  lphpr(s,0,N,fc) = _;
  lphpr(s,1,N,fc) = tf1s(s,1-s,1,2*ma.PI*fc);
  lphpr(s,O,N,fc) = lphpr(s,(O-2),N,fc) : tf2s(s,0,1-s,a1s,1,w1) with {
    parity = N % 2;
    S = (O-parity)/2; // current section number
    a1s = -2*cos((ma.PI)*-1 + (1-parity)*ma.PI/(2*N) + (S-1+parity)*ma.PI/N);
    w1 = 2*ma.PI*fc;
  };
};

//================Special Filter-Bank Delay-Equalizing Allpass Filters====================
// These special allpass filters are needed by filterbank et al. below.
// They are equivalent to (`lowpass(N,fc)` +|- `highpass(N,fc))/2`, but with 
// canceling pole-zero pairs removed (which occurs for odd N).
//========================================================================================

//--------------------`lowpass_plus`|`minus_highpass`----------------
// TODO
//-------------------------------------------------------------------
// TODO: author JOS, revised by RM
highpass_plus_lowpass(1,fc) = _;
highpass_plus_lowpass(3,fc) = tf2s(1,-1,1,1,1,w1) with { w1 = 2*ma.PI*fc; };
highpass_minus_lowpass(3,fc) = tf1s(-1,1,1,w1) with { w1 = 2*ma.PI*fc; };
highpass_plus_lowpass(5,fc) = tf2s(1,-a11,1,a11,1,w1)
with { 
  a11 = 1.618033988749895;
  w1 = 2*ma.PI*fc; 
};
highpass_minus_lowpass(5,fc) = tf1s(1,-1,1,w1) : tf2s(1,-a12,1,a12,1,w1)
with { 
  a12 = 0.618033988749895;
  w1 = 2*ma.PI*fc; 
};

// Catch-all definitions for generality - even order is done:

highpass_plus_lowpass(N,fc) = switch_odd_even(N%2,N,fc) with {
  switch_odd_even(0,N,fc) = highpass_plus_lowpass_even(N,fc);
  switch_odd_even(1,N,fc) = highpass_plus_lowpass_odd(N,fc);
};

highpass_minus_lowpass(N,fc) = switch_odd_even(N%2,N,fc) with {
  switch_odd_even(0,N,fc) = highpass_minus_lowpass_even(N,fc);
  switch_odd_even(1,N,fc) = highpass_minus_lowpass_odd(N,fc);
};

highpass_plus_lowpass_even(N,fc) = highpass(N,fc) + lowpass(N,fc);
highpass_minus_lowpass_even(N,fc) = highpass(N,fc) - lowpass(N,fc);

// FIXME: Rewrite the following, as for orders 3 and 5 above,
//        to eliminate pole-zero cancellations:
highpass_plus_lowpass_odd(N,fc) = highpass(N,fc) + lowpass(N,fc);
highpass_minus_lowpass_odd(N,fc) = highpass(N,fc) - lowpass(N,fc);

//==========================Elliptic (Cauer) Lowpass Filters==============================
// Elliptic (Cauer) Lowpass Filters
//
// #### References
//
// * <http://en.wikipedia.org/wiki/Elliptic_filter
// * functions `ncauer` and `ellip` in Octave
//========================================================================================

//-----------------------------`lowpass3e`-----------------------------
// Third-order Elliptic (Cauer) lowpass filter.
//
// #### Usage
//
// ```
// _ : lowpass3e(fc) : _
// ```
//
// Where: 
//
// * `fc`: -3dB frequency in Hz
//
// #### Design
// 
// For spectral band-slice level display (see `octave_analyzer3e`):
//
// ```
// [z,p,g] = ncauer(Rp,Rs,3);  % analog zeros, poles, and gain, where
// Rp = 60  % dB ripple in stopband
// Rs = 0.2 % dB ripple in passband
// ```
//---------------------------------------------------------------------
// TODO: author JOS, revised by RM
lowpass3e(fc) = tf2s(b21,b11,b01,a11,a01,w1) : tf1s(0,1,a02,w1)
with { 
  a11 = 0.802636764161030; // format long; poly(p(1:2)) % in octave
  a01 = 1.412270893774204;
  a02 = 0.822445908998816; // poly(p(3)) % in octave
  b21 = 0.019809144837789; // poly(z)
  b11 = 0;
  b01 = 1.161516418982696;
  w1 = 2*ma.PI*fc; 
};

//-----------------------------`lowpass6e`-----------------------------
// Sixth-order Elliptic/Cauer lowpass filter.
//
// #### Usage
//
// ```
// _ : lowpass6e(fc) : _
// ```
//
// Where: 
//
// * `fc`: -3dB frequency in Hz
//
// #### Design 
//
// For spectral band-slice level display (see octave_analyzer6e):
//
// ```
// [z,p,g] = ncauer(Rp,Rs,6);  % analog zeros, poles, and gain, where
//  Rp = 80  % dB ripple in stopband
//  Rs = 0.2 % dB ripple in passband
// ```
//----------------------------------------------------------------------
// TODO: author JOS, revised by RM
lowpass6e(fc) = 
              tf2s(b21,b11,b01,a11,a01,w1) :
              tf2s(b22,b12,b02,a12,a02,w1) :
              tf2s(b23,b13,b03,a13,a03,w1)
with { 
  b21 = 0.000099999997055;
  a21 = 1;
  b11 = 0;
  a11 = 0.782413046821645;
  b01 = 0.000433227200555;
  a01 = 0.245291508706160;
  b22 = 1;
  a22 = 1;
  b12 = 0;
  a12 = 0.512478641889141;
  b02 = 7.621731298870603;
  a02 = 0.689621364484675;
  b23 = 1;
  a23 = 1;
  b13 = 0;
  a13 = 0.168404871113589;
  b03 = 53.536152954556727;
  a03 = 1.069358407707312;
  w1 = 2*ma.PI*fc; 
};

//=========================Elliptic Highpass Filters======================================
//========================================================================================

//-----------------------------`highpass3e`-----------------------------
// Third-order Elliptic (Cauer) highpass filter. Inversion of `lowpass3e` wrt unit 
// circle in s plane (s <- 1/s)
//
// #### Usage
//
// ```
// _ : highpass3e(fc) : _
// ```
//
// Where: 
//
// * `fc`: -3dB frequency in Hz
//-------------------------------------------------------------------------
// TODO: author JOS, revised by RM
highpass3e(fc) = tf2s(b01/a01,b11/a01,b21/a01,a11/a01,1/a01,w1) : 
                 tf1s(1/a02,0,1/a02,w1)
with { 
  a11 = 0.802636764161030;
  a01 = 1.412270893774204;
  a02 = 0.822445908998816;
  b21 = 0.019809144837789;
  b11 = 0;
  b01 = 1.161516418982696;
  w1 = 2*ma.PI*fc; 
};

//-----------------------------`highpass6e`-----------------------------
// Sixth-order Elliptic/Cauer highpass filter. Inversion of lowpass3e wrt unit 
// circle in s plane (s <- 1/s)
//
// #### Usage
//
// ```
// _ : highpass6e(fc) : _
// ```
//
// Where: 
//
// * `fc`: -3dB frequency in Hz
//-------------------------------------------------------------------------
// TODO: author JOS, revised by RM
highpass6e(fc) = 
              tf2s(b01/a01,b11/a01,b21/a01,a11/a01,1/a01,w1) :
              tf2s(b02/a02,b12/a02,b22/a02,a12/a02,1/a02,w1) :
              tf2s(b03/a03,b13/a03,b23/a03,a13/a03,1/a03,w1)
with { 
  b21 = 0.000099999997055;
  a21 = 1;
  b11 = 0;
  a11 = 0.782413046821645;
  b01 = 0.000433227200555;
  a01 = 0.245291508706160;
  b22 = 1;
  a22 = 1;
  b12 = 0;
  a12 = 0.512478641889141;
  b02 = 7.621731298870603;
  a02 = 0.689621364484675;
  b23 = 1;
  a23 = 1;
  b13 = 0;
  a13 = 0.168404871113589;
  b03 = 53.536152954556727;
  a03 = 1.069358407707312;
  w1 = 2*ma.PI*fc; 
};

//========================Butterworth Bandpass/Bandstop Filters===========================
//========================================================================================

//--------------------`bandpass`----------------
// Order 2*Nh Butterworth bandpass filter made using the transformation
// `s <- s + wc^2/s` on `lowpass(Nh)`, where `wc` is the desired bandpass center 
// frequency.  The `lowpass(Nh)` cutoff `w1` is half the desired bandpass width.
// `bandpass` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : bandpass(Nh,fl,fu) : _
// ```
//
// Where:
// 
// * `Nh`: HALF the desired bandpass order (which is therefore even)
// * `fl`: lower -3dB frequency in Hz
// * `fu`: upper -3dB frequency in Hz
// Thus, the passband width is `fu-fl`, 
//       and its center frequency is `(fl+fu)/2`.
//
// #### Reference
//
// <http://cnx.org/content/m16913/latest/>
//-------------------------------------------------------------------------
// TODO: author JOS, revised by RM
bandpass(Nh,fl,fu) = bandpass0_bandstop1(0,Nh,fl,fu);


//--------------------`bandstop`----------------
// Order 2*Nh Butterworth bandstop filter made using the transformation
// `s <- s + wc^2/s` on `highpass(Nh)`, where `wc` is the desired bandpass center 
// frequency.  The `highpass(Nh)` cutoff `w1` is half the desired bandpass width.
// `bandstop` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : bandstop(Nh,fl,fu) : _
// ```
// Where:
// 
// * `Nh`: HALF the desired bandstop order (which is therefore even)
// * `fl`: lower -3dB frequency in Hz
// * `fu`: upper -3dB frequency in Hz
// Thus, the passband (stopband) width is `fu-fl`, 
//       and its center frequency is `(fl+fu)/2`.
//
// #### Reference
//
// <http://cnx.org/content/m16913/latest/>
//-------------------------------------------------------------------------
// TODO: author JOS, revised by RM
bandstop(Nh,fl,fu) = bandpass0_bandstop1(1,Nh,fl,fu);
bandpass0_bandstop1(s,Nh,fl,fu) = bpbsr(s,Nh,Nh,fl,fu) 
with {
  wl = 2*ma.PI*fl; // digital (z-plane) lower passband edge
  wu = 2*ma.PI*fu; // digital (z-plane) upper passband edge

  c = 2.0*ma.SR; // bilinear transform scaling used in tf2sb, tf1sb
  wla = c*tan(wl/c); // analog (s-splane) lower cutoff
  wua = c*tan(wu/c); // analog (s-splane) upper cutoff

  wc = sqrt(wla*wua); // s-plane center frequency
  w1 = wua - wc^2/wua; // s-plane lowpass prototype cutoff

  bpbsr(s,0,Nh,fl,fu) = _;
  bpbsr(s,1,Nh,fl,fu) = tf1sb(s,1-s,1,w1,wc);
  bpbsr(s,O,Nh,fl,fu) = bpbsr(s,O-2,Nh,fl,fu) : tf2sb(s,0,(1-s),a1s,1,w1,wc)
  with {
    parity = Nh % 2;
    S = (O-parity)/2; // current section number
    a1s = -2*cos(-1*ma.PI + (1-parity)*ma.PI/(2*Nh) + (S-1+parity)*ma.PI/Nh);
  };
};

//===========================Elliptic Bandpass Filters====================================
//========================================================================================

//---------------------`bandpass6e`-----------------------------
// Order 12 elliptic bandpass filter analogous to `bandpass(6)`.
//--------------------------------------------------------------
// TODO: author JOS, revised by RM
bandpass6e(fl,fu) = tf2sb(b21,b11,b01,a11,a01,w1,wc) : tf1sb(0,1,a02,w1,wc)
with { 
  a11 = 0.802636764161030; // In octave: format long; poly(p(1:2))
  a01 = 1.412270893774204;
  a02 = 0.822445908998816; // poly(p(3))
  b21 = 0.019809144837789; // poly(z)
  b11 = 0;
  b01 = 1.161516418982696;

  wl = 2*ma.PI*fl; // digital (z-plane) lower passband edge
  wu = 2*ma.PI*fu; // digital (z-plane) upper passband edge

  c = 2.0*ma.SR; // bilinear transform scaling used in tf2sb, tf1sb
  wla = c*tan(wl/c); // analog (s-splane) lower cutoff
  wua = c*tan(wu/c); // analog (s-splane) upper cutoff

  wc = sqrt(wla*wua); // s-plane center frequency
  w1 = wua - wc^2/wua; // s-plane lowpass cutoff
};

//----------------------`bandpass12e`---------------------------
// Order 24 elliptic bandpass filter analogous to `bandpass(6)`.
//--------------------------------------------------------------
// TODO: author JOS, revised by RM
bandpass12e(fl,fu) = 
              tf2sb(b21,b11,b01,a11,a01,w1,wc) :
              tf2sb(b22,b12,b02,a12,a02,w1,wc) :
              tf2sb(b23,b13,b03,a13,a03,w1,wc)
with { // octave script output:
  b21 = 0.000099999997055;
  a21 = 1;
  b11 = 0;
  a11 = 0.782413046821645;
  b01 = 0.000433227200555;
  a01 = 0.245291508706160;
  b22 = 1;
  a22 = 1;
  b12 = 0;
  a12 = 0.512478641889141;
  b02 = 7.621731298870603;
  a02 = 0.689621364484675;
  b23 = 1;
  a23 = 1;
  b13 = 0;
  a13 = 0.168404871113589;
  b03 = 53.536152954556727;
  a03 = 1.069358407707312;

  wl = 2*ma.PI*fl; // digital (z-plane) lower passband edge
  wu = 2*ma.PI*fu; // digital (z-plane) upper passband edge

  c = 2.0*ma.SR; // bilinear transform scaling used in tf2sb, tf1sb
  wla = c*tan(wl/c); // analog (s-splane) lower cutoff
  wua = c*tan(wu/c); // analog (s-splane) upper cutoff

  wc = sqrt(wla*wua); // s-plane center frequency
  w1 = wua - wc^2/wua; // s-plane lowpass cutoff
};

//=================Parametric Equalizers (Shelf, Peaking)=================================
// Parametric Equalizers (Shelf, Peaking)
//
// #### References
//
// * <http://en.wikipedia.org/wiki/Equalization>
// * <http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt>
// * Digital Audio Signal Processing, Udo Zolzer, Wiley, 1999, p. 124
// * https://ccrma.stanford.edu/~jos/filters/Low_High_Shelving_Filters.html>
// * https://ccrma.stanford.edu/~jos/filters/Peaking_Equalizers.html>
// * maxmsp.lib in the Faust distribution
// * bandfilter.dsp in the faust2pd distribution 
//========================================================================================

//----------------------`low_shelf`----------------------
// First-order "low shelf" filter (gain boost|cut between dc and some frequency)
// `low_shelf` is a standard Faust function.
//
// #### Usage
//
// ``` 
// _ : lowshelf(N,L0,fx) : _
// _ : low_shelf(L0,fx) : _ // default case (order 3)
// _ : lowshelf_other_freq(N,L0,fx) : _
// ```
//
// Where:
// * `N`: filter order 1, 3, 5, ... (odd only). (default should be 3)
// * `L0`: desired level (dB) between dc and fx (boost `L0>0` or cut `L0<0`)
// * `fx`: -3dB frequency of lowpass band (`L0>0`) or upper band (`L0<0`)
//       (see "SHELF SHAPE" below).
//
// The gain at SR/2 is constrained to be 1.
// The generalization to arbitrary odd orders is based on the well known
// fact that odd-order Butterworth band-splits are allpass-complementary
// (see filterbank documentation below for references).
//
// #### Shelf Shape
// The magnitude frequency response is approximately piecewise-linear
// on a log-log plot ("BODE PLOT").  The Bode "stick diagram" approximation
// L(lf) is easy to state in dB versus dB-frequency lf = dB(f):
//
// * L0 > 0:
// 	* L(lf) = L0, f between 0 and fx = 1st corner frequency;
// 	* L(lf) = L0 - N * (lf - lfx), f between fx and f2 = 2nd corner frequency;
// 	* L(lf) = 0, lf > lf2.
// 	* lf2 = lfx + L0/N = dB-frequency at which level gets back to 0 dB.
// * L0 < 0:
// 	* L(lf) = L0, f between 0 and f1 = 1st corner frequency;
// 	* L(lf) = - N * (lfx - lf), f between f1 and lfx = 2nd corner frequency;
// 	* L(lf) = 0, lf > lfx.
// 	* lf1 = lfx + L0/N = dB-frequency at which level goes up from L0.
//
//  See `lowshelf_other_freq`.
//--------------------------------------------------------------
// TODO: author JOS, revised by RM
lowshelf(N,L0,fx) = filterbank(N,(fx)) : _, *(ba.db2linear(L0)) :> _;
// Special cases and optimization:
low_shelf  = lowshelf(3); // default = 3rd order Butterworth
low_shelf1(L0,fx,x) = x + (ba.db2linear(L0)-1)*lowpass(1,fx,x); // optimized
low_shelf1_l(G0,fx,x) = x + (G0-1)*lowpass(1,fx,x); // optimized
lowshelf_other_freq(N, L0, fx) = ba.db2linear(ba.linear2db(fx) + L0/N); // convenience

//-------------`high_shelf`--------------
// First-order "high shelf" filter (gain boost|cut above some frequency).
// `high_shelf` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : highshelf(N,Lpi,fx) : _
// _ : high_shelf(L0,fx) : _ // default case (order 3)
// _ : highshelf_other_freq(N,Lpi,fx) : _
// ```
//
// Where:
//
// * `N`: filter order 1, 3, 5, ... (odd only).
// * `Lpi`: desired level (dB) between fx and SR/2 (boost Lpi>0 or cut Lpi<0)
// * `fx`: -3dB frequency of highpass band (L0>0) or lower band (L0<0)
//        (Use highshelf_other_freq() below to find the other one.)
//
// The gain at dc is constrained to be 1.
// See `lowshelf` documentation above for more details on shelf shape.
//--------------------------------------------------------------
// TODO: author JOS, revised by RM
highshelf(N,Lpi,fx) = filterbank(N,(fx)) : *(ba.db2linear(Lpi)), _ :> _;
// Special cases and optimization:
high_shelf = highshelf(3); // default = 3rd order Butterworth
high_shelf1(Lpi,fx,x) = x + (ba.db2linear(Lpi)-1)*highpass(1,fx,x); // optimized
high_shelf1_l(Gpi,fx,x) = x + (Gpi-1)*highpass(1,fx,x); //optimized

// shelf transitions between frequency fx and this one:
highshelf_other_freq(N, Lpi, fx) = ba.db2linear(ba.linear2db(fx) - Lpi/N);


//-------------------`peak_eq`------------------------------
// Second order "peaking equalizer" section (gain boost or cut near some frequency)
// Also called a "parametric equalizer" section.
// `peak_eq` is a standard Faust function.
//
// #### Usage 
//
// ```
// _ : peak_eq(Lfx,fx,B) : _;
// ```
//
// Where:
//
// * `Lfx`: level (dB) at fx (boost Lfx>0 or cut Lfx<0)
// * `fx`: peak frequency (Hz)
// * `B`: bandwidth (B) of peak in Hz
//--------------------------------------------------------------
// TODO: author JOS, revised by RM
peak_eq(Lfx,fx,B) = tf2s(1,b1s,1,a1s,1,wx) with {
  T = float(1.0/ma.SR);
  Bw = B*T/sin(wx*T); // prewarp s-bandwidth for more accuracy in z-plane
  a1 = ma.PI*Bw;
  b1 = g*a1;
  g = ba.db2linear(abs(Lfx));
  b1s = select2(Lfx>0,a1,b1); // When Lfx>0, pole dominates bandwidth
  a1s = select2(Lfx>0,b1,a1); // When Lfx<0, zero dominates
  wx = 2*ma.PI*fx;
};

//--------------------`peak_eq_cq`----------------------------
// Constant-Q second order peaking equalizer section. 
//
// #### Usage
//
// ```
// _ : peak_eq_cq(Lfx,fx,Q) : _;
// ```
//
// Where:
//
// * `Lfx`: level (dB) at fx
// * `fx`: boost or cut frequency (Hz)
// * `Q`: "Quality factor" = fx/B where B = bandwidth of peak in Hz
//------------------------------------------------------------
// TODO: author JOS, revised by RM
peak_eq_cq(Lfx,fx,Q) = peak_eq(Lfx,fx,fx/Q);

//-------------------`peak_eq_rm`--------------------------
// Regalia-Mitra second order peaking equalizer section
//
// #### Usage
//
// ```
// _ : peak_eq_rm(Lfx,fx,tanPiBT) : _;
// ```
//
// Where:
//
// * `Lfx`: level (dB) at fx
// * `fx`: boost or cut frequency (Hz)
// * `tanPiBT`: `tan(PI*B/SR)`, where B = -3dB bandwidth (Hz) when 10^(Lfx/20) = 0
//         ~ PI*B/SR for narrow bandwidths B
//
// #### Reference 
//
// P.A. Regalia, S.K. Mitra, and P.P. Vaidyanathan,
// "The Digital All-Pass Filter: A Versatile Signal Processing Building Block"
// Proceedings of the IEEE, 76(1):19-37, Jan. 1988.  (See pp. 29-30.)
//------------------------------------------------------------
// TODO: author JOS, revised by RM
peak_eq_rm(Lfx,fx,tanPiBT) = _ <: _,A,_ : +,- : *(0.5),*(K/2.0) : + with {
  A = tf2(k2, k1*(1+k2), 1, k1*(1+k2), k2) <: _,_; // allpass
  k1 = 0.0 - cos(2.0*ma.PI*fx/ma.SR);
  k2 = (1.0 - tanPiBT)/(1.0 + tanPiBT);
  K = ba.db2linear(Lfx);
};


//---------------------`spectral_tilt`-------------------------
// Spectral tilt filter, providing an arbitrary spectral rolloff factor
// alpha in (-1,1), where
//  -1 corresponds to one pole (-6 dB per octave), and
//  +1 corresponds to one zero (+6 dB per octave).
// In other words, alpha is the slope of the ln magnitude versus ln frequency.
// For a "pinking filter" (e.g., to generate 1/f noise from white noise),
// set alpha to -1/2.
//
// #### Usage
//
// ```
// _ : spectral_tilt(N,f0,bw,alpha) : _
// ```
// Where:
//
// * `N`: desired integer filter order (fixed at compile time)
// * `f0`: lower frequency limit for desired roll-off band
// * `bw`: bandwidth of desired roll-off band
// * `alpha`: slope of roll-off desired in nepers per neper
//         (ln mag / ln radian freq)
//
// #### Examples
// See `spectral_tilt_demo`.
//
// #### Reference
// Link to appear here when write up is done
//------------------------------------------------------------
// TODO: author JOS, revised by RM
spectral_tilt(N,f0,bw,alpha) = seq(i,N,sec(i)) with {
  sec(i) = g * tf1s(b1,b0,a0,1) with {
    g = a0/b0; // unity dc-gain scaling
    b1 = 1.0;
    b0 = mzh(i);
    a0 = mph(i);
    mzh(i) = prewarp(mz(i),ma.SR,w0); // prewarping for bilinear transform
    mph(i) = prewarp(mp(i),ma.SR,w0);
    prewarp(w,SR,wp) = wp * tan(w*T/2) / tan(wp*T/2) with { T = 1/ma.SR; };
    mz(i) = w0 * r ^ (-alpha+i); // minus zero i in s plane
    mp(i) = w0 * r ^ i; // minus pole i in s plane
    w0 = 2 * ma.PI * f0; // radian frequency of first pole
    f1 = f0 + bw; // upper band limit
    r = (f1/f0)^(1.0/float(N-1)); // pole ratio (2 => octave spacing)
  };
};


//----------------------`levelfilter`----------------------
// Dynamic level lowpass filter.
// `levelfilter` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : levelfilter(L,freq) : _
// ```
//
// Where:
//
// * `L`: desired level (in dB) at Nyquist limit (SR/2), e.g., -60
// * `freq`: corner frequency (-3dB point) usually set to fundamental freq
// * `N`: Number of filters in series where L = L/N
//
// #### Reference
//
// <https://ccrma.stanford.edu/realsimple/faust_strings/Dynamic_Level_Lowpass_Filter.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
levelfilter(L,freq,x) = (L * L0 * x) + ((1.0-L) * lp2out(x))
with {
  L0 = pow(L,1/3);
  Lw = ma.PI*freq/ma.SR; // = w1 T / 2
  Lgain = Lw / (1.0 + Lw);
  Lpole2 = (1.0 - Lw) / (1.0 + Lw);
  lp2out = *(Lgain) : + ~ *(Lpole2);
};


//----------------------`levelfilterN`----------------------
// Dynamic level lowpass filter.
//
// #### Usage
//
// ```
// _ : levelfilterN(N,freq,L) : _
// ```
//
// Where:
//
// * `L`: desired level (in dB) at Nyquist limit (SR/2), e.g., -60
// * `freq`: corner frequency (-3dB point) usually set to fundamental freq
// * `N`: Number of filters in series where L = L/N
//
// #### Reference
//
// <https://ccrma.stanford.edu/realsimple/faust_strings/Dynamic_Level_Lowpass_Filter.html>
//------------------------------------------------------------
// TODO: author JOS, revised by RM
levelfilterN(N,freq,L) = seq(i,N,levelfilter((L/N),freq));


//=================================Mth-Octave Filter-Banks================================
// Mth-octave filter-banks split the input signal into a bank of parallel signals, one 
// for each spectral band. They are related to the Mth-Octave Spectrum-Analyzers in 
// `analysis.lib`.
// The documentation of this library contains more details about the implementation.
// The parameters are:
//
// * `M`: number of band-slices per octave (>1)
// * `N`: total number of bands (>2)
// * `ftop`: upper bandlimit of the Mth-octave bands (<SR/2)
//
// In addition to the Mth-octave output signals, there is a highpass signal
// containing frequencies from ftop to SR/2, and a "dc band" lowpass signal 
// containing frequencies from 0 (dc) up to the start of the Mth-octave bands.
// Thus, the N output signals are
//
// ```
// highpass(ftop), MthOctaveBands(M,N-2,ftop), dcBand(ftop*2^(-M*(N-1)))
// ```
//
// A Filter-Bank is defined here as a signal bandsplitter having the
// property that summing its output signals gives an allpass-filtered
// version of the filter-bank input signal.  A more conventional term for
// this is an "allpass-complementary filter bank".  If the allpass filter
// is a pure delay (and possible scaling), the filter bank is said to be
// a "perfect-reconstruction filter bank" (see Vaidyanathan-1993 cited
// below for details).  A "graphic equalizer", in which band signals
// are scaled by gains and summed, should be based on a filter bank.
//
// The filter-banks below are implemented as Butterworth or Elliptic
// spectrum-analyzers followed by delay equalizers that make them 
// allpass-complementary.
//
// #### Increasing Channel Isolation
//
// Go to higher filter orders - see Regalia et al. or Vaidyanathan (cited 
// below) regarding the construction of more aggressive recursive 
// filter-banks using elliptic or Chebyshev prototype filters.
//   
// #### References
//
// * "Tree-structured complementary filter banks using all-pass sections",
//   Regalia et al., IEEE Trans. Circuits & Systems, CAS-34:1470-1484, Dec. 1987
// * "Multirate Systems and Filter Banks", P. Vaidyanathan, Prentice-Hall, 1993
// * Elementary filter theory: https://ccrma.stanford.edu/~jos/filters/
//========================================================================================

//------------------------`mth_octave_filterbank[n]`-------------------------
// Allpass-complementary filter banks based on Butterworth band-splitting.
// For Butterworth band-splits, the needed delay equalizer is easily found.
//
// #### Usage
//
// ```
// _ : mth_octave_filterbank(O,M,ftop,N) : par(i,N,_); // Oth-order
// _ : mth_octave_filterbank_alt(O,M,ftop,N) : par(i,N,_); // dc-inverted version
// ```
//
// Also for convenience:
//
// ```
// _ : mth_octave_filterbank3(M,ftop,N) : par(i,N,_); // 3d-order Butterworth
// _ : mth_octave_filterbank5(M,ftop,N) : par(i,N,_); // 5th-roder Butterworth
// mth_octave_filterbank_default = mth_octave_analyzer6e;
// ```
//
// Where: 
//
// * `O`: order of filter used to split each frequency band into two
// * `M`: number of band-slices per octave
// * `ftop`: highest band-split crossover frequency (e.g., 20 kHz)
// * `N`: total number of bands (including dc and Nyquist)
//------------------------------------------------------------
// TODO: author JOS and Orlarey, revised by RM 
mth_octave_filterbank(O,M,ftop,N) = 
    an.mth_octave_analyzer(O,M,ftop,N) : 
    delayeq(N) with {
  fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
  ap(n) = highpass_plus_lowpass(O,fc(n));   // delay-equalizing allpass
  delayeq(N) = par(i,N-2,apchain(i+1)), _, _;
  apchain(i) = seq(j,N-1-i,ap(j+1));
};

// dc-inverted version.  This reduces the delay-equalizer order for odd O.
// Negating the input signal makes the dc band noninverting
//   and all higher bands sign-inverted (if preferred).
mth_octave_filterbank_alt(O,M,ftop,N) = 
    an.mth_octave_analyzer(O,M,ftop,N) : delayeqi(O,N) with {
  fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
  ap(n) = highpass_minus_lowpass(O,fc(n)); // half the order of 'plus' case
  delayeqi(N) = par(i,N-2,apchain(i+1)), _, *(-1.0);
  apchain(i) = seq(j,N-1-i,ap(j+1));
};

// Note that even-order cases require complex coefficients.
// See Vaidyanathan 1993 and papers cited there for more info.
mth_octave_filterbank3(M,ftop,N) = mth_octave_filterbank_alt(3,M,ftop,N);
mth_octave_filterbank5(M,ftop,N) = mth_octave_filterbank(5,M,ftop,N);

mth_octave_filterbank_default = mth_octave_filterbank5;


//===============Arbritary-Crossover Filter-Banks and Spectrum Analyzers==================
// These are similar to the Mth-octave analyzers above, except that the
// band-split frequencies are passed explicitly as arguments. 
//========================================================================================

// ACKNOWLEDGMENT
// Technique for processing a variable number of signal arguments due
// to Yann Orlarey (as is the entire Faust framework!)

//---------------`filterbank`--------------------------
// Filter bank.
// `filterbank` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : filterbank (O,freqs) : par(i,N,_); // Butterworth band-splits
// ```
// Where: 
//
// * `O`: band-split filter order (ODD integer required for filterbank[i])
// * `freqs`: (fc1,fc2,...,fcNs) [in numerically ascending order], where
//           Ns=N-1 is the number of octave band-splits 
//           (total number of bands N=Ns+1). 
// 
// If frequencies are listed explicitly as arguments, enclose them in parens:
//
// ```
// _ : filterbank(3,(fc1,fc2)) : _,_,_
// ```
//---------------------------------------------------
// TODO: author JOS, revised by RM
filterbank(O,lfreqs) = an.analyzer(O,lfreqs) : delayeq(nb) with
{
   nb = ba.count(lfreqs);
   fc(n) = ba.take(n, lfreqs);
   ap(n) = highpass_plus_lowpass(O,fc(n));
   delayeq(1) = _,_; // par(i,0,...) does not fly
   delayeq(nb) = par(i,nb-1,apchain(nb-1-i)),_,_;
   apchain(0) = _;
   apchain(i) =  ap(i) : apchain(i-1);
};

//-----------------`filterbanki`----------------------
// Inverted-dc filter bank.
//
// #### Usage
//
// ```
// _ : filterbanki(O,freqs) : par(i,N,_); // Inverted-dc version 
// ```
//
// Where: 
//
// * `O`: band-split filter order (ODD integer required for `filterbank[i]`)
// * `freqs`: (fc1,fc2,...,fcNs) [in numerically ascending order], where
//           Ns=N-1 is the number of octave band-splits 
//           (total number of bands N=Ns+1). 
// 
// If frequencies are listed explicitly as arguments, enclose them in parens:
//
// ```
// _ : filterbanki(3,(fc1,fc2)) : _,_,_
// ```
//---------------------------------------------------
// TODO: author JOS, revised by RM 
filterbanki(O,lfreqs) = _ <: bsplit(nb) with
{
   nb = ba.count(lfreqs);
   fc(n) = ba.take(n, lfreqs);
   lp(n) = lowpass(O,fc(n));
   hp(n) = highpass(O,fc(n));
   ap(n) = highpass_minus_lowpass(O,fc(n));
   bsplit(0) = *(-1.0);
   bsplit(i) = (hp(i) : delayeq(i-1)), (lp(i) <: bsplit(i-1));
   delayeq(0) = _; // moving the *(-1) here inverts all outputs BUT dc
   delayeq(i) =  ap(i) : delayeq(i-1);
};