// filter.lib - digital filters of various types useful in audio and beyond

declare name "Faust Filter Library";
declare author "Julius O. Smith (jos at ccrma.stanford.edu)";
declare copyright "Julius O. Smith III";
declare version "1.29";
declare license "STK-4.3"; // Synthesis Tool Kit 4.3 (MIT style license)
declare reference "https://ccrma.stanford.edu/~jos/filters/";

import("music.lib"); // delay, frac and, from math.lib, SR and PI

//---------------------- zero(z) --------------------------
// z = location of zero along real axis in z-plane
// Difference equation: y(n) = x(n) - z * x(n-1)
// Reference: https://ccrma.stanford.edu/~jos/filters/One_Zero.html

zero(z) =  _ <: _,mem : _,*(z) : -;

//------------------------ pole(p) ---------------------------
// p = pole location = feedback coefficient
// Could also be called a "leaky integrator".
// Difference equation: y(n) = x(n) + p * y(n-1)
// Reference: https://ccrma.stanford.edu/~jos/filters/One_Pole.html

pole(p) = + ~ *(p);

//---------------------- integrator --------------------------
// pole(1) [implemented separately for block-diagram clarity]

integrator = + ~ _ ;

//----------------------- tau2pole ---------------------------
// tau2pole(tau) returns a real pole giving exponential decay with
// tau = time-constant in seconds
//
tau2pole(tau) = exp(-1.0/(tau*SR));

//---------------------- smooth(s) --------------------------
// Exponential smoothing by a unity-dc-gain one-pole lowpass
//
// USAGE: smooth(tau2pole(tau)), where 
//   tau = desired smoothing time constant in seconds,
// or
//   smooth(s), where s = smoothness between 0 and 1.
//    s=0 for no smoothing
//    s=0.999 is "very smooth"
// s>1 is unstable, and s=1 yields the zero signal for all inputs.
// The exponential time-constant is approximately 
// 1/(1-s) samples, when s is close to (but less than) 1.
// Reference: 
// https://ccrma.stanford.edu/~jos/mdft/Convolution_Example_2_ADSR.html

smooth(s) = *(1.0 - s) : + ~ *(s);

//------------------- dcblockerat(fb) -----------------------
// fb = "break frequency" in Hz, i.e., -3 dB gain 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).
// Reference:
// https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html

dcblockerat(fb) = *(b0) : zero(1) : pole(p)
with {
  wn = PI*fb/SR;
  b0 = 1.0 / (1 + wn);
  p = (1 - wn) * b0;
};

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

//------------ notchw(width,freq), notch(freq) --------------
// 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;
};

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

//----------------------- ff_comb, ff_fcomb ----------------------------
// Feed-Forward Comb Filter
// 
// USAGE:
//   _ : ff_comb(maxdel,intdel,b0,bM) : _
//   _ : 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
//
// Note that ff_comb requires integer delays  (uses delay()  internally)
// while ff_fcomb takes floating-point delays (uses fdelay() internally).
//
// REFERENCE:
//   https://ccrma.stanford.edu/~jos/pasp/Feedforward_Comb_Filters.html

ff_comb (maxdel,M,b0,bM) = _ <: *(b0), bM *  delay(maxdel,M) : + ;
ff_fcomb(maxdel,M,b0,bM) = _ <: *(b0), bM * fdelay(maxdel,M) : + ;

// Typical special case:
ffcombfilter(maxdel,del,g) = ff_comb(maxdel,del,1,g);

//----------------------- fb_comb, fb_fcomb, rev1 -----------------------
// Feed-Back Comb Filter
//
// USAGE: 
//   _ : fb_comb(maxdel,intdel,b0,aN) : _
//   _ : fb_fcomb(maxdel,del,b0,aN) : _
//   _ : rev1(maxdel,del,-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) = (+ <:  delay(maxdel,N),_) ~ *(-aN) : !,*(b0);
fb_fcomb(maxdel,N,b0,aN) = (+ <: fdelay(maxdel,N),_) ~ *(-aN) : !,*(b0);

// The "rev1 section" dates back to the 1960s in computer-music reverberation.
// See the jcrev and brassrev in effect.lib for usage examples.
rev1(maxdel,N,g) = fb_comb (maxdel,N,1,-g);

// Typical special case:
fbcombfilter(maxdel,intdel,g) = (+ : delay(maxdel,intdel)) ~ *(g);
ffbcombfilter(maxdel,del,g) = (+ : fdelay(maxdel,del)) ~ *(g);

//------------------- allpass_comb, allpass_fcomb, rev2 -----------------
// Schroeder Allpass Comb Filter
//
// USAGE: 
//   _ : allpass_comb (maxdel,intdel,aN) : _
//   _ : allpass_fcomb(maxdel,del,aN) : _
//   _ : rev2(maxdel,del,-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
//
// 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.
// 
// 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) = (+ <: 
   delay(maxdel,N-1),*(aN)) ~ *(-aN)
   : mem,_ : + ;

// The "rev2 section" dates back to the 1960s in computer-music reverberation:
rev2(maxlen,len,g) = allpass_comb(maxlen,len,-g);

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

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

//---------------------------- iir (tfN) -------------------------------
// Nth-order Infinite-Impulse-Response (IIR) digital filter,
// implemented in terms of the Transfer-Function (TF) coefficients.
// Such filter structures are termed "direct form".
//
// 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) = sub ~ fir(av) : fir(bv);

//----------------------------- sub ---------------------------------
sub(x,y) = y-x; // move to math.lib?

//----------------------------- fir ---------------------------------
// FIR filter (convolution of FIR filter coefficients with a signal)
//
// USAGE: 
//   _ : fir(bv) : _
// 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);

//--------------------------- conv, 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

conv(kv,x) = sum(i,count(kv),take(i+1,kv) * x@i); // count() from math.lib

// Named special cases:
//----------------------------- tf1, tf2 ---------------------------------
// tfN = N'th-order direct-form digital filter
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)

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

// REFERENCES: 
//  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

//------------------------------ block, crossn,crossn1 ----------------------------------
// signal block/crossing utilities
// (move to math.lib?)

// block - terminate n signals (goes with bus(n) in math.lib)
block(n) = par(i,n,!);

// crossnn - cross two bus(n)s:
crossnn(n) = bus(n),bus(n) <: block(n),bus(n),bus(n),block(n);

// crossn1 - cross bus(n) and bus(1):
crossn1(n) = bus(n),(bus(1)<:bus(n)) <: block(n),bus(n),bus(n),block(n):bus(1),block(n-1),bus(n);

//------------------------------- 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 si = 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 = 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) = 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

bvav2nuv(bv,av) = par(m,M+1,nu(m)) with {
  M = count(av);
  nu(m) = 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) = take(i+1,Pr(m-1));
  Ari(m,i) = 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, allpassnt -----------------------

iir_lat2(bv,av) = allpassnt(M,sv) : sum(i,M+1,*(take(M-i+1,tg)))
with {
  M = count(av);
  sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
  tg = bvav2nuv(bv,av); // tg = vector of tap gains
};

// two-multiply lattice allpass (nested order-1 direct-form-ii allpasses):
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) = crossnn(1) : _, (_<:mem,_) : +,_;
  fsec(n) = crossn1(n) : _, (_<:mem,_),par(i,n-1,_) : +, par(i,n,_);
  innertaps(n) = par(i,n,_);
  s = take(n,sv); // reflection coefficient s = sin(theta)
};

//------------------------------- iir_kl, allpassnklt -------------------------
iir_kl(bv,av) = allpassnklt(M,sv) : sum(i,M+1,*(tghr(i)))
with {
  M = 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) = take(M+1-i,tg);
  tghr(n) = tgr(n)/pi(n);
  pi(0) = 1;
  pi(n) = pi(n-1)*(1+take(M-n+1,sv)); // all sign parameters '+'
};

// Kelly-Lochbaum ladder allpass with tap lines:
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 = take(n,sv);
  sumandtaps(n) = +,par(i,n,_);
};


//------------------------------- iir_lat1, allpassn1mt -------------------------
iir_lat1(bv,av) = allpassn1mt(M,sv) : sum(i,M+1,*(tghr(i+1)))
with {
  M = count(av);
  sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
  tg = bvav2nuv(bv,av); // tg = vector of tap gains
  tgr(i) = 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+take(M-n+2,sv)); // all sign parameters '+'
};

// one-multiply lattice allpass with tap lines:
allpassn1mt(0,sv) = _;
allpassn1mt(n,sv)= _ <: _,_ : ((+:*(s) <: _,_),_ : _,+ : crossnn(1)
		  : allpassn1mt(n-1,sv),_)~(*(-1)) : fsec(n)
with {
//0:  fsec(n) = _',_ : +
  fsec(1) = crossnn(1) : _, (_<:mem,_) : +,_;
  fsec(n) = crossn1(n) : _, (_<:mem,_),par(i,n-1,_) : +, par(i,n,_);
  innertaps(n) = par(i,n,_);
  s = take(n,sv); // reflection coefficient s = sin(theta)
};

//------------------------------- iir_nl, allpassnnlt -------------------------
// Normalized ladder filter
//
// 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 = 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) = take(M+1-i,tg);
  tghr(n) = tgr(n)/pi(n);
  pi(0) = 1;
  s(n) = take(M-n+1,sv);
  c(n) = sqrt(max(0,1-s(n)*s(n))); // compiler crashes on sqrt(-)
  pi(n) = pi(n-1)*c(n);
};

// Normalized ladder allpass with tap lines:
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 = take(n,sv);
  c = sqrt(max(0,1-s*s));
  sumandtaps(n) = +,par(i,n,_);
};

//========================= Useful special cases ============================

//-------------------------------- tf2np ------------------------------------
// tf2np - biquad based on a stable second-order Normalized Ladder Filter
// (more robust to modulation than tf2 and protected against instability)
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) = 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) = 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
//
wgr(f,r,x) = (*(G),_<:_,((+:*(C))<:_,_),_:+,_,_:+(x),-) ~ cross : _,*(0-gi)
with { 
  C = cos(2*PI*f/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
//
nlf2(f,r,x) = ((_<:_,_),(_<:_,_) : (*(s),*(c),*(c),*(0-s)) :> 
              (*(r),+(x))) ~ cross
with { 
  th = 2*PI*f/SR;
  c = cos(th);
  s = sin(th);
  cross = _,_ <: !,_,_,!;
};

//===================== 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
//   1. https://ccrma.stanford.edu/~jos/pasp/Conventional_Ladder_Filters.html
//      https://ccrma.stanford.edu/~jos/pasp/Nested_Allpass_Filters.html
//   2. Linear Prediction of Speech, Markel and Gray, Springer Verlag, 1976
//
// QUICK GUIDE
//   allpassn   - two-multiply lattice - each section is two multiply-adds
//   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).
//   allpassnkl - Kelly-Lochbaum form - four multiplies and two adds per
//                section, but all signals have an immediate physical
//                interpretation as traveling pressure waves, etc.
//   allpassn1m - One-multiply form - one multiply and three adds per section.
//                Normally the most efficient in special-purpose hardware.
//
// TYPICAL USAGE
//   _ : allpassn(N,sv) : _
// where
//   N  = allpass order (number of ladder or lattice sections)
//   sv = (s1,s2,...,sN) = reflection coefficients (between -1 and 1).
//        For allpassnn only, sv is replaced by tv, where sv(i) = sin(tv(i)),
//        where tv(i) may range between -PI and PI.
//
// two-multiply:
allpassn(0,sv) = _;
allpassn(n,sv) = _ <: ((+ <: (allpassn(n-1,sv)),*(s))~(*(-s))) : _',_ :+
with { s = take(n,sv); };

// 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(take(n,tv));  s=sin(take(n,tv)); };

// power-normalized with sparse delays dv(n)>1:
allpassnns(0,tv,dmax,dv) = _;
allpassnns(n,tv,dmax,dv) = _ <: *(s), (*(c) : (+ : dl
           : allpassnns(n-1,tv,dmax,dv))~(*(-s))) : _, mem*c : +
with { c=cos(take(n,tv));  s=sin(take(n,tv));
       dl=delay(dmax,(take(n,dv)-1)); };

// Kelly-Lochbaum:
allpassnkl(0,sv) = _;
allpassnkl(n,sv) = _ <: *(s),(*(1+s) : (+
                   : allpassnkl(n-1,sv))~(*(-s))) : _, mem*(1-s) : +
with { s = take(n,sv); };

// one-multiply:
allpassn1m(0,sv) = _;
allpassn1m(n,sv)= _ <: _,_ : ((+:*(s) <: _,_),_ : _,+ : cross
		  : allpassn1m(n-1,sv),_)~(*(-1)) : _',_ : +
with {s = take(n,sv); cross = _,_ <: !,_,_,!; };

//===== Digital Filter Sections Specified as Analog Filter Sections =====
//
//------------------------- tf2s, 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/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/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;
};

//----------------------------- 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/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.
//
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(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.
//
tf1sb(b1,b0,a0,w1,wc) = tf2(b0d/a0d,b1d/a0d,b2d/a0d,a1d/a0d,a2d/a0d) with {
  T = 1.0/float(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 = 2nd-order lowpass with corner resonance:
resonlp(fc,Q,gain) = tf2s(b2,b1,b0,a1,a0,wc)
with {
     wc = 2*PI*fc;
     a1 = 2/Q;
     a0 = 1;
     b2 = 0;
     b1 = 0;
     b0 = gain;
};

// resonhp = 2nd-order highpass with corner resonance:
resonhp(fc,Q,gain,x) = gain*x-resonlp(fc,Q,gain,x);

// resonbp = 2nd-order bandpass
resonbp(fc,Q,gain) = tf2s(b2,b1,b0,a1,a0,wc) 
with {
     wc = 2*PI*fc;
     a1 = 2/Q;
     a0 = 1;
     b2 = 0;
     b1 = gain;
     b0 = 0;
};

//================ Butterworth Lowpass/Highpass Filters ======================
// Nth-order Butterworth lowpass or highpass filters
//
// USAGE: 
//   _ : lowpass(N,fc) : _
//   _ : highpass(N,fc) : _
// where
//   N = filter order (number of poles) [nonnegative integer]
//   fc = desired cut-off frequency (-3dB frequency) in Hz
// REFERENCE: 
//  https://ccrma.stanford.edu/~jos/filters/Butterworth_Lowpass_Design.html
//  'butter' function in Octave ("[z,p,g] = butter(N,1,'s');")
// ACKNOWLEDGMENT
//  Generalized recursive formulation initiated by Yann Orlarey.

lowpass(N,fc) = lowpass0_highpass1(0,N,fc);
highpass(N,fc) = lowpass0_highpass1(1,N,fc);
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*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(-PI + (1-parity)*PI/(2*N) + (S-1+parity)*PI/N);
    w1 = 2*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 ------------------

highpass_plus_lowpass(1,fc) = _;
highpass_plus_lowpass(3,fc) = tf2s(1,-1,1,1,1,w1) with { w1 = 2*PI*fc; };
highpass_minus_lowpass(3,fc) = tf1s(-1,1,1,w1) with { w1 = 2*PI*fc; };
highpass_plus_lowpass(5,fc) = tf2s(1,-a11,1,a11,1,w1)
with { 
  a11 = 1.618033988749895;
  w1 = 2*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*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 ===================
// USAGE: 
//   _ : lowpass3e(fc) : _
//   _ : lowpass6e(fc) : _
// where fc = -3dB frequency in Hz
//
// REFERENCES:
//   http://en.wikipedia.org/wiki/Elliptic_filter
//   functions 'ncauer' and 'ellip' in Octave

//----------------------------- lowpass3e -----------------------------
// Third-order Elliptic (Cauer) lowpass filter
// 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*PI*fc; 
};

//----------------------------- lowpass6e -----------------------------
// Sixth-order Elliptic/Cauer lowpass filter
// 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*PI*fc; 
};

//===================== Elliptic Highpass Filters =====================
// USAGE: 
//   _ : highpass3e(fc) : _
//   _ : highpass6e(fc) : _
// where fc = -3dB frequency in Hz

//----------------------------- highpass3e -----------------------------
// Third-order Elliptic (Cauer) highpass filter
// DESIGN: Inversion of lowpass3e wrt unit circle in s plane (s <- 1/s)
//
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*PI*fc; 
};

//----------------------------- highpass6e -----------------------------
// Sixth-order Elliptic/Cauer highpass filter
// METHOD: Inversion of lowpass3e wrt unit circle in s plane (s <- 1/s)
//
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*PI*fc; 
};

//================== Butterworth Bandpass/Bandstop Filters =====================
// 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.
// A notch-like "bandstop" filter is similarly made from highpass(Nh).
//
// USAGE: 
//   _ : bandpass(Nh,fl,fu) : _
//   _ : bandstop(Nh,fl,fu) : _
// where 
//   Nh = HALF the desired bandpass/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/
//
bandpass(Nh,fl,fu) = bandpass0_bandstop1(0,Nh,fl,fu);
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*PI*fl; // digital (z-plane) lower passband edge
  wu = 2*PI*fu; // digital (z-plane) upper passband edge

  c = 2.0*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(-PI + (1-parity)*PI/(2*Nh) + (S-1+parity)*PI/Nh);
  };
};

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

//----------------------------- bandpass6e -----------------------------
// Order 12 elliptic bandpass filter analogous to bandpass(6) above.
//
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*PI*fl; // digital (z-plane) lower passband edge
  wu = 2*PI*fu; // digital (z-plane) upper passband edge

  c = 2.0*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 -----------------------------

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*PI*fl; // digital (z-plane) lower passband edge
  wu = 2*PI*fu; // digital (z-plane) upper passband edge

  c = 2.0*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)  ==================
// REFERENCES
// - http://en.wikipedia.org/wiki/Equalization
// - Digital Audio Signal Processing, Udo Zolzer, Wiley, 1999, p. 124
// - http://www.harmony-central.com/Computer/Programming/Audio-EQ-Cookbook.txt
//   http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
// - 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)
// USAGE: lowshelf(L0,fx), where
//  L0 = desired boost (dB) between dc and fx
//  fx = desired transition frequency (Hz) from boost to unity gain
// The gain at SR/2 is constrained to be 1.
//
low_shelf = low_shelf3; // default
low_shelf1(L0,fx,x) = x + (db2linear(L0)-1)*lowpass(1,fx,x);
low_shelf1_l(G0,fx,x) = x + (G0-1)*lowpass(1,fx,x);
low_shelf3(L0,fx,x) = x + (db2linear(L0)-1)*lowpass(3,fx,x);
low_shelf5(L0,fx,x) = x + (db2linear(L0)-1)*lowpass(5,fx,x);

//----------------------------- high_shelf -----------------------------------
// First-order "high shelf" filter (gain boost|cut above some frequency)
//
// USAGE: high_shelf(Lpi,fx), where
//  Lpi = desired boost or cut (dB) between fx and SR/2
//  fx = desired transition frequency in Hz
// The gain at dc is constrained to be 1
//
high_shelf=high_shelf3; //default
high_shelf1(Lpi,fx,x) = x + (db2linear(Lpi)-1)*highpass(1,fx,x);
high_shelf1_l(Gpi,fx,x) = x + (Gpi-1)*highpass(1,fx,x);
high_shelf3(Lpi,fx,x) = x + (db2linear(Lpi)-1)*highpass(3,fx,x);
high_shelf5(Lpi,fx,x) = x + (db2linear(Lpi)-1)*highpass(5,fx,x);

//-------------------------------- peak_eq -----------------------------------
// Second order "peaking equalizer" section 
// (gain boost or cut near some frequency)
//   Also called a "parametric equalizer" section
// USAGE: _ : peak_eq(Lfx,fx,B) : _;
// where
//    Lfx = level (dB) at fx
//     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/SR);
  Bw = B*T/sin(wx*T); // prewarp s-bandwidth for more accuracy in z-plane
  a1 = PI*Bw;
  b1 = g*a1;
  g = 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*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
//
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.)
//
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*PI*fx/SR);
  k2 = (1.0 - tanPiBT)/(1.0 + tanPiBT);
  K = db2linear(Lfx);
};

//-------------------------- parametric_eq_demo ------------------------------
// USAGE: _ : parametric_eq_demo: _ ;
parametric_eq_demo = // input_signal : 
	low_shelf(LL,FL) :
        peak_eq(LP,FP,BP) :
	high_shelf(LH,FH)
// Recommended: 
//       : mth_octave_spectral_level_demo(2) // half-octave spectrum analyzer
with {
  eq_group(x) = hgroup("[0] PARAMETRIC EQ SECTIONS 
                [tooltip: See Faust's filter.lib for info and pointers]",x);

  ls_group(x) = eq_group(vgroup("[1] Low Shelf",x));
  LL = ls_group(hslider("[0] Low Boost|Cut [unit:dB] [style:knob]
                [tooltip: Amount of low-frequency boost or cut in decibels]",
                0,-40,40,0.1));
  FL = ls_group(hslider("[1] Transition Frequency [unit:Hz] [style:knob]
                [tooltip: Transition-frequency from boost (cut) to unity gain]",
                200,1,5000,1));

  pq_group(x) = eq_group(vgroup("[2] Peaking Equalizer
                [tooltip: Parametric Equalizer sections from filter.lib]",x));
  LP = pq_group(hslider("[0] Peak Boost|Cut [unit:dB] [style:knob]
                [tooltip: Amount of local boost or cut in decibels]",
                0,-40,40,0.1));
  FP = pq_group(hslider("[1] Peak Frequency [unit:PK] [style:knob]
       [tooltip: Peak Frequency in Piano Key (PK) units (A-440= 49 PK)]",
                49,1,100,1)) : smooth(0.999) : pianokey2hz 
		with { pianokey2hz(x) = 440.0*pow(2.0, (x-49.0)/12); };

  Q = pq_group(hslider("[2] Peak Q [style:knob]
              [tooltip: Quality factor (Q) of the peak = center-frequency/bandwidth]",
              40,1,50,0.1));

  BP = FP/Q;

  hs_group(x) = eq_group(vgroup("[3] High Shelf 
                [tooltip: A high shelf provides a boost or cut 
                          above some frequency]",x));
  LH = hs_group(hslider("[0] High Boost|Cut [unit:dB] [style:knob]
                [tooltip: Amount of high-frequency boost or cut in decibels]",
                0,-40,40,.1));
  FH = hs_group(hslider("[1] Transition Frequency [unit:Hz] [style:knob]
                [tooltip: Transition-frequency from boost (cut) to unity gain]",
                8000,20,10000,1));
};

//========================= Lagrange Interpolation ========================
// Reference:
// https://ccrma.stanford.edu/~jos/pasp/Lagrange_Interpolation.html
// 
//------------------ fdelay1, fdelay2, fdelay3, fdelay4 ---------------
// Delay lines interpolated using Lagrange interpolation
// USAGE: _ : fdelayN(maxdelay, delay, inputsignal) : _
//        (exactly like fdelay in music.lib)
// where N=1,2,3, or 4 is the order of the Lagrange interpolation polynomial.
//
// NOTE: requested delay should not be less than (N-1)/2.
//
// NOTE: While the implementations below appear to use multiple delay lines,
//  they in fact use only one thanks to optimization by the Faust compiler.

// first-order case (linear interpolation) - equivalent to fdelay in music.lib
// delay d in [0,1] 
fdelay1(n,d,x)  = delay(n,id,x)*(1 - fd) + delay(n,id+1,x)*fd
with {
  id = int(d);
  fd = frac(d);
};

// second-order (quadratic) case, delay in [0.5,1.5]
// delay d should be at least 0.5
fdelay2(n,d,x) = delay(n,id,x)*(1-fd)*(2-fd)/2 
               + delay(n,id+1,x)*(2-fd)*fd
               + delay(n,id+2,x)*(fd-1)*fd/2
with {
  o = 0.49999; // offset to make life easy for interpolator
  dmo = d - o; // assumed nonnegative
  id = int(dmo);
  fd = o + frac(dmo);
};

// third-order (cubic) case, delay in [1,2]
// delay d should be at least 1
fdelay3(n,d,x) = delay(n,id,x) * (0-fdm1*fdm2*fdm3)/6
               + delay(n,id+1,x) * fd*fdm2*fdm3/2
               + delay(n,id+2,x) * (0-fd*fdm1*fdm3)/2
               + delay(n,id+3,x) * fd*fdm1*fdm2/6
with {
  id = int(d-1); 
  fd = 1+frac(d);
  fdm1 = fd-1;
  fdm2 = fd-2;
  fdm3 = fd-3;
};

// fourth-order (quartic) case, delay in [1.5,2.5]
// delay d should be at least 1.5
fdelay4(n,d,x) = delay(n,id,x)   * fdm1*fdm2*fdm3*fdm4/24 
               + delay(n,id+1,x) * (0-fd*fdm2*fdm3*fdm4)/6
               + delay(n,id+2,x) * fd*fdm1*fdm3*fdm4/4
               + delay(n,id+3,x) * (0-fd*fdm1*fdm2*fdm4)/6
               + delay(n,id+4,x) * fd*fdm1*fdm2*fdm3/24
with {
//v1:  o = 1;
  o = 1.49999;
  dmo = d - o; // assumed nonnegative
  id = int(dmo);
  fd = o + frac(dmo);
  fdm1 = fd-1;
  fdm2 = fd-2;
  fdm3 = fd-3;
  fdm4 = fd-4;
};

// fifth-order case, delay in [2,3]
// delay d should be at least 2
fdelay5(n,d,x) = 
   delay(n,id,x)   *    -fdm1*fdm2*fdm3*fdm4*fdm5/120
 + delay(n,id+1,x) *  fd*     fdm2*fdm3*fdm4*fdm5/24
 + delay(n,id+2,x) * -fd*fdm1*     fdm3*fdm4*fdm5/12
 + delay(n,id+3,x) *  fd*fdm1*fdm2*     fdm4*fdm5/12
 + delay(n,id+4,x) * -fd*fdm1*fdm2*fdm3*     fdm5/24
 + delay(n,id+5,x) *  fd*fdm1*fdm2*fdm3*fdm4     /120
with {
//v1:  o = 1;
  o = 1.99999;
  dmo = d - o; // assumed nonnegative
  id = int(dmo);
  fd = o + frac(dmo);
  fdm1 = fd-1;
  fdm2 = fd-2;
  fdm3 = fd-3;
  fdm4 = fd-4;
  fdm5 = fd-5;
};

//====================== Thiran Allpass Interpolation =====================
// Reference:
// https://ccrma.stanford.edu/~jos/pasp/Thiran_Allpass_Interpolators.html
// 
//---------------- fdelay1a, fdelay2a, fdelay3a, fdelay4a -------------
// Delay lines interpolated using Thiran allpass interpolation
// USAGE: fdelayNa(maxdelay, delay, inputsignal) 
//        (exactly like fdelay in music.lib)
// where N=1,2,3, or 4 is the order of the Thiran interpolation filter,
// and the delay argument is at least N - 1/2.
//
// (Move the following and similar notes above to filter-lib-doc.txt?)
//
// NOTE: The interpolated delay should not be less than N - 1/2.
//            (The allpass delay ranges from N - 1/2 to N + 1/2.)
//            This constraint can be alleviated by altering the code,
//            but be aware that allpass filters approach zero delay
//            by means of pole-zero cancellations.
//            The delay range [N-1/2,N+1/2] is not optimal. What is?
//
// NOTE: Delay arguments too small will produce an UNSTABLE allpass!
//
// NOTE: Because allpass interpolation is recursive, it is not as robust
//       as Lagrange interpolation under time-varying conditions.
//       (You may hear clicks when changing the delay rapidly.)
//
// first-order allpass interpolation, delay d in [0.5,1.5]
fdelay1a(n,d,x)  = delay(n,id,x) : tf1(eta,1,eta)
with {
  o = 0.49999; // offset to make life easy for allpass
  dmo = d - o; // assumed nonnegative
  id = int(dmo);
  fd = o + frac(dmo);
  eta = (1-fd)/(1+fd); // allpass coefficient
};

// second-order allpass delay in [1.5,2.5]
fdelay2a(n,d,x) = delay(n,id,x) : tf2(a2,a1,1,a1,a2)
with {
  o = 1.49999;
  dmo = d - o; // delay range is [order-1/2, order+1/2]
  id = int(dmo);
  fd = o + frac(dmo);
  a1o2 = (2-fd)/(1+fd); // share some terms (the compiler does this anyway)
  a1 = 2*a1o2;
  a2 = a1o2*(1-fd)/(2+fd);
};

// third-order allpass delay in [2.5,3.5]
// delay d should be at least 2.5
fdelay3a(n,d,x) = delay(n,id,x) : iir((a3,a2,a1,1),(a1,a2,a3))
with {
  o = 2.49999;
  dmo = d - o;
  id = int(dmo);
  fd = o + frac(dmo);
  a1o3 = (3-fd)/(1+fd);
  a2o3 = a1o3*(2-fd)/(2+fd);
  a1 = 3*a1o3;
  a2 = 3*a2o3;
  a3 = a2o3*(1-fd)/(3+fd);
};

// fourth-order allpass delay in [3.5,4.5]
// delay d should be at least 3.5
fdelay4a(n,d,x) = delay(n,id,x) : tf4(a4,a3,a2,a1,1,a1,a2,a3,a4)
with {
  o = 3.49999;
  dmo = d - o;
  id = int(dmo);
  fd = o + frac(dmo);
  a1o4 = (4-fd)/(1+fd);
  a2o6 = a1o4*(3-fd)/(2+fd);
  a3o4 = a2o6*(2-fd)/(3+fd);
  a1 = 4*a1o4;
  a2 = 6*a2o6;
  a3 = 4*a3o4;
  a4 = a3o4*(1-fd)/(4+fd);
};

//================ Mth-Octave Filter-Banks and Spectrum-Analyzers ============
// Mth-octave filter-banks and spectrum-analyzers split the input signal into a
// bank of parallel signals, one for each spectral band. 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.
//
// A SPECTRUM-ANALYZER is defined here as any band-split whose bands span
// the relevant spectrum, but whose band-signals do not
// necessarily sum to the original signal, either exactly or to within an
// allpass filtering. Spectrum analyzer outputs are normally at least nearly
// "power complementary", i.e., the power spectra of the individual bands
// sum to the original power spectrum (to within some negligible tolerance).
//
// 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_analyzer ----------------------------
//
// USAGE
//  _ : mth_octave_analyzer(O,M,ftop,N) : par(i,N,_); // Oth-order Butterworth
//  _ : mth_octave_analyzer6e(M,ftop,N) : par(i,N,_); // 6th-order elliptic
//
// 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)
//
// ACKNOWLEDGMENT
//  Recursive band-splitting formulation improved by Yann Orlarey.

mth_octave_analyzer6e(M,ftop,N) = _ <: bsplit(N-1) with {
  fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
  lp(n) = lowpass6e(fc(n));  // 6th-order elliptic - see other choices above
  hp(n) = highpass6e(fc(n)); //   (search for lowpass* and highpass*)
  bsplit(0)  = _; 
  bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
};

// Butterworth analyzers may be cascaded with allpass
// delay-equalizers to make (allpass-complementary) filter banks:

mth_octave_analyzer(O,M,ftop,N) = _ <: bsplit(N-1) with {
  fc(n) = ftop * 2^(float(n-N+1)/float(M));
  lp(n) = lowpass(O,fc(n)); // Order O Butterworth
  hp(n) = highpass(O,fc(n));
  bsplit(0)  = _; 
  bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
};

mth_octave_analyzer3(M,ftop,N) = mth_octave_analyzer(3,M,ftop,N);
mth_octave_analyzer5(M,ftop,N) = mth_octave_analyzer(5,M,ftop,N);
mth_octave_analyzer_default = mth_octave_analyzer6e; // default analyzer

//------------------------ mth_octave_filterbank -------------------------
// Allpass-complementary filter banks based on Butterworth band-splitting.
// For Butterworth band-splits, the needed delay equalizer is easily found.

mth_octave_filterbank(O,M,ftop,N) = 
    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) = 
    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;

//======================= Mth-Octave Spectral Level =========================
// Spectral Level: Display (in bar graphs) the average signal level in each 
//                 spectral band.
//
//------------------------ mth_octave_spectral_level -------------------------
// USAGE: _ : mth_octave_spectral_level(M,ftop,NBands,tau,dB_offset);
// where 
//    M = bands per octave
//    ftop = lower edge frequency of top band
//    NBands = number of passbands (including highpass and dc bands),
//    tau = spectral display averaging-time (time constant) in seconds,
//    dB_offset = constant dB offset in all band level meters.
//
mth_octave_spectral_level6e(M,ftop,N,tau,dB_offset) = _<:
    _,mth_octave_analyzer6e(M,ftop,N) : 
    _,(display:>_):attach with {
  display = par(i,N,dbmeter(i));
  dbmeter(i) = abs : smooth(tau2pole(tau)) : linear2db : +(dB_offset) : 
     meter(N-i-1);
    meter(i) = speclevel_group(vbargraph("[%2i] [unit:dB] 
     [tooltip: Spectral Band Level in dB]", -50, 10));
  // Can M be included in the label string somehow?
  speclevel_group(x)  = hgroup("[0] CONSTANT-Q SPECTRUM ANALYZER (6E)
     [tooltip: See Faust's filter.lib for documentation and references]", x);
};

mth_octave_spectral_level_default = mth_octave_spectral_level6e;
spectral_level = mth_octave_spectral_level(2,10000,20);  // simplest case

//---------------------- mth_octave_spectral_level_demo ----------------------
// Demonstrate mth_octave_spectral_level in a standalone GUI.
//
// USAGE:  _ : mth_octave_spectral_level_demo(BandsPerOctave);

mth_octave_spectral_level_demo(M) = 
  mth_octave_spectral_level_default(M,ftop,N,tau,dB_offset) 
with {
  // Span nearly 10 octaves so that lowest band-edge is at 
  // ftop*2^(-Noct+2) = 40 Hz when ftop=10 kHz:
  N = int(10*M); // without 'int()', segmentation fault observed for M=1.67
  ftop = 10000;
  ctl_group(x)  = hgroup("[1] SPECTRUM ANALYZER CONTROLS", x);
  tau = ctl_group(hslider("[0] Level Averaging Time [unit:sec]
        [tooltip: band-level averaging time in seconds]",
        0.1,0,1,0.01)); 
  dB_offset = ctl_group(hslider("[1] Level dB Offset [unit:dB]
        [tooltip: Level offset in decibels]",
        50,0,100,1)); 
};

spectral_level_demo = mth_octave_spectral_level_demo(1.5); // 2/3 octave

//---------------- (third|half)_octave_(analyzer|filterbank) -----------------

// Named special cases of mth_octave_* with defaults filled in:

third_octave_analyzer(N) = mth_octave_analyzer_default(3,10000,N);
third_octave_filterbank(N) = mth_octave_filterbank_default(3,10000,N);
// Third-Octave Filter-Banks have been used in audio for over a century.
// See, e.g.,
//   Acoustics [the book], by L. L. Beranek
//   Amer. Inst. Physics for the Acoustical Soc. America,
//   http://asa.aip.org/publications.html, 1986 (1st ed.1954)

// Third-octave bands across the audio spectrum are too wide for current
// typical computer screens, so half-octave bands are the default:
half_octave_analyzer(N) = mth_octave_analyzer_default(2,10000,N);
half_octave_filterbank(N) = mth_octave_filterbank_default(2,10000,N);

octave_filterbank(N) = mth_octave_filterbank_default(1,10000,N);
octave_analyzer(N) = mth_octave_analyzer_default(1,10000,N);

//=========================== Filter-Bank Demos ==============================
// Graphic Equalizer: Each filter-bank output signal routes through a fader.
//
// USAGE: _ : mth_octave_filterbank_demo(M) : _
// where 
//   M = number of bands per octave

mth_octave_filterbank_demo(M) = bp1(bp,mthoctavefilterbankdemo) with {
   bp1 = component("effect.lib").bypass1;
   mofb_group(x)  = vgroup("CONSTANT-Q FILTER BANK (Butterworth dyadic tree)
     [tooltip: See Faust's filter.lib for documentation and references]", x);
   bypass_group(x) = mofb_group(hgroup("[0]", x));
   slider_group(x)  = mofb_group(hgroup("[1]", x));
   N = 10*M; // total number of bands (highpass band, octave-bands, dc band)
   ftop = 10000;
   mthoctavefilterbankdemo = chan;
   chan = mth_octave_filterbank_default(M,ftop,N) : 
          sum(i,N,(*(db2linear(fader(N-i)))));
   fader(i) = slider_group(vslider("[%2i] [unit:dB]
              [tooltip: Bandpass filter gain in dB]", -10, -70, 10, 0.1)) : 
              smooth(0.999);
   bp = bypass_group(checkbox("[0] Bypass  
              [tooltip: When this is checked, the filter-bank has no effect]"));
};

filterbank_demo = mth_octave_filterbank_demo(1); // octave-bands = default

//=========== Arbritary-Crossover Filter-Banks and Spectrum Analyzers ========
// These are similar to the Mth-octave filter-banks above, except that the
// band-split frequencies are passed explicitly as arguments. 
//
// USAGE:
//    _ : filterbank (O,freqs) : par(i,N,_); // Butterworth band-splits
//    _ : filterbanki(O,freqs) : par(i,N,_); // Inverted-dc version 
//    _ : analyzer   (O,freqs) : par(i,N,_); // No delay equalizer
//
// 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)) : _,_,_
//
// ACKNOWLEDGMENT
//    Technique for processing a variable number of signal arguments due
//    to Yann Orlarey (as is the entire Faust framework!)
//
//------------------------------ analyzer --------------------------------------
analyzer(O,lfreqs) = _ <: bsplit(nb) with
{
   nb = count(lfreqs);
   fc(n) = take(n, lfreqs);
   lp(n) = lowpass(O,fc(n));
   hp(n) = highpass(O,fc(n));
   bsplit(0) = _;
   bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
};

//----------------------------- filterbank -------------------------------------
filterbank(O,lfreqs) = analyzer(O,lfreqs) : delayeq with
{
   nb = count(lfreqs);
   fc(n) = take(n, lfreqs);
   ap(n) = highpass_plus_lowpass(O,fc(n));
   delayeq = par(i,nb-1,apchain(nb-1-i)),_,_;
   apchain(0) = _;
   apchain(i) =  ap(i) : apchain(i-1);
};

//----------------------------- filterbanki ------------------------------------
filterbanki(O,lfreqs) = _ <: bsplit(nb) with
{
   fc(n) = 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);
};