filters.lib

//##################################### filters.lib ########################################
// Filters library. Its official prefix is `fi`.
//
// The Filters library is organized into 23 sections:
//
// * [Basic Filters](#basic-filters)
// * [Comb Filters](#comb-filters)
// * [Direct-Form Digital Filter Sections](#direct-form-digital-filter-sections)
// * [Direct-Form Second-Order Biquad Sections](#direct-form-second-order-biquad-sections)
// * [Ladder/Lattice Digital Filters](#ladderlattice-digital-filters)
// * [Useful Special Cases](#useful-special-cases)
// * [Ladder/Lattice Allpass Filters](#ladderlattice-allpass-filters)
// * [Digital Filter Sections Specified as Analog Filter Sections](#digital-filter-sections-specified-as-analog-filter-sections)
// * [Simple Resonator Filters](#simple-resonator-filters)
// * [Butterworth Lowpass/Highpass Filters](#butterworth-lowpasshighpass-filters)
// * [Special Filter-Bank Delay-Equalizing Allpass Filters](#special-filter-bank-delay-equalizing-allpass-filters)
// * [Elliptic (Cauer) Lowpass Filters](#elliptic-cauer-lowpass-filters)
// * [Elliptic Highpass Filters](#elliptic-highpass-filters)
// * [Butterworth Bandpass/Bandstop Filters](#butterworth-bandpassbandstop-filters)
// * [Elliptic Bandpass Filters](#elliptic-bandpass-filters)
// * [Parametric Equalizers (Shelf, Peaking)](#parametric-equalizers-shelf-peaking)
// * [Mth-Octave Filter-Banks](#mth-octave-filter-banks)
// * [Arbitrary-Crossover Filter-Banks and Spectrum Analyzers](#arbitrary-crossover-filter-banks-and-spectrum-analyzers)
// * [State Variable Filters (SVF)](#state-variable-filters)
// * [Linkwitz-Riley 4th-order 2-way, 3-way, and 4-way crossovers](#linkwitz-riley-4th-order-2-way-3-way-and-4-way-crossovers)
// * [Standardized Filters](#standardized-filters)
// * [Averaging Functions](#averaging-functions)
// * [Kalman Filters](#kalman-filters)
//
// #### References
// * <https://github.com/grame-cncm/faustlibraries/blob/master/filters.lib>
//
//########################################################################################

// NOTE ABOUT LICENSES:
// Each function in this library has its own license. Licenses are declared
// before each function. Corresponding license terms can be found at the
// bottom of this file or in the Faust libraries documentation.

ma = library("maths.lib");
ba = library("basics.lib");
ro = library("routes.lib");
de = library("delays.lib");
an = library("analyzers.lib");
ef = library("misceffects.lib");
si = library("signals.lib");
fi = library("filters.lib"); // for compatible copy/paste out of this file
la = library("linearalgebra.lib");

declare name "Faust Filters Library";
declare version "1.7.1";

//===============================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>
//----------------------------------------------------------
declare zero author "Julius O. Smith III";
declare zero copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare zero license "MIT-style STK-4.3 license";
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>
//------------------------------------------------------------
declare pole author "Julius O. Smith III";
declare pole copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare pole license "MIT-style STK-4.3 license";
pole(p) = + ~ *(p);

//----------------------`(fi.)integrator`--------------------------
// Same as `pole(1)` [implemented separately for block-diagram clarity].
//------------------------------------------------------------
declare integrator author "Julius O. Smith III";
declare integrator copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare integrator license "MIT-style STK-4.3 license";
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 f_b)}$$
// (which can be seen as a 1st-order Butterworth highpass filter)
// by the low-frequency-matching bilinear transform method
// (i.e., using the typical frequency-scaling constant `2*SR`).
//
// #### Usage
//
// ```
// _ : dcblockerat(fb) : _
// ```
//
// Where:
//
// * `fb`: "break frequency" in Hz, i.e., -3 dB gain frequency (see 2nd reference below)
//
// #### References
// * <https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html>
// * <https://ccrma.stanford.edu/~jos/spectilt/Bode_Plots.html>
//------------------------------------------------------------
declare dcblockerat author "Julius O. Smith III";
declare dcblockerat copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare dcblockerat license "MIT-style STK-4.3 license";
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 : _
// ```
//------------------------------------------------------------
declare dcblocker author "Julius O. Smith III";
declare dcblocker copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare dcblocker license "MIT-style STK-4.3 license";
dcblocker = zero(1) : pole(0.995);

//----------------------------`(fi.)lptN`--------------------------------------
// One-pole lowpass filter with arbitrary dis/charging factors set in dB and
// times set in seconds.
//
// #### Usage
//
// ```
// _ : lptN(N, tN) : _
// ```
//
// Where:
//
// * `N`: is the attenuation factor in dB
// * `tN`: is the filter period in seconds, that is, the time for the
// impulse response to decay by `N` dB
//
// #### Reference
// <https://ccrma.stanford.edu/~jos/mdft/Exponentials.html>
//----------------------------------------------------------
declare lptN author "Julius O. Smith III";
declare lptN copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare lptN license "MIT-style STK-4.3 license";
lptN(N, tN, x) = x : si.smooth(ba.tau2pole(tN / log(10.0^(float(N)/20.0))));
// Special cases of lptN
lptau(tN, x) = lptN(8.6858896381, tN, x); // Tau time constant, i.e., 1/e atten. after tN secs
lpt60(tN, x) = lptN(60, tN, x); // T60 constant, i.e., 1/1000 atten. after tN secs
lpt19(tN, x) = lptN(19, tN, x); // T19 constant, i.e., 1/e^2.2 atten. after tN secs

//=======================================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>
//------------------------------------------------------------
declare ff_comb author "Julius O. Smith III";
declare ff_comb copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare ff_comb license "MIT-style STK-4.3 license";
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)
// * `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>
//------------------------------------------------------------
declare ff_fcomb author "Julius O. Smith III";
declare ff_fcomb copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare ff_fcomb license "MIT-style STK-4.3 license";
ff_fcomb(maxdel,del,b0,bM) = _ <: *(b0), bM * de.fdelay(maxdel,del) : +;

//-----------`(fi.)ffcombfilter`-------------------
// Typical special case of `ff_comb()` where: `b0 = 1`.
//------------------------------------------------------------
declare ff_combfilter author "Julius O. Smith III";
declare ff_combfilter copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare ff_combfilter license "MIT-style STK-4.3 license";
ffcombfilter(maxdel,del,g) = ff_comb(maxdel,del,1,g);

//---------------------`(fi.)fb_comb_common`---------------------
// A generic feedback comb filter.
//
// #### Usage
//
// ```
// _ : fb_comb_common(dop,N,b0,aN) : _
// ```
//
// Where
//
// * `dop`: delay operator, e.g. `@` or `de.fdelay4a(2048)`
// * `N`: current delay
// * `b0`: gain applied to input
// * `aN`: gain applied to delay-line output
//
// #### Example test program
//
// ```
// process = fb_comb_common(@,N,b0,aN);
// ```
// implements the following difference equation:
// ```
// y[n] = b0 x[n] + aN y[n - N]
// ```
//
// See more examples in `filters.lib` below.
// --------------------------------------------------------
declare fb_comb_common author "Oleg Nesterov";
fb_comb_common(dop,N,b0,aN) = + ~ aN * dop(N-1) : *(b0);

//-----------------------`(fi.)fb_comb`-----------------------
// Feed-Back Comb Filter (integer delay).
//
// #### Usage
//
// ```
// _ : fb_comb(maxdel,del,b0,aN) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `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>
//------------------------------------------------------------
declare fb_comb author "Julius O. Smith III";
declare fb_comb copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>, revised by Oleg Nesterov";
declare fb_comb license "MIT-style STK-4.3 license";
fb_comb(maxdel,del,b0,aN) = fb_comb_common(de.delay(maxdel),del,b0,-aN) : 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)
// * `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>
//------------------------------------------------------------
declare fb_fcomb author "Julius O. Smith III";
declare fb_fcomb copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>, revised by Oleg Nesterov";
declare fb_fcomb license "MIT-style STK-4.3 license";
fb_fcomb(maxdel,del,b0,aN) = fb_comb_common(de.fdelay(maxdel),del,b0,-aN) : 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.
//------------------------------------------------------------
declare rev1 author "Julius O. Smith III";
declare rev1 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare rev1 license "MIT-style STK-4.3 license";
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>
//------------------------------------------------------------
declare fbcombfilter author "Julius O. Smith III";
declare fbcombfilter copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare fbcombfilter license "MIT-style STK-4.3 license";
fbcombfilter(maxdel,intdel,g) = (+ : de.delay(maxdel,intdel)) ~ *(g);

declare ffbcombfilter author "Julius O. Smith III";
declare ffbcombfilter copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare ffbcombfilter license "MIT-style STK-4.3 license";
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>
//------------------------------------------------------------
declare allpass_comb author "Julius O. Smith III";
declare allpass_comb copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpass_comb license "MIT-style STK-4.3 license";
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>
//------------------------------------------------------------
declare allpass_fcomb author "Julius O. Smith III";
declare allpass_fcomb copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpass_fcomb license "MIT-style STK-4.3 license";
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.
//------------------------------------------------------------
declare rev2 author "Julius O. Smith III";
declare rev2 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare rev2 license "MIT-style STK-4.3 license";
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);`).
// ```
//------------------------------------------------------------
declare allpass_fcomb5 author "Julius O. Smith III";
declare allpass_fcomb5 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpass_fcomb5 license "MIT-style STK-4.3 license";
allpass_fcomb5(maxdel,N,aN) = (+ <: de.fdelay5(maxdel,N-1),*(aN)) ~ *(-aN) : mem,_ : +;

declare allpass_fcomb1a author "Julius O. Smith III";
declare allpass_fcomb1a copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpass_fcomb1a license "MIT-style STK-4.3 license";
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:
//
// * `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>
//------------------------------------------------------------
declare iir author "Julius O. Smith III";
declare iir copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare iir license "MIT-style STK-4.3 license";
iir(bv,av) = ma.sub ~ fir(av) : fir(bv);

//-----------------------------`(fi.)fir`---------------------------------
// FIR filter (convolution of FIR filter coefficients with a signal). `fir` is standard Faust function.
//
// #### 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 test program
//
// 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);
declare fir author "Julius O. Smith III";
declare fir copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare fir license "MIT-style STK-4.3 license";
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

declare convN author "Julius O. Smith III";
declare convN copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare convN license "MIT-style STK-4.3 license";
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

declare conv author "Julius O. Smith III";
declare conv copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare conv license "MIT-style STK-4.3 license";
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:
//
// * `b`: transfer-function numerator
// * `a`: transfer-function denominator (monic)
//
// #### Reference
// <https://ccrma.stanford.edu/~jos/fp/Direct_Form_I.html>
//------------------------------------------------------------
declare tf1 author "Julius O. Smith III";
declare tf1 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf1 license "MIT-style STK-4.3 license";
tf1(b0,b1,a1) = _ <: *(b0), (mem : *(b1)) :> + ~ *(0-a1);

declare tf2 author "Julius O. Smith III";
declare tf2 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf2 license "MIT-style STK-4.3 license";
tf2(b0,b1,b2,a1,a2) = iir((b0,b1,b2),(a1,a2));
// tf2 is a variant of tf22 below with duplicated mems

declare tf3 author "Julius O. Smith III";
declare tf3 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf3 license "MIT-style STK-4.3 license";
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>
//------------------------------------------------------------
declare notchw author "Julius O. Smith III";
declare notchw copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare notchw license "MIT-style STK-4.3 license";
notchw(width,freq) = tf2(b0,b1,b2,a1,a2)
with {
  fb = 0.5*width; // First design a dcblockerat(width/2)
  wn = ma.PI*fb/ma.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*ma.PI*freq/ma.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:
//
// * `b`: transfer-function numerator
// * `a`: transfer-function denominator (monic)
//
// #### Reference
// <https://ccrma.stanford.edu/~jos/fp/Direct_Form_I.html>
//------------------------------------------------------------
declare tf21 author "Julius O. Smith III";
declare tf21 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf21 license "MIT-style STK-4.3 license";
tf21(b0,b1,b2,a1,a2) = // tf2, direct-form 1:
    _ <:(mem<:((mem:*(b2)),*(b1))),*(b0) :>_
    : ((_,_,_:>_) ~(_<:*(-a1),(mem:*(-a2))));

declare tf22 author "Julius O. Smith III";
declare tf22 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf22 license "MIT-style STK-4.3 license";
tf22(b0,b1,b2,a1,a2) = // tf2, direct-form 2:
    _ : (((_,_,_:>_)~*(-a1)<:mem,*(b0))~*(-a2))
      : (_<:mem,*(b1)),_ : *(b2),_,_ :> _;

declare tf22t author "Julius O. Smith III";
declare tf22t copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf22t license "MIT-style STK-4.3 license";
tf22t(b0,b1,b2,a1,a2) = // tf2, direct-form 2 transposed:
    _ : (_,_,(_ <: *(b2)',*(b1)',*(b0))
      : _,+',_,_ :> _)~*(-a1)~*(-a2) : _;

declare tf21t author "Julius O. Smith III";
declare tf21t copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf21t license "MIT-style STK-4.3 license";
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).
//------------------------------------------------------------
declare av2sv author "Julius O. Smith III";
declare av2sv copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare av2sv license "MIT-style STK-4.3 license";
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
//------------------------------------------------------------
declare bvav2nuv author "Julius O. Smith III";
declare bvav2nuv copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare bvav2nuv license "MIT-style STK-4.3 license";
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 lattice IIR filter of arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_lat2(bv,av) : _
// ```
//
// Where:
//
// * `bv`: transfer-function numerator
// * `av`: transfer-function denominator (monic)
//------------------------------------------------------------
declare iir_lat2 author "Julius O. Smith III";
declare iir_lat2 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare iir_lat2 license "MIT-style STK-4.3 license";
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), with taps.
//
// #### Usage
//
// ```
// _ : allpassnt(n,sv) : si.bus(n+1)
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflection coefficients (-1 1)
//
// The first output is the n-th order allpass output,
// while the remaining outputs are taps taken from the
// input of each delay element from the input to the output.
// See (fi.)allpassn for the single-output case.
//------------------------------------------------------------
declare allpassnt author "Julius O. Smith III";
declare allpassnt copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpassnt license "MIT-style STK-4.3 license";
allpassnt(0,sv) = _;
allpassnt(n,sv) = _ : ((+ <: (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`: transfer-function numerator
// * `av`: transfer-function denominator (monic)
//------------------------------------------------------------
declare iir_kl author "Julius O. Smith III";
declare iir_kl copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare iir_kl license "MIT-style STK-4.3 license";
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:
//
// ```
// _ : allpassnklt(n,sv) : _
// ```
//
// Where:
//
// * `n`: the order of the filter
// * `sv`: the reflection coefficients (-1 1)
//------------------------------------------------------------
declare allpassnklt author "Julius O. Smith III";
declare allpassnklt copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpassnklt license "MIT-style STK-4.3 license";
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 lattice IIR filter of arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_lat1(bv,av) : _
// ```
//
// Where:
//
// * bv: transfer-function numerator as a bank of parallel signals
// * av: transfer-function denominator as a bank of parallel signals
//------------------------------------------------------------
declare iir_lat1 author "Julius O. Smith III";
declare iir_lat1 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare iir_lat1 license "MIT-style STK-4.3 license";
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 (fixed at compile time)
// * `sv`: the reflection coefficients (-1 1)
//------------------------------------------------------------
declare allpassn1mt author "Julius O. Smith III";
declare allpassn1mt copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpassn1mt license "MIT-style STK-4.3 license";
allpassn1mt(0,sv) = _;
allpassn1mt(n,sv) = _ <: _,_ : ((+:*(s) <: _,_),_ : _,+ : ro.crossnn(1)
          : allpassn1mt(n-1,sv),_)~(*(-1)) : 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_nl`-------------------------
// Normalized ladder filter of arbitrary order.
//
// #### Usage
//
// ```
// _ : iir_nl(bv,av) : _
// ```
//
// Where:
//
// * `bv`: transfer-function numerator
// * `av`: transfer-function denominator (monic)
//
// #### 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>
//------------------------------------------------------------
declare iir_nl author "Julius O. Smith III";
declare iir_nl copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare iir_nl license "MIT-style STK-4.3 license";
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 (fixed at compile time)
// * `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>
//------------------------------------------------------------
declare allpassnnlt author "Julius O. Smith III";
declare allpassnnlt copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpassnnlt license "MIT-style STK-4.3 license";
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:
//
// * `b`: transfer-function numerator
// * `a`: transfer-function denominator (monic)
//------------------------------------------------------------
declare tf2np author "Julius O. Smith III";
declare tf2np copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf2np license "MIT-style STK-4.3 license";
tf2np(b0,b1,b2,a1,a2) = allpassnnlt(M,sv) : sum(i,M+1,*(tghr(i)))
with {
  smax = 1.0-ma.EPSILON; // 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>
//------------------------------------------------------------
declare wgr author "Julius O. Smith III";
declare wgr copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare wgr license "MIT-style STK-4.3 license";
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>
//------------------------------------------------------------
declare nlf2 author "Julius O. Smith III";
declare nlf2 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare nlf2 license "MIT-style STK-4.3 license";
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
//------------------------------------------------------------
declare apnl author "Julius O. Smith III";
declare apnl copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare apnl license "MIT-style STK-4.3 license";
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.)scatN`--------------------------
// N-port scattering junction.
//
// #### Usage
//
// ```
// si.bus(N) : scatN(N,av,filter) : si.bus(N)
// ```
//
// Where:
//
// * `N`: number of incoming/outgoing waves
// * `av`: vector (list) of `N` alpha parameters (each between 0 and 2, and normally summing to 2): <https://ccrma.stanford.edu/~jos/pasp/Alpha_Parameters.html>
// * `filter` : optional junction filter to apply (`_` for none, see below)
//
// With no filter:
//
// - The junction is _lossless_ when the alpha parameters sum to 2 ("allpass").
// - The junction is _passive_ but lossy when the alpha parameters sum to less than 2 ("resistive loss").
// - Dynamic and reactive junctions are obtained using the `filter` argument.
//   For guaranteed stability, the filter should be _positive real_. (See 2nd ref. below).
//
// For \(N=2\) (two-port scattering), the reflection coefficient \(\rho\) corresponds
// to alpha parameters \(1\pm\rho\).
//
// #### Example: Whacky echo chamber made of 16 lossless "acoustic tubes":
//
// ```
// process = _ : *(1.0/sqrt(N)) <: daisyRev(16,2,0.9999) :> _,_ with { 
//   daisyRev(N,Dp2,G) = si.bus(N) : (si.bus(2*N) :> si.bus(N)
//     : fi.scatN(N, par(i,N,2*G/float(N)), fi.lowpass(1,5000.0))
//     : par(i,N,de.delay(DS(i),DS(i)-1))) ~ si.bus(N) with { DS(i) = 2^(Dp2+i); };
// };
// ```
//
// #### References
// * <https://ccrma.stanford.edu/~jos/pasp/Loaded_Waveguide_Junctions.html>
// * <https://ccrma.stanford.edu/~jos/pasp/Passive_String_Terminations.html>
// * <https://ccrma.stanford.edu/~jos/pasp/Unloaded_Junctions_Alpha_Parameters.html>
//------------------------------------------------------------
declare scatN author "Julius O. Smith III";
declare scatN copyright "Copyright (C) 2024 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare scatN license "MIT-style STK-4.3 license";

scatN(0,av,filter) = !;
scatN(N,av,filter) = incomingWaves <: (junctionSum : filter <: si.bus(N)), par(i,N,*(-1)) :> si.bus(N)
with {
  incomingWaves = si.bus(N);
  alpha(i) = ba.take(i+1,av);
  junctionSum = par(i,N,*(alpha(i))) :> _; // Junction velocity/pressure for series/parallel junction
  outgoingWaves = junctionSum <: si.bus(N), negatedIncoming :> si.bus(N);
  alphaSum = sum(i,N,alpha(i));
};

//---------------`(fi.)scat`-----------------
// Scatter off of reflectance r with reflection coefficient s.
//
// #### Usage:
//
// ```
// _ : scat(s,r) : _
// ```
// #### Where:
//
// * `s`: reflection coefficient between -1 and 1 for stability
// * `r`: single-input, single-output block diagram,
//        having gain less than 1 at all frequencies for stability.
//
// #### Example:  The following program should produce all zeros:
//
// ```
// process = fi.allpassn(3,(.3,.2,.1)), fi.scat(.1, fi.scat(.2, fi.scat(.3, _)))
//           :> - : ^(2) : +~_;
// ```
//
// #### Reference:
// * <https://ccrma.stanford.edu/~jos/pasp/Scattering_Impedance_Changes.html>
//----------------------------------------------
declare scat author "Julius O. Smith III";
declare scat copyright "Copyright (C) 2024 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare scat license "MIT-style STK-4.3 license";

scat(s,r) = _ <: ((+ <: (r,*(s)))~(*(-s))) : _',_ :+;

//---------------`(fi.)allpassn`-----------------
// Two-multiply lattice filter.
//
// #### Usage:
//
// ```
// _ : allpassn(n,sv) : _
// ```
// #### Where:
//
// * `n`: the order of the filter
// * `sv`: the reflection coefficients (-1 1)
// * `sv`: the reflection coefficients  (s1,s2,...,sN), each between -1 and 1.
//
// Equivalent to `fi.allpassnt(n,sv) : _, par(i,n,!);`
// Equivalent to `fi.scat( s(n), fi.scat( s(n-1), ..., fi.scat( s(1), _ )))
//               with { s(k) = ba.take(k,sv); } ;`
// Identical to `allpassn` in `old/filter.lib`.
//
// #### References
// * 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>
//----------------------------------------------
declare allpassn author "Julius O. Smith III";
declare allpassn copyright "Copyright (c) 2003-2024 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpassn license "MIT-style STK-4.3 license";
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)):
declare allpassnn author "Julius O. Smith III";
declare allpassnn copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpassnn license "MIT-style STK-4.3 license";
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.)allpassnkl`-----------------
// 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:
declare allpassnnkl author "Julius O. Smith III";
declare allpassnnkl copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpassnnkl license "MIT-style STK-4.3 license";
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:
declare allpassn1m author "Julius O. Smith III";
declare allpassn1m copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpassn1m license "MIT-style STK-4.3 license";
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 test program
//
// 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>
//----------------------------------------------
declare tf2s author "Julius O. Smith III";
declare tf2s copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf2s license "MIT-style STK-4.3 license";
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.)tf1snp`-------------------------------
// First-order special case of tf2snp above.
//
// #### Usage
//
// ```
// _ : tf1snp(b1,b0,a0) : _
// ```
//----------------------------------------------
declare tf1snp author "Julius O. Smith III";
declare tf1snp copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf1snp license "MIT-style STK-4.3 license";
tf1snp(b1,b0,a0,w1) = fi.tf2snp(b1,b0,0,a0,0,w1); // FIXME: Faust compiler does not fully optimize - does C++?

//-----------------------------`(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) : _
// ```
//----------------------------------------------
declare tf3slf author "Julius O. Smith III";
declare tf3slf copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf3slf license "MIT-style STK-4.3 license";
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 test program
//
// 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>
//----------------------------------------------
declare tf1s author "Julius O. Smith III";
declare tf1s copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf1s license "MIT-style STK-4.3 license";
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 prototype
// 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) : _
// ```
//----------------------------------------------
declare tf2sb author "Julius O. Smith III";
declare tf2sb copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf2sb license "MIT-style STK-4.3 license";
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) : _
// ```
//----------------------------------------------
declare tf1sb author "Julius O. Smith III";
declare tf1sb copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare tf1sb license "MIT-style STK-4.3 license";
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:
declare resonlp author "Julius O. Smith III";
declare resonlp copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare resonlp license "MIT-style STK-4.3 license";
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:
declare resonhp author "Julius O. Smith III";
declare resonhp copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare resonhp license "MIT-style STK-4.3 license";
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
declare resonbp author "Julius O. Smith III";
declare resonbp copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare resonbp license "MIT-style STK-4.3 license";
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 numerical expression
// * `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');")`
//------------------------------
declare lowpass author "Julius O. Smith III";
declare lowpass copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare lowpass license "MIT-style STK-4.3 license";
lowpass(N,fc) = lowpass0_highpass1(0,N,fc);


//----------------`(fi.)highpass`--------------------
// Nth-order Butterworth highpass filter.
// `highpass` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : highpass(N,fc) : _
// ```
//
// Where:
//
// * `N`: filter order (number of poles), nonnegative constant numerical expression
// * `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');")`
//------------------------------
declare highpass author "Julius O. Smith III";
declare highpass copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare higpass license "MIT-style STK-4.3 license";
highpass(N,fc) = lowpass0_highpass1(1,N,fc);


//-------------`(fi.)lowpass0_highpass1`--------------
declare lowpass0_highpass1 author "Julius O. Smith III";
declare lowpass0_highpass1 "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare lowpass0_highpass1 "MIT-style STK-4.3 license";
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`----------------
declare highpass_plus_lowpass author "Julius O. Smith III";
declare highpass_plus_lowpass copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highpass_plus_lowpass license "MIT-style STK-4.3 license";
highpass_plus_lowpass(1,fc) = _;
highpass_plus_lowpass(3,fc) = tf2s(1,-1,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;
};

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

declare highpass_minus_lowpass author "Julius O. Smith III";
declare highpass_minus_lowpass copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highpass_minus_lowpass license "MIT-style STK-4.3 license";
highpass_minus_lowpass(3,fc) = tf1s(-1,1,1,w1) with { 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_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);
};

declare highpass_plus_lowpass_even author "Julius O. Smith III";
declare highpass_plus_lowpass_even copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highpass_plus_lowpass_even license "MIT-style STK-4.3 license";
highpass_plus_lowpass_even(N,fc) = highpass(N,fc) + lowpass(N,fc);

declare highpass_minus_lowpass_even author "Julius O. Smith III";
declare highpass_minus_lowpass_even copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highpass_plus_lowpass_even license "MIT-style STK-4.3 license";
highpass_minus_lowpass_even(N,fc) = highpass(N,fc) - lowpass(N,fc);

declare highpass_plus_lowpass_odd author "Julius O. Smith III";
declare highpass_plus_lowpass_odd copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highpass_plus_lowpass_odd license "MIT-style STK-4.3 license";
// 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);

declare highpass_minus_lowpass_odd author "Julius O. Smith III";
declare highpass_minus_lowpass_odd copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highpass_plus_lowpass_odd license "MIT-style STK-4.3 license";
// FIXME: Rewrite the following, as for orders 3 and 5 above,
//        to eliminate pole-zero cancellations:
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
// ```
//---------------------------------------------------------------------
declare lowpass3e author "Julius O. Smith III";
declare lowpass3e copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare lowpass3e license "MIT-style STK-4.3 license";
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
// ```
//----------------------------------------------------------------------
declare lowpass6e author "Julius O. Smith III";
declare lowpass6e copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare lowpass6e license "MIT-style STK-4.3 license";
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
//-------------------------------------------------------------------------
declare highpass3e author "Julius O. Smith III";
declare highpass3e copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highpass3e license "MIT-style STK-4.3 license";
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
//-------------------------------------------------------------------------
declare highpass6e author "Julius O. Smith III";
declare highpass6e copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highpass6e license "MIT-style STK-4.3 license";
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`.
//
//-------------------------------------------------------------------------
declare bandpass author "Julius O. Smith III";
declare bandpass copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare bandpass license "MIT-style STK-4.3 license";
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`.
//
//-------------------------------------------------------------------------
declare bandstop author "Julius O. Smith III";
declare bandstop copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare bandstop license "MIT-style STK-4.3 license";
bandstop(Nh,fl,fu) = bandpass0_bandstop1(1,Nh,fl,fu);

declare bandpass0_bandstop1 author "Julius O. Smith III";
declare bandpass0_bandstop1 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare bandpass0_bandstop1 license "MIT-style STK-4.3 license";
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-plane) lower cutoff
  wua = c*tan(wu/c); // analog (s-plane) 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)`.
//--------------------------------------------------------------
declare bandpass6e author "Julius O. Smith III";
declare bandpass6e copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare bandpass6e license "MIT-style STK-4.3 license";
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-plane) lower cutoff
  wua = c*tan(wu/c); // analog (s-plane) 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)`.
//--------------------------------------------------------------
declare bandpass12e author "Julius O. Smith III";
declare bandpass12e copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare bandpass12e license "MIT-style STK-4.3 license";
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-plane) lower cutoff
  wua = c*tan(wu/c); // analog (s-plane) upper cutoff

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

//------------------------`(fi.)pospass`---------------------------
// Positive-Pass Filter (single-side-band filter).
//
// #### Usage
//
// ```
// _ : pospass(N,fc) : _,_
// ```
//
// where
//
// * `N`: filter order (Butterworth bandpass for positive frequencies).
// * `fc`: lower bandpass cutoff frequency in Hz.
//   - Highpass cutoff frequency at ma.SR/2 - fc Hz.
//
// #### Example test program
//
// * See `dm.pospass_demo`
// * Look at frequency response
//
// #### Method
//
// A filter passing only positive frequencies can be made from a
// half-band lowpass by modulating it up to the positive-frequency range.
// Equivalently, down-modulate the input signal using a complex sinusoid at -SR/4 Hz,
// lowpass it with a half-band filter, and modulate back up by SR/4 Hz.
// In Faust/math notation:
// $$pospass(N) = \ast(e^{-j\frac{\pi}{2}n}) : \mbox{lowpass(N,SR/4)} : \ast(e^{j\frac{\pi}{2}n})$$
//
// An approximation to the Hilbert transform is given by the
// imaginary output signal:
//
// ```
// hilbert(N) = pospass(N) : !,*(2);
// ```
//
// #### References
// * <https://ccrma.stanford.edu/~jos/mdft/Analytic_Signals_Hilbert_Transform.html>
// * <https://ccrma.stanford.edu/~jos/sasp/Comparison_Optimal_Chebyshev_FIR_I.html>
// * <https://ccrma.stanford.edu/~jos/sasp/Hilbert_Transform.html>
//------------------------------------------------------------
declare pospass author "Julius O. Smith III";
declare pospass copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare pospass license "MIT-style STK-4.3 license";
pospass(N,fc) = fi.pospass0(lpf) with {
  fcs = ma.SR/4 - fc; // Upper lowpass cutoff = (SR/2 - fc) - SR/4
  lpf = fi.lowpass(N,fcs); // Butterworth lowpass
};

declare pospass6e author "Julius O. Smith III";
declare pospass6e copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare pospass6e license "MIT-style STK-4.3 license";
pospass6e(fc) = fi.pospass0(lpf) with {
  lpf = fi.lowpass6e(ma.SR/4 - fc); // Elliptic lowpass, order 6
};

declare pospass0 author "Julius O. Smith III";
declare pospass0 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare pospass0 license "MIT-style STK-4.3 license";
pospass0(lpf) = unmodulate : lpf, lpf : modulate with {
  c = 1-1' : +~(*(-1):mem); // complex sinusoid rotating at SR/4
  s = c'; // ||: 0, 1, 0, -1 :||
  unmodulate = _ <: *(c),*(-s); // subtract SR/4 from all input frequencies
  modulate(x,y) = c*x-s*y, c*y + s*x; // add SR/4 to all frequencies
};


//=================Parametric Equalizers (Shelf, Peaking)=================================
// Parametric Equalizers (Shelf, Peaking).
//
// #### References
// * <http://en.wikipedia.org/wiki/Equalization>
// * <https://webaudio.github.io/Audio-EQ-Cookbook/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, a constant numerical expression)
// * `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`.
//
// #### References
// See "Parametric Equalizers" above for references regarding
// `low_shelf`, `high_shelf`, and `peak_eq`.
//
//--------------------------------------------------------------
declare lowshelf author "Julius O. Smith III";
declare lowshelf copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare lowshelf license "MIT-style STK-4.3 license";
lowshelf(N,L0,fx) = filterbank(N,(fx)) : _, *(ba.db2linear(L0)) :> _;

// Special cases and optimization:
declare low_shelf author "Julius O. Smith III";
declare low_shelf copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare low_shelf license "MIT-style STK-4.3 license";
low_shelf  = lowshelf(3); // default = 3rd order Butterworth

declare low_shelf1 author "Julius O. Smith III";
declare low_shelf1 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare low_shelf1 license "MIT-style STK-4.3 license";
low_shelf1(L0,fx,x) = x + (ba.db2linear(L0)-1)*lowpass(1,fx,x); // optimized

declare low_shelf1_l author "Julius O. Smith III";
declare low_shelf1_l copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare low_shelf1_l license "MIT-style STK-4.3 license";
low_shelf1_l(G0,fx,x) = x + (G0-1)*lowpass(1,fx,x); // optimized

declare lowshelf_other_freq author "Julius O. Smith III";
declare lowshelf_other_freq copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare lowshelf_other_freq license "MIT-style STK-4.3 license";
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, a constant numerical expression).
// * `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.
//
// #### References
// See "Parametric Equalizers" above for references regarding
// `low_shelf`, `high_shelf`, and `peak_eq`.
//
//--------------------------------------------------------------
declare highshelf author "Julius O. Smith III";
declare highshelf copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highshelf license "MIT-style STK-4.3 license";
highshelf(N,Lpi,fx) = filterbank(N,(fx)) : *(ba.db2linear(Lpi)), _ :> _;
// Special cases and optimization:
high_shelf = highshelf(3); // default = 3rd order Butterworth

declare high_shelf1 author "Julius O. Smith III";
declare high_shelf1 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare high_shelf1 license "MIT-style STK-4.3 license";
high_shelf1(Lpi,fx,x) = x + (ba.db2linear(Lpi)-1)*highpass(1,fx,x); // optimized

declare high_shelf1_l author "Julius O. Smith III";
declare high_shelf1_l copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare high_shelf1_l license "MIT-style STK-4.3 license";
high_shelf1_l(Gpi,fx,x) = x + (Gpi-1)*highpass(1,fx,x); //optimized

// shelf transitions between frequency fx and this one:
declare highshelf_other_freq author "Julius O. Smith III";
declare highshelf_other_freq copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare highshelf_other_freq license "MIT-style STK-4.3 license";
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
//
// #### References
// See "Parametric Equalizers" above for references regarding
// `low_shelf`, `high_shelf`, and `peak_eq`.
//
//--------------------------------------------------------------
declare peak_eq author "Julius O. Smith III";
declare peak_eq copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare peak_eq license "MIT-style STK-4.3 license";
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
//
// #### References
// See "Parametric Equalizers" above for references regarding
// `low_shelf`, `high_shelf`, and `peak_eq`.
//
//------------------------------------------------------------
declare peak_eq_cq author "Julius O. Smith III";
declare peak_eq_cq copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare peak_eq_cq license "MIT-style STK-4.3 license";
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.)
// See also "Parametric Equalizers" above for references on shelf
// and peaking equalizers in general.
//
//------------------------------------------------------------
declare peak_eq_rm author "Julius O. Smith III";
declare peak_eq_rm copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare peak_eq_rm license "MIT-style STK-4.3 license";
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)
//
// #### Example test program
//
// See `dm.spectral_tilt_demo` and the documentation for `no.pink_noise`.
//
// #### 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>
//
//------------------------------------------------------------
declare spectral_tilt author "Julius O. Smith III";
declare spectral_tilt copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare spectral_tilt license "MIT-style STK-4.3 license";
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,ma.EPSILON); // 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>
//------------------------------------------------------------
declare levelfilter author "Julius O. Smith III";
declare levelfilter copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare levelfilter license "MIT-style STK-4.3 license";
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:
//
// * `N`: Number of filters in series where L = L/N, a constant numerical expression
// * `freq`: corner frequency (-3dB point) usually set to fundamental freq
// * `L`: desired level (in dB) at Nyquist limit (SR/2), e.g., -60
//
// #### Reference
// <https://ccrma.stanford.edu/realsimple/faust_strings/Dynamic_Level_Lowpass_Filter.html>
//------------------------------------------------------------
declare levelfilterN author "Julius O. Smith III";
declare levelfilterN copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare levelfilterN license "MIT-style STK-4.3 license";
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), a constant numerical expression
// * `N`: total number of bands (>2), a constant numerical expression
// * `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, a constant numerical expression
// * `M`: number of band-slices per octave, a constant numerical expression
// * `ftop`: highest band-split crossover frequency (e.g., 20 kHz)
// * `N`: total number of bands (including dc and Nyquist), a constant numerical expression
//------------------------------------------------------------
declare mth_octave_filterbank author "Julius O. Smith III";
declare mth_octave_filterbank copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare mth_octave_filterbank license "MIT-style STK-4.3 license";
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).
declare mth_octave_filterbank_alt author "Julius O. Smith III";
declare mth_octave_filterbank_alt copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare mth_octave_filterbank_alt license "MIT-style STK-4.3 license";
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.
declare mth_octave_filterbank3 author "Julius O. Smith III";
declare mth_octave_filterbank3 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare mth_octave_filterbank3 license "MIT-style STK-4.3 license";
mth_octave_filterbank3(M,ftop,N) = mth_octave_filterbank_alt(3,M,ftop,N);

declare mth_octave_filterbank5 author "Julius O. Smith III";
declare mth_octave_filterbank5 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare mth_octave_filterbank5 license "MIT-style STK-4.3 license";
mth_octave_filterbank5(M,ftop,N) = mth_octave_filterbank(5,M,ftop,N);

declare mth_octave_filterbank_default author "Julius O. Smith III";
declare mth_octave_filterbank_default copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare mth_octave_filterbank_default license "MIT-style STK-4.3 license";
mth_octave_filterbank_default = mth_octave_filterbank5;


//===============Arbitrary-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], a constant numerical expression)
// * `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)) : _,_,_
// ```
//---------------------------------------------------
declare filterbank author "Julius O. Smith III";
declare filterbank copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare filterbank license "MIT-style STK-4.3 license";
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]`, a constant numerical expression)
// * `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)) : _,_,_
// ```
//---------------------------------------------------
declare filterbanki author "Julius O. Smith III";
declare filterbanki copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare filterbanki license "MIT-style STK-4.3 license";
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);
};


//===============State Variable Filters=========================================================
// #### References
// Solving the continuous SVF equations using trapezoidal integration
//
// <https://cytomic.com/files/dsp/SvfLinearTrapOptimised2.pdf>
//========================================================================================

//-----------------`(fi.)svf`----------------------
// An environment with `lp`, `bp`, `hp`, `notch`, `peak`, `ap`, `bell`, `ls`, `hs` SVF based filters.
// All filters have `freq` and `Q` parameters, the `bell`, `ls`, `hs` ones also have a `gain` third parameter.
//
// #### Usage
//
// ```
// _ : svf.xx(freq, Q, [gain]) : _
// ```
//
// Where:
//
// * `freq`: cut frequency
// * `Q`: quality factor
// * `[gain]`: gain in dB
//
/// ```
//---------------------------------------------------
declare svf author "Oleg Nesterov";
declare svf copyright "Copyright (C) 2020 Oleg Nesterov <oleg@redhat.com>";
declare svf license "MIT-style STK-4.3 license";

svf = environment {

    // Internal implementation
    svf(T,F,Q,G) = tick ~ (_,_) : !,!,si.dot(3, mix)
    with {
        tick(ic1eq, ic2eq, v0) =
            2*v1 - ic1eq,
            2*v2 - ic2eq,
            v0, v1, v2
        with {
            v1 = ic1eq + g *(v0-ic2eq) : /(1 + g*(g+k));
            v2 = ic2eq + g * v1;
        };

        A = pow(10.0, G/40.0);

        g = tan(F * ma.PI/ma.SR) : case {
            (7) => /(sqrt(A));
            (8) => *(sqrt(A));
            (t) => _;
        } (T);

        k = case {
            (6) => 1/(Q*A);
            (t) => 1/Q;
        } (T);

        mix = case {
            (0) => 0, 0, 1;
            (1) => 0, 1, 0;
            (2) => 1, -k, -1;
            (3) => 1, -k, 0;
            (4) => 1, -k, -2;
            (5) => 1, -2*k, 0;
            (6) => 1, k*(A*A-1), 0;
            (7) => 1, k*(A-1), A*A-1;
            (8) => A*A, k*(1-A)*A, 1-A*A;
        } (T);
    };

    // External API
    lp(f,q)     = svf(0, f, q, 0);
    bp(f,q)     = svf(1, f, q, 0);
    hp(f,q)     = svf(2, f, q, 0);
    notch(f,q)  = svf(3, f, q, 0);
    peak(f,q)   = svf(4, f, q, 0);
    ap(f,q)     = svf(5, f, q, 0);
    bell(f,q,g) = svf(6, f, q, g);
    ls(f,q,g)   = svf(7, f, q, g);
    hs(f,q,g)   = svf(8, f, q, g);
};


//-----------------`(fi.)svf_morph`--------------------
// An SVF-based filter that can smoothly morph between
// being lowpass, bandpass, and highpass.
//
// #### Usage
//
// ```
// _ : svf_morph(freq, Q, blend) : _
// ```
//
// Where:
//
// * `freq`: cutoff frequency
// * `Q`: quality factor
// * `blend`: [0..2] continuous, where 0 is `lowpass`, 1 is `bandpass`, and 2 is `highpass`. For performance, the value is not clamped to [0..2].
//
//
// #### Example test program
// 
// ```
// process = no.noise : svf_morph(freq, q, blend)
// with {
//   blend = hslider("Blend", 0, 0, 2, .01) : si.smoo;
//   q = hslider("Q", 1, 0.1, 10, .01) : si.smoo;
//   freq = hslider("freq", 5000, 100, 18000, 1) : si.smoo;
// };
// ```
//
// #### Reference
// <https://github.com/mtytel/vital/blob/636ca0ef517a4db087a6a08a6a8a5e704e21f836/src/synthesis/filters/digital_svf.cpp#L292-L295>
//-----------------------------------------------------
declare svf_morph author "David Braun and Clarence W. Rowley";
declare svf_morph copyright "Copyright (C) 2024 David Braun <braun@ccrma.stanford.edu>";
declare svf_morph license "MIT-style STK-4.3 license";
svf_morph(f, q, _b, x) = svf.lp(f,q,x)*w_LP + svf.bp(f,q,x)*w_BP + svf.hp(f,q,x)*w_HP
with {
  b = _b - 1; // b is now -1 to 1.
  w_LP = max(-b, 0);
  w_BP = sqrt(1 - b^2);
  w_HP = max(b, 0);
};


//-----------------`(fi.)svf_notch_morph`--------------------
// An SVF-based notch-filter that can smoothly morph between
// being lowpass, notch, and highpass.
//
// #### Usage
//
// ```
// _ : svf_notch_morph(freq, Q, blend) : _
// ```
//
// Where:
//
// * `freq`: cutoff frequency
// * `Q`: quality factor
// * `blend`: [0..2] continuous, where 0 is `lowpass`, 1 is `notch`, and 2 is `highpass`. For performance, the value is not clamped to [0..2].
//
//
// #### Example test program
// 
// ```
// process = no.noise : svf_notch_morph(freq, q, blend)
// with {
//   blend = hslider("Blend", 0, 0, 2, .01) : si.smoo;
//   q = hslider("Q", 1, 0.1, 10, .01) : si.smoo;
//   freq = hslider("freq", 5000, 100, 18000, 1) : si.smoo;
// };
// ```
//
// #### Reference
// <https://github.com/mtytel/vital/blob/636ca0ef517a4db087a6a08a6a8a5e704e21f836/src/synthesis/filters/digital_svf.cpp#L256C36-L263>
//-----------------------------------------------------------
declare svf_notch_morph author "David Braun and Clarence W. Rowley";
declare svf_notch_morph copyright "Copyright (C) 2024 David Braun <braun@ccrma.stanford.edu>";
declare svf_notch_morph license "MIT-style STK-4.3 license";
svf_notch_morph(f, q, _b, x) = svf.lp(f,q,x)*w_LP + svf.hp(f,q,x)*w_HP
with {
  b = _b - 1; // b is now -1 to 1.
  w_LP = min(1-b, 1);
  w_HP = min(1+b, 1);
};


//----------`(fi.)SVFTPT`---------------------------------------------------------
//
// Topology-preserving transform implementation following Zavalishin's method.
//
// Outputs: lowpass, highpass, bandpass, normalised bandpass, notch, allpass, 
// peaking.
//
// Each individual output can be recalled with its name in the environment as in:
//      `SVFTPT.LP2(1000.0, .707)`.
//
// The 7 outputs can be recalled by using `SVF` name as in:
//      `SVFTPT.SVF(1000.0, .707)`.
//
// Even though the implementation is different, the characteristics of this
// filter are comparable to those of the `svf` environment in this library.
//
// #### Usage:
//
// ```
// _ : SVFTPT.xxx(CF, Q) : _
// ```
//
// Where:
//
// * `xxx` can be one of the following: `LP2`, `HP2`, `BP2`, `BP2Norm`, `Notch2`, `AP2`, `Peaking2`
// * `CF`: cutoff in Hz
// * `Q`: resonance
// 
//------------------------------------------------------------------------------
declare SVFTPT author "Dario Sanfilippo";
declare SVFTPT copyright
    "Copyright (C) 2024 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare SVFTPT license "MIT License";
SVFTPT = environment {
    SVF(CF, Q, x) = f ~ si.bus(2) : (! , ! , _ , _ , _ , _ , _ , _ , _)
        with {
            g = tan(CF * ma.PI * ma.T);
            R2 = 1.0 / Q;
            gPlusR2 = g + R2;
            f(s0, s1) = u0 , u1 , LP , HP , BP , BPNorm , Notch , AP , Peaking
                with {
                    HP = (x - s0 * gPlusR2 - s1) / (1.0 + g * gPlusR2);
                    v0 = HP * g;
                    BP = s0 + v0;
                    v1 = BP * g;
                    LP = s1 + v1;
                    BPNorm = BP * R2;
                    Notch = x - BPNorm;
                    AP = x - BP * (R2 + R2);
                    Peaking = LP - HP;
                    u0 = v0 + BP;
                    u1 = v1 + LP;
                };
        };
    LP2(CF, Q, x) = SVF(CF, Q, x) : ba.selectn(7, 0);
    HP2(CF, Q, x) = SVF(CF, Q, x) : ba.selectn(7, 1);
    BP2(CF, Q, x) = SVF(CF, Q, x) : ba.selectn(7, 2);
    BP2Norm(CF, Q, x) = SVF(CF, Q, x) : ba.selectn(7, 3);
    Notch2(CF, Q, x) = SVF(CF, Q, x) : ba.selectn(7, 4);
    AP2(CF, Q, x) = SVF(CF, Q, x) : ba.selectn(7, 5);
    Peaking2(CF, Q, x) = SVF(CF, Q, x) : ba.selectn(7, 6);
};

//----------`(fi.)dynamicSmoother`------------------------------------------------
//
// Adaptive smoother based on Andy Simper's paper.
// 
// This filter uses both the lowpass and bandpass outputs of a 
// state-variable filter. The lowpass is used to smooth out the input signal,
// the bandpass, which is a smoothed out version of the highpass, provides
// information on the rate of change of the input. Hence, the bandpass signal
// can be used to adjust the cutoff of the filter to quickly follow the input's
// fast and large variations while effectively filtering out local 
// perturbations.
//
// This implementation does not use an approximation for the CF computation,
// and it deploys guards to prevent overshooting with extreme sensitivity 
// values.
//
// #### Usage:
//
// ```
// _ : dynamicSmoother(sensitivity, baseCF) : _
// ```
//
// Where:
//
// * `sensitivity`: sensitivity to changes in the input signal.
//      The range is, theoretically, from 0 to INF, though anything between
//      0.0 and 1.0 should be reasonable
// * `baseCF`: cutoff frequency, in Hz, when there is no variation in the 
//      input signal
//
// #### Reference
//  <https://cytomic.com/files/dsp/DynamicSmoothing.pdf>
//
//------------------------------------------------------------------------------
declare dynamicSmoothing author "Dario Sanfilippo";
declare dynamicSmoothing copyright
    "Copyright (C) 2024 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare dynamicSmoothing license "MIT License";
dynamicSmoothing(sensitivity, baseCF, x) = f ~ _ : ! , _
    with {
        f(s) = SVFTPT.BP2(CF, .5 , x) , SVFTPT.LP2(CF, .5, x)
            with {
                CF = min(ma.SR * .125, baseCF + sensitivity * abs(s) * ma.SR * .5);
            };
    };


//-----------------`(fi.)oneEuro`----------------------------------
// The One Euro Filter (1€ Filter) is an adaptive lowpass filter.
// This kind of filter is commonly used in object-tracking,
// not necessarily audio processing.
// 
// #### Usage
//
// ```
// _ : oneEuro(derivativeCutoff, beta, minCutoff) : _
// ```
//
// Where:
//
// * `derivativeCutoff`: Used to filter the first derivative of the input. 1 Hz is a good default.
// * `beta`: "Speed" parameter where higher values reduce latency.
// * `minCutoff`: Minimum cutoff frequency in Hz. Lower values remove more jitter.
//
// #### References
// * <https://gery.casiez.net/1euro/>
//------------------------------------------------------------------------
declare oneEuro author "David Braun";
declare oneEuro copyright "Copyright (C) 2024 by David Braun <braun@ccrma.stanford.edu>";
declare oneEuro license "MIT";
oneEuro(derivativeCutoff, beta, minCutoff) = _oneEuro ~ _
with {
    // exponential moving average that calculates its own alpha from a cutoff in Hz
    // _ : ema(cutoff) : _
    ema(cutoff) = tick ~ _
    with {
        alpha = 1.0 / (1.0 + ma.SR / (2 * ma.PI * cutoff));
        tick(prev, x) = prev*(1-alpha) + x*alpha;
    };

    _oneEuro(prev, x) = x : ema(adaptiveCutoff)
    with {
        derivative = (x - prev)*ma.SR;
        derivativeFiltered = derivative : ema(derivativeCutoff);
        adaptiveCutoff =  minCutoff + beta * abs(derivativeFiltered);  
    };
};


//===========Linkwitz-Riley 4th-order 2-way, 3-way, and 4-way crossovers=====
//
// The Linkwitz-Riley (LR) crossovers are designed to produce a fully-flat
// magnitude response when their outputs are combined. The 4th-order
// LR filters (LR4) have a 24dB/octave slope and they are rather popular audio
// crossovers used in multi-band processing.
//
// The LR4 can be constructed by cascading two second-order Butterworth
// filters. For the second-order Butterworth filters, we will use the SVF
// filter implemented above by setting the Q-factor to 1.0 / sqrt(2.0).
// These will be cascaded in pairs to build the LR4 highpass and lowpass.
// For the phase correction, we will use the 2nd-order Butterworth allpass.
//
// #### Reference
// Zavalishin, Vadim. "The art of VA filter design." Native Instruments, Berlin, Germany (2012).
//=============================================================================


//----------`(fi.)lowpassLR4`---------------------------------------------------
// 4th-order Linkwitz-Riley lowpass.
//
// #### Usage
//
// ```
// _ : lowpassLR4(cf) : _
// ```
//
// Where:
//
// * `cf` is the lowpass cutoff in Hz
//------------------------------------------------------------------------------
declare lowpassLR4 author "Dario Sanfilippo";
declare lowpassLR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare lowpassLR4 license "MIT-style STK-4.3 license";
lowpassLR4(cf, x) = x : seq(i, 2, svf.lp(cf, 1.0 / sqrt(2.0)));


//----------`(fi.)highpassLR4`--------------------------------------------------
// 4th-order Linkwitz-Riley highpass.
//
// #### Usage
//
// ```
// _ : highpassLR4(cf) : _
// ```
//
// Where:
//
// * `cf` is the highpass cutoff in Hz
//------------------------------------------------------------------------------
declare highpassLR4 author "Dario Sanfilippo";
declare highpassLR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare highpassLR4 license "MIT-style STK-4.3 license";
highpassLR4(cf, x) = x : seq(i, 2, svf.hp(cf, 1.0 / sqrt(2.0)));


//----------`(fi.)crossover2LR4`------------------------------------------------
// Two-way 4th-order Linkwitz-Riley crossover.
//
// #### Usage
//
// ```
// _ : crossover2LR4(cf) : si.bus(2)
// ```
//
// Where:
//
// * `cf` is the crossover split cutoff in Hz
//------------------------------------------------------------------------------
declare crossover2LR4 author "Dario Sanfilippo";
declare crossover2LR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare crossover2LR4 license "MIT-style STK-4.3 license";
crossover2LR4(cf, x) = lowpassLR4(cf, x) , highpassLR4(cf, x);


//----------`(fi.)crossover3LR4`------------------------------------------------
// Three-way 4th-order Linkwitz-Riley crossover.
//
// #### Usage
//
// ```
// _ : crossover3LR4(cf1, cf2) : si.bus(3)
// ```
//
// Where:
//
// * `cf1` is the crossover lower split cutoff in Hz
// * `cf2` is the crossover upper split cutoff in Hz
//------------------------------------------------------------------------------
declare crossover3LR4 author "Dario Sanfilippo";
declare crossover3LR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare crossover3LR4 license "MIT-style STK-4.3 license";
crossover3LR4(cf1, cf2, x) =
    crossover2LR4(cf1, x) : svf.ap(cf2, 1.0 / sqrt(2.0)) , crossover2LR4(cf2);


//----------`(fi.)crossover4LR4`------------------------------------------------
// Four-way 4th-order Linkwitz-Riley crossover.
//
// #### Usage
//
// ```
// _ : crossover4LR4(cf1, cf2, cf3) : si.bus(4)
// ```
//
// Where:
//
// * `cf1` is the crossover lower split cutoff in Hz
// * `cf2` is the crossover mid split cutoff in Hz
// * `cf3` is the crossover upper split cutoff in Hz
//------------------------------------------------------------------------------
declare crossover4LR4 author "Dario Sanfilippo";
declare crossover4LR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare crossover4LR4 license "MIT-style STK-4.3 license";
crossover4LR4(cf1, cf2, cf3, x) =
    crossover2LR4(cf2, x) :
        svf.ap(cf3, 1.0 / sqrt(2.0)) ,
        svf.ap(cf1, 1.0 / sqrt(2.0)) :
            crossover2LR4(cf1) ,
            crossover2LR4(cf3);


//----------`(fi.)crossover8LR4`------------------------------------------------
// Eight-way 4th-order Linkwitz-Riley crossover.
//
// #### Usage
//
// ```
// _ : crossover8LR4(cf1, cf2, cf3, cf4, cf5, cf6, cf7) : si.bus(8)
// ```
//
// Where:
//
// * `cf1-cf7` are the crossover cutoff frequencies in Hz
//------------------------------------------------------------------------------
declare crossover8LR4 author "Dario Sanfilippo";
declare crossover8LR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare crossover8LR4 license "MIT-style STK-4.3 license";
crossover8LR4(cf1, cf2, cf3, cf4, cf5, cf6, cf7, x) =
    crossover2LR4(cf4, x) :
        (ap(cf6) : ap(cf5) : ap(cf7)) ,
        (ap(cf2) : ap(cf1) : ap(cf3)) :
            crossover2LR4(cf2) ,
            crossover2LR4(cf6) :
                ap(cf3) ,
                ap(cf1) ,
                ap(cf7) ,
                ap(cf5) :
                    crossover2LR4(cf1) ,
                    crossover2LR4(cf3) ,
                    crossover2LR4(cf5) ,
                    crossover2LR4(cf7)
    with {
        ap(cf) = svf.ap(cf, 1.0 / sqrt(2.0));
    };


//=========================== Standardized Filters ============================
//=============================================================================
//
// This section provides filters that are defined by national or
// international standards, e.g. for measurement applications.

//----------------------`(fi.)itu_r_bs_1770_4_kfilter`-------------------------
// The prefilter from Recommendation ITU-R BS.1770-4 for loudness
// measurement. Also known as "K-filter". The recommendation defines
// biquad filter coefficients for a fixed sample rate of 48kHz (page
// 4-5). Here, we construct biquads for arbitrary samplerates.  The
// resulting filter is normalized, such that the magnitude at 997Hz is
// unity gain 1.0.
//
// Please note, the ITU-recommendation handles the normalization in
// equation (2) by subtracting 0.691dB, which is not needed with
// `itu_r_bs_1770_4_kfilter`.
//
// One option for future improvement might be, to round those filter
// coefficients, that are almost equal to one. Second, the maximum
// magnitude difference at 48kHz between the ITU-defined filter and
// `itu_r_bs_1770_4_kfilter` is 0.001dB, which obviously could be
// less.
//
// #### Usage
//
// ```
// _ : itu_r_bs_1770_4_kfilter : _
// ```
//
// #### Reference
// <https://www.itu.int/rec/R-REC-BS.1770>
// <https://gist.github.com/jkbd/07521a98f7873a2dc3dbe16417930791>
//-----------------------------------------------------------------------------
declare itu_r_bs_1770_4_kfilter author "Jakob Dübel";
declare itu_r_bs_1770_4_kfilter copyright "Copyright (C) 2022 Jakob Dübel";
declare itu_r_bs_1770_4_kfilter license "ISC license";

itu_r_bs_1770_4_kfilter = stage1 : stage2 : normalize997Hz
with {
  freq2k(f_c) = tan((ma.PI * f_c)/ma.SR);

  stage1 = tf22t(b0,b1,b2,a1,a2)
  with {
    f_c = 1681.7632251028442; // Hertz
    gain = 3.9997778685513232; // Decibel
    K = freq2k(f_c);
    V_0 = pow(10, (gain/20.0));

    denominator = 1.0 + sqrt(2.0)*K + K^2;
    b0 = (V_0 + sqrt((2.0*V_0))*K + K^2) / denominator;
    b1 = 2.0*(K^2 - V_0) / denominator;
    b2 = (V_0 - sqrt(2.0*V_0)*K + K^2) / denominator;

    a1 = 2*(K^2 - 1) / denominator;
    a2 = (1 - sqrt(2.0)*K + K^2) / denominator;
  };

  stage2 = tf22t(b0,b1,b2,a1,a2)
  with {
    f_c = 38.135470876002174; // Hertz
    Q = 0.5003270373223665;
    K = freq2k(f_c);

    denominator = (K^2) * Q + K + Q;
    b0 = Q / denominator;
    b1 = -2*Q / denominator;
    b2 = b0;

    a1 = (2*Q * (K^2 - 1)) / denominator;
    a2 = ((K^2) * Q - K + Q) / denominator;
  };

  normalize997Hz = *(0.9273671710547968);
};


//============================Averaging Functions==============================
//=============================================================================
//
// These are a set of samplerate independent averaging functions based on
// moving-average and one-pole filters with specific response characteristics.

//----------------------------`(fi.)avg_rect`----------------------------------
// Moving average.
//
// #### Usage
//
// ```
// _ : avg_rect(period) : _
// ```
//
// Where:
//
// * `period` is the averaging frame in seconds
//-----------------------------------------------------------------------------
declare avg_rect author "Dario Sanfilippo and Julius O. Smith III";
declare avg_rect copyright "Copyright (C) 2020 Dario Sanfilippo
      <sanfilippo.dario@gmail.com> and
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare avg_rect license "MIT-style STK-4.3 license";
avg_rect(period, x) = x : ba.slidingMean(rint(period * ma.SR));


//----------------------------`(fi.)avg_tau`-------------------------------------
// Averaging function based on a one-pole filter and the tau response time.
// Tau represents the effective length of the one-pole impulse response,
// that is, tau is the integral of the filter's impulse response. This
// response is slower to reach the final value but has less ripples in
// non-steady signals.
//
// #### Usage
//
// ```
// _ : avg_tau(period) : _
// ```
//
// Where:
//
// * `period` is the time, in seconds, for the system to decay by 1/e,
// or to reach 1-1/e of its final value.
//
// #### Reference
// <https://ccrma.stanford.edu/~jos/mdft/Exponentials.html>
//-----------------------------------------------------------------------------
declare avg_tau author "Dario Sanfilippo and Julius O. Smith III";
declare avg_tau copyright "Copyright (C) 2020 Dario Sanfilippo
      <sanfilippo.dario@gmail.com> and
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare avg_tau license "MIT-style STK-4.3 license";
avg_tau(period, x) = fi.lptau(period, x);


//----------------------------`(fi.)avg_t60`-------------------------------------
// Averaging function based on a one-pole filter and the t60 response time.
// This response is particularly useful when the system is required to
// reach the final value after about `period` seconds.
//
// #### Usage
//
// ```
// _ : avg_t60(period) : _
// ```
//
// Where:
//
// * `period` is the time, in seconds, for the system to decay by 1/1000,
// or to reach 1-1/1000 of its final value.
//
// #### Reference
// <https://ccrma.stanford.edu/~jos/mdft/Audio_Decay_Time_T60.html>
//-----------------------------------------------------------------------------
declare avg_t60 author "Dario Sanfilippo and Julius O. Smith III";
declare avg_t60 copyright "Copyright (C) 2020 Dario Sanfilippo
      <sanfilippo.dario@gmail.com> and
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare avg_t60 license "MIT-style STK-4.3 license";
avg_t60(period, x) = fi.lpt60(period, x);


//----------------------------`(fi.)avg_t19`-------------------------------------
// Averaging function based on a one-pole filter and the t19 response time.
// This response is close to the moving-average algorithm as it roughly reaches
// the final value after `period` seconds and shows about the same
// oscillations for non-steady signals.
//
// #### Usage
//
// ```
// _ : avg_t19(period) : _
// ```
//
// Where:
//
// * `period` is the time, in seconds, for the system to decay by 1/e^2.2,
// or to reach 1-1/e^2.2 of its final value.
//
// #### Reference
// Zölzer, U. (2008). Digital audio signal processing (Vol. 9). New York: Wiley.
//-----------------------------------------------------------------------------
declare avg_t19 author "Dario Sanfilippo and Julius O. Smith III";
declare avg_t19 copyright "Copyright (C) 2020 Dario Sanfilippo
      <sanfilippo.dario@gmail.com> and
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare avg_t19 license "MIT-style STK-4.3 license";
avg_t19(period, x) = fi.lpt19(period, x);


//============================Kalman Filters===================================
//=============================================================================
//
// Functions related to the Kalman filter.

kalmanEnv = environment {

    matMul = la.matMul;

    // Predicts the next state.
    // * `N`: State size
    // * `M`: Measurement size
    // * `F`: State transition matrix (NxN)
    // * `B`: Control input matrix (NxM)
    // Implicit arguments:
    // * `x`: Current state (Nx1)
    // * `u`: Control input (Mx1)
    predictState(N, M, F, B) = si.vecOp((Fx, Bu), +)
    with {
        // Fx = F @ x, where F is NxN and x is Nx1
        Fx = matMul(N, N, N, 1, F, si.bus(N));
        
        // Bu = B @ u, where B is NxM and u is Mx1
        Bu = matMul(N, M, M, 1, B, si.bus(M));
    };

    // Predicts the covariance matrix for the next state.
    // * `N`: State size
    // * `F`: State transition matrix (NxN)
    // * `Q`: Process noise covariance matrix (NxN)
    // * `P`: Current state covariance matrix (NxN)
    // No implict arguments.
    predictCovariance(N, F, Q, P) = si.vecOp((FPF_T, Q), +)
    with {
        // FP = F @ P
        FP = matMul(N, N, N, N, F, P);
        
        // FPF_T = FP @ F^T, where FP is NxN and F^T is NxN
        FPF_T = matMul(N, N, N, N, FP, la.transpose2(N, N, F));
    };

    // Calculates the Kalman gain (NxM).
    // * `N`: State size
    // * `M`: Measurement size
    // * `H`: Observation matrix (MxN)
    // * `R`: Measurement noise covariance matrix (MxM)
    // Implict arguments:
    // * `P_pred`: Predicted covariance matrix (NxN)
    kalmanGain(N, M, H, R) = P_pred <: 
        matMul(N, M, M, M, matMul(N, N, N, M, P_pred, H_T), S_inv)
    with {
        P_pred = si.bus(N*N);
        // H_T = Transpose of H
        H_T = la.transpose2(M, N, H);
        
        // HP = H * P_pred, where H is MxN and P_pred is NxN
        HP = matMul(M, N, N, N, H, P_pred);
        
        // HPH_T = HP @ H^T, where HP is MxN and H^T is NxM, resulting in MxM
        HPH_T = matMul(M, N, N, M, HP, H_T);
        
        // S_inv = inverse(HPH_T + R)
        S_inv = si.vecOp((HPH_T, R), +) : la.inverse(M);
    };

    // Updates the state estimate based on the Kalman gain and measurement.
    // * `N`: State size
    // * `M`: Measurement size
    // * `H`: Observation matrix (MxN)
    // Implicit arguments:
    // * `K`: Kalman gain matrix (NxM)
    // * `x_pred`: (Nx1)
    // * `z`: (Mx1)
    updateState(N, M, H) = si.bus(N*M), si.bus(N), si.bus(M) <: 
        si.vecOp((get_x_pred, K_mul_residual(N, M, H)), +)
    with {
        get_x_pred = par(i, N*M, !), par(i, N, _), par(i, M, !);
    };

    // todo: find a way to not use mySwap.
    // Its purpose is to route K, x, z into K, z, x
    mySwap(N, M) = si.bus(N*M), si.bus(N), si.bus(M) <: 
        select_K, select_z, select_x_pred
    with {
        select_K = par(i, N*M, _), par(i, N, !), par(i, M, !);
        select_x_pred = par(i, N*M, !), par(i, N, _), par(i, M, !);
        select_z = par(i, N*M, !), par(i, N, !), par(i, M, _);
    };

    K_mul_residual(N, M, H) = mySwap(N, M) : matMul(N, M, M, 1, si.bus(N*M), residual)
    with {
        // residual = z - H @ x_pred, where z is Mx1 and H @ x_pred is Mx1
        residual = si.vecOp((si.bus(M), matMul(M, N, N, 1, H, si.bus(N))), -);
    };

    // Updates the covariance matrix based on the Kalman gain and observation matrix.
    // * `N`: State size
    // * `M`: Measurement size
    // * `H`: Observation matrix (MxN)
    // * `K`: Kalman gain matrix (NxM)
    // * `P_pred`: Predicted covariance matrix (NxN)
    // No implict arguments.
    updateCovariance(N, M, H, K, P_pred) = 
        matMul(N, N, N, N, I_minus_KH, P_pred)
    with {
        // KH = K @ H, where K is NxM and H is MxN
        KH = matMul(N, M, M, N, K, H);
            
        // I_minus_KH = I - K @ H, where both are NxN
        I_minus_KH = si.vecOp((la.identity(N), KH), -);
    };

    // implict arguments:
    // * `u`: Control input (Mx1)
    // * `z`: Measurement signal (Mx1)
    kalmanFilter(N, M, B, R, H, Q, F, reset) = (p, x, u, z <: pNew, xNew) ~ (pBus, xBus) : pCut, xBus
    with {
        doReset = reset > 0;

        pBus = si.bus(N*N);
        pCut = par(i, N*N, !);
        pInit = la.diag(N, par(i, N, 100));  // large value for initial uncertainty in P.
        p = pBus, pInit : ba.selectbus(N*N, 2, doReset);

        xBus = si.bus(N);
        xCut = par(i, N, !);
        xInit = par(i, N, 0);  // initialize with zeros.
        x = xBus, xInit : ba.selectbus(N, 2, doReset);

        u = si.bus(M);
        z = si.bus(M);

        uCut = par(i, M, !);
        zCut = par(i, M, !);

        predCov = predictCovariance(N, F, Q, pBus);
        gain = kalmanGain(N, M, H, R);

        xNew = pBus, xBus, u, z : updateState(N, M, H, gain(predCov), predictState(N, M, F, B, xBus, u), z);        
        pNew = pBus, xCut, uCut, zCut : predCov <: updateCovariance(N, M, H, gain(pBus), pBus);
    };
};

//-----------`(fi.)kalman`----------------------------------------------------
// The Kalman filter. It returns the state (a bus of size `N`).
// Note that the only compile-time constant arguments are `N` and `M`.
// Other arguments are capitalized because they're matrices, and it makes
// reading them much easier.
//
// #### Usage
// ```
// kalman(N, M, B, R, H, Q, F, reset, u, z) : si.bus(N)
// ```
//
// Where:
//
// * `N`: State size (constant int)
// * `M`: Measurement size (constant int)
// * `B`: Control input matrix (NxM)
// * `R`: Measurement noise covariance matrix (MxM)
// * `H`: Observation matrix (MxN)
// * `Q`: Process noise covariance matrix (NxN)
// * `F`: State transition matrix (NxN)
// * `reset`: Reset trigger. Whenever `reset>0`, the internal state `x` and covariance matrix `P` are reset.
// * `u`: Control input (Mx1)
// * `z`: Measurement signal (Mx1)
//
// #### Example test programs
// Demo 1 `(N=1, M=1)` (don't listen, just use oscilloscope):
//
// ```
// process = fi.kalman(N, M, B, R, H, Q, F, reset, u, z) : it.interpolate_linear(filteredAmt, z)
// with {
//    B = 1.;
//    R = 0.1;
//    H = 1;
//    Q = .01; 
//    F = la.identity(N);
//    reset = button("reset");
//
//    // Dimensions
//    N = 1; // State size
//    M = 1; // Measurement size
//
//    freq = hslider("Freq", 1, 0.01, 10, .01);
//    u = 0.; // constant input
//    trueState = os.osc(freq)*.5 + u;
//    noiseGain = hslider("Noise Gain", .1, 0, 1, .01);
//
//    filteredAmt = hslider("Filter Amount", 1, 0, 1, .01) : si.smoo;
//
//    measurementNoise = no.noise*noiseGain;
//    z = trueState + measurementNoise; // Observed state
// };
// ```
//
// Demo 2 `(N=2, M=1)` (don't listen, just use oscilloscope)
//
// ```
// process = fi.kalman(N, M, B, R, H, Q, F, reset, u, z)
// with {
//     B = par(i, N, 0);
//     R = (0.1);
//     H = (1, 0);
//     Q = la.diag(2, par(i, N, .1));
//     F = la.identity(N);
//     reset = 0;
//     u = si.bus(M);
//     z = si.bus(M);
// 
//     // Dimensions
//     N = 2; // State size
//     M = 1; // Measurement size
// };
// ```
// #### References
// * <https://en.wikipedia.org/wiki/Kalman_filter>
// * <https://www.cs.unc.edu/~welch/kalman/index.html>
//----------------------------------------------------------------------------
declare kalman author "David Braun";
declare kalman license "MIT License";
kalman = kalmanEnv.kalmanFilter;


/*******************************************************************************
# Licenses

## STK 4.3 License

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
of the Software, and to permit persons to whom the Software is furnished to do
so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

Any person wishing to distribute modifications to the Software is asked to send
the modifications to the original developer so that they can be incorporated
into the canonical version.  For software copyrighted by Julius O. Smith III,
email your modifications to <jos@ccrma.stanford.edu>.  This is, however, not a
binding provision of this license.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

--------------------------------------------------------------------------------

## LGPL 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.

--------------------------------------------------------------------------------

## ISC License

Permission to use, copy, modify, and/or distribute this software for any
purpose with or without fee is hereby granted, provided that the above
copyright notice and this permission notice appear in all copies.

THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.

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