filters.lib

//##################################### filters.lib ########################################
// Faust Filters library; Its official prefix is `fi`.
//
// The Filters library is organized into 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
//
// For more information, see ../documentation/library.pdf
//########################################################################################

ma = library("maths.lib");
ba = library("basics.lib");
de = library("delays.lib");
an = library("analyzers.lib");
fi = library("filters.lib"); // for compatible copy/paste out of this file

declare name "Faust Filters 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!)

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

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

//----------------------`(fi.)zero`--------------------------
// One zero filter. Difference equation: $y(n) = x(n) - zx(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>
//----------------------------------------------------------
zero(z) =  _ <: _,mem : _,*(z) : -;

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

//----------------------`(fi.)integrator`--------------------------
// Same as `pole(1)` [implemented separately for block-diagram clarity].
//------------------------------------------------------------
integrator = + ~ _ ;

//-------------------`(fi.)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) = \frac{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>
//------------------------------------------------------------
dcblockerat(fb) = *(b0) : zero(1) : pole(p)
with {
  wn = ma.PI*fb/ma.SR;
  b0 = 1.0 / (1 + wn);
  p = (1 - wn) * b0;
};

//----------------------`(fi.)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 : _
// ```
//------------------------------------------------------------
dcblocker = zero(1) : pole(0.995);

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

//------`(fi.)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>
//------------------------------------------------------------
ff_comb (maxdel,M,b0,bM) = _ <: *(b0), bM *  de.delay(maxdel,M) : + ;

//------`(fi.)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>
//------------------------------------------------------------
ff_fcomb(maxdel,M,b0,bM) = _ <: *(b0), bM * de.fdelay(maxdel,M) : + ;

//-----------`(fi.)ffcombfilter`-------------------
// Typical special case of `ff_comb()` where: `b0 = 1`.
//------------------------------------------------------------
ffcombfilter(maxdel,del,g) = ff_comb(maxdel,del,1,g);


//-----------------------`(fi.)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>
//------------------------------------------------------------
fb_comb (maxdel,N,b0,aN) = (+ <:  de.delay(maxdel,N-1),_) ~ *(-aN) : !,*(b0):mem;


//-----------------------`(fi.)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>
//------------------------------------------------------------
fb_fcomb(maxdel,N,b0,aN) = (+ <: de.fdelay(maxdel,float(N)-1.0),_) ~ *(-aN) : !,*(b0):mem;

//-----------------------`(fi.)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 `reverbs.lib` for usage examples.
//------------------------------------------------------------
rev1(maxdel,N,g) = fb_comb (maxdel,N,1,-g);

//-----`(fi.)fbcombfilter` and `(fi.)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>
//------------------------------------------------------------
fbcombfilter(maxdel,intdel,g) = (+ : de.delay(maxdel,intdel)) ~ *(g);
ffbcombfilter(maxdel,del,g) = (+ : de.fdelay(maxdel,del)) ~ *(g);


//-------------------`(fi.)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>
//------------------------------------------------------------
allpass_comb(maxdel,N,aN) = (+ <: de.delay(maxdel,N-1),*(aN)) ~ *(-aN) : mem,_ : + ;


//-------------------`(fi.)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>
//------------------------------------------------------------
allpass_fcomb(maxdel,N,aN) = (+ <: de.fdelay(maxdel,N-1),*(aN)) ~ *(-aN) : mem,_ : + ;


//-----------------------`(fi.)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 `reverbs.lib` for usage examples.
//------------------------------------------------------------
rev2(maxlen,len,g) = allpass_comb(maxlen,len,-g);

//-------------------`(fi.)allpass_fcomb5` and `(fi.)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);`).
// ```
//------------------------------------------------------------
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

//----------------------------`(fi.)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>
//------------------------------------------------------------
iir(bv,av) = ma.sub ~ fir(av) : fir(bv);

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

//---------------`(fi.)conv` and `(fi.)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,...))
// ```
//------------------------------------------------------------
//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);

//----------------`(fi.)tf1`, `(fi.)tf2` and `(fi.)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>
//------------------------------------------------------------
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;
};

//------------`(fi.)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>
//------------------------------------------------------------
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>
//========================================================================================

//----------------`(fi.)tf21`, `(fi.)tf22`, `(fi.)tf22t` and `(fi.)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>
//------------------------------------------------------------
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>
//========================================================================================

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

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

//--------------------`(fi.)iir_lat2`-----------------------
// Two-multiply latice IIR filter of arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_lat2(bv,av) : _
// ```
//
// Where:
//
// * bv: zeros as a bank of parallel signals
// * av: poles as a bank of parallel signals
//------------------------------------------------------------
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
};

//-----------------------`(fi.)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 reflection coefficients (-1 1)
//------------------------------------------------------------
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)
};

//--------------------`(fi.)iir_kl`-----------------------
// Kelly-Lochbaum ladder IIR filter of arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_kl(bv,av) : _
// ```
//
// Where:
//
// * bv: zeros as a bank of parallel signals
// * av: poles as a bank of parallel signals
//------------------------------------------------------------
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 '+'
};

//-----------------------`(fi.)allpassnklt`--------------------------
// Kelly-Lochbaum ladder allpass.
//
// #### Usage:
//
// ```
// _ : allpassklt(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflection coefficients (-1 1)
//------------------------------------------------------------
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,_);
};

//--------------------`(fi.)iir_lat1`-----------------------
// One-multiply latice IIR filter of arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_lat1(bv,av) : _
// ```
//
// Where:
//
// * bv: zeros as a bank of parallel signals
// * av: poles as a bank of parallel signals
//------------------------------------------------------------
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 '+'
};

//-----------------------`(fi.)allpassn1mt`--------------------------
// One-multiply lattice allpass with tap lines.
//
// #### Usage
//
// ```
// _ : allpassn1mt(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflection coefficients (-1 1)
//------------------------------------------------------------
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)
};

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

//-------------------------------`(fi.)allpassnnlt`-------------------------
// Normalized ladder allpass filter of arbitrary order.
//
// #### Usage:
//
// ```
// _ : allpassnnlt(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflection 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>
//------------------------------------------------------------
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=======================================
//========================================================================================

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

//-----------------------------`(fi.)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>
//------------------------------------------------------------
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 = _,_ <: !,_,_,!;
};

//-----------------------------`(fi.)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>
//------------------------------------------------------------
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 = _,_ <: !,_,_,!;
};


//------------`(fi.)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
//------------------------------------------------------------
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
//========================================================================================

//---------------`(fi.)allpassn`-----------------
// Two-multiply lattice - each section is two multiply-adds.
//
// #### Usage:
//
// ```
// _ : allpassn(n,sv) : _
// ```
// #### Where:
//
// * `n`: the order of the filter
// * `sv`: the reflection 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.
//----------------------------------------------
allpassn(0,sv) = _;
allpassn(n,sv) = _ <: ((+ <: (allpassn(n-1,sv)),*(s))~(*(-s))) : _',_ :+
with { s = ba.take(n,sv); };

//---------------`(fi.)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 reflection coefficients (-PI PI)
//----------------------------------------------
// 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)); };

//---------------`(fi.)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 reflection coefficients (-1 1)
//----------------------------------------------
// 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); };

//---------------`(fi.)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 reflection coefficients (-1 1)
//----------------------------------------------
// 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==================
//========================================================================================

//-------------------------`(fi.)tf2s` and `(fi.)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>
//----------------------------------------------
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;
};

//-----------------------------`(fi.)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) : _
// ```
//----------------------------------------------
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);
};

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

//-----------------------------`(fi.)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) : _
// ```
//----------------------------------------------
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;
};

//-----------------------------`(fi.)tf1sb`--------------------------------
// First-to-second-order lowpass-to-bandpass section mapping,
// analogous to tf2sb above.
//
// #### Usage
//
// ```
// _ : tf1sb(b1,b0,a0,w1,wc) : _
// ```
//----------------------------------------------
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==================================
//========================================================================================

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


//------------------`(fi.)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)
//---------------------------------------------------------------------
// resonhp = 2nd-order highpass with corner resonance:
resonhp(fc,Q,gain,x) = gain*x-resonlp(fc,Q,gain,x);


//------------------`(fi.)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)
//---------------------------------------------------------------------
// 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==============================
//========================================================================================

//----------------`(fi.)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');")`
//------------------------------
lowpass(N,fc) = lowpass0_highpass1(0,N,fc);


//----------------`(fi.)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');")`
//------------------------------
highpass(N,fc) = lowpass0_highpass1(1,N,fc);


//-------------`(fi.)lowpass0_highpass1`--------------
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).
//========================================================================================

//--------------------`(fi.)lowpass_plus`|`minus_highpass`----------------
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
//========================================================================================

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

//-----------------------------`(fi.)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
// ```
//----------------------------------------------------------------------
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======================================
//========================================================================================

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

//-----------------------------`(fi.)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
//-------------------------------------------------------------------------
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===========================
//========================================================================================

//--------------------`(fi.)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/>
//-------------------------------------------------------------------------
bandpass(Nh,fl,fu) = bandpass0_bandstop1(0,Nh,fl,fu);


//--------------------`(fi.)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/>
//-------------------------------------------------------------------------
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====================================
//========================================================================================

//---------------------`(fi.)bandpass6e`-----------------------------
// Order 12 elliptic bandpass filter analogous to `bandpass(6)`.
//--------------------------------------------------------------
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
};

//----------------------`(fi.)bandpass12e`---------------------------
// Order 24 elliptic bandpass filter analogous to `bandpass(6)`.
//--------------------------------------------------------------
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
//========================================================================================

//----------------------`(fi.)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`.
//--------------------------------------------------------------
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

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


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

//--------------------`(fi.)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
//------------------------------------------------------------
peak_eq_cq(Lfx,fx,Q) = peak_eq(Lfx,fx,fx/Q);

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


//---------------------`(fi.)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 > 0
// * `bw`: bandwidth of desired roll-off band
// * `alpha`: slope of roll-off desired in nepers per neper,
//         between -1 and 1 (ln mag / ln radian freq)
//
// #### Examples
// See `spectral_tilt_demo`.
//
// #### Reference
// J.O. Smith and H.F. Smith,
// "Closed Form Fractional Integration and Differentiation via Real Exponentially Spaced Pole-Zero Pairs",
// arXiv.org publication arXiv:1606.06154 [cs.CE], June 7, 2016,
// http://arxiv.org/abs/1606.06154
//------------------------------------------------------------
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
    f0p = max(f0,1.0e-7); // cannot go to zero
    w0 = 2 * ma.PI * f0p; // radian frequency of first pole
    f1 = f0p + bw; // upper band limit
    r = (f1/f0p)^(1.0/float(N-1)); // pole ratio (2 => octave spacing)
  };
};


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


//----------------------`(fi.)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>
//------------------------------------------------------------
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/
//========================================================================================

//------------------------`(fi.)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,_); // 3rd-order Butterworth
// _ : mth_octave_filterbank5(M,ftop,N) : par(i,N,_); // 5th-order Butterworth
// mth_octave_filterbank_default = mth_octave_filterbank5;
// ```
//
// 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)
//------------------------------------------------------------
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!)

//---------------`(fi.)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)) : _,_,_
// ```
//---------------------------------------------------
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);
};

//-----------------`(fi.)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)) : _,_,_
// ```
//---------------------------------------------------
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);
};

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

// TODO: Add GRAME functions here

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