//##################################### filters.lib ########################################
// Filters library. Its official prefix is `fi`.
//
// This library provides a comprehensive collection of linear and nonlinear
// filters used in audio and signal processing. It includes low-pass, high-pass,
// band-pass, allpass, shelving, equalizer, and crossover filters, as well as advanced
// analog and digital filter design sections for both educational and production use.
//
// The Filters library is organized into 24 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)
// * [Adaptive Filters](#adaptive-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.9.0";

//===============================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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// zero_test = src : fi.zero(0.5);
// ```
//
// #### References
// <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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// pole_test = src : fi.pole(0.9);
// ```
//
// #### References
// <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 "LicenseRef-STK-4.3";
pole(p) = + ~ *(p);

//----------------------`(fi.)integrator`--------------------------
// Same as `pole(1)` [implemented separately for block-diagram clarity].
//
// #### Usage
//
// ```
// _ : integrator : _
// ```
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// integrator_test = src : fi.integrator;
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// dcblockerat_test = src : fi.dcblockerat(30);
// ```
//
// #### 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 "LicenseRef-STK-4.3";
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 a standard Faust function.
//
// #### Usage
//
// ```
// _ : dcblocker : _
// ```
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// dcblocker_test = src : fi.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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lptN_test = src : fi.lptN(60, 0.1);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
lptN(N, tN, x) = x : si.smooth(ba.tau2pole(tN / log(10.0^(float(N)/20.0))));
// Special cases of lptN
//------------------------`(fi.)lptau`--------------------------
// One-pole lowpass with a tau time constant (1/e attenuation after `tN` seconds).
//
// #### Usage
//
// ```
// _ : lptau(tN) : _
// ```
//
// Where:
//
// * `tN`: tau time constant in seconds (1/e attenuation after `tN` seconds)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lptau_test = src : fi.lptau(0.1);
// ```
//------------------------------------------------------------
lptau(tN, x) = lptN(8.6858896381, tN, x); // Tau time constant, i.e., 1/e atten. after tN secs

//------------------------`(fi.)lpt60`--------------------------
// One-pole lowpass with a T60 time constant (60 dB attenuation after `tN` seconds).
//
// #### Usage
//
// ```
// _ : lpt60(tN) : _
// ```
//
// Where:
//
// * `tN`: T60 time constant in seconds (60 dB attenuation after `tN` seconds)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lpt60_test = src : fi.lpt60(0.3);
// ```
//------------------------------------------------------------
lpt60(tN, x) = lptN(60, tN, x); // T60 constant, i.e., 1/1000 atten. after tN secs

//------------------------`(fi.)lpt19`--------------------------
// One-pole lowpass with a T19 time constant (approx. 19 dB attenuation after `tN` seconds).
//
// #### Usage
//
// ```
// _ : lpt19(tN) : _
// ```
//
// Where:
//
// * `tN`: T19 time constant in seconds (approximately 19 dB attenuation after `tN` seconds)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lpt19_test = src : fi.lpt19(0.2);
// ```
//------------------------------------------------------------
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// ff_comb_test = src : fi.ff_comb(2048, 64, 1, 0.7);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// ff_fcomb_test = src : fi.ff_fcomb(2048, 64.5, 1, 0.7);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
ff_fcomb(maxdel,del,b0,bM) = _ <: *(b0), bM * de.fdelay(maxdel,del) : +;

//-----------`(fi.)ffcombfilter`-------------------
// Typical special case of `ff_comb()` where: `b0 = 1`.
//
// #### Usage
//
// ```
// _ : ffcombfilter(maxdel,del,g) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `del`: current (float) comb-filter delay between 0 and maxdel
// * `g`: gain applied to the delayed tap
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// ffcombfilter_test = src : fi.ffcombfilter(2048, 64, 0.7);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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.
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// fb_comb_common_test = src : fi.fb_comb_common(@, 64, 0.8, 0.6);
// ```
// --------------------------------------------------------
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// fb_comb_test = src : fi.fb_comb(2048, 64, 0.7, 0.6);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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).
// `fb_fcomb` is a standard Faust function.
//
// #### 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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// fb_fcomb_test = src : fi.fb_fcomb(2048, 64.5, 0.7, 0.6);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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.
//
// #### Usage
//
// ```
// _ : rev1(maxdel,N,g) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `N`: current feedback comb-filter delay in samples
// * `g`: feedback gain (its sign is negated internally)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// rev1_test = src : fi.rev1(2048, 64, 0.6);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// fbcombfilter_test = src : fi.fbcombfilter(2048, 64, 0.6);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// allpass_comb_test = src : fi.allpass_comb(2048, 64, 0.6);
// ```
//
// #### 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 "LicenseRef-STK-4.3";
allpass_comb(maxdel,N,aN) = (+ <: de.delay(maxdel,N-1),*(aN)) ~ *(-aN) : mem,_ : +;


//-------------------`(fi.)allpass_fcomb`-----------------
// Schroeder Allpass Comb Filter.
// `allpass_fcomb` is a standard Faust function.
// 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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// allpass_fcomb_test = src : fi.allpass_fcomb(2048, 64.5, 0.6);
// ```
//
// #### 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 "LicenseRef-STK-4.3";
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.
//
// #### Usage
//
// ```
// _ : rev2(maxlen,len,g) : _
// ```
//
// Where:
//
// * `maxlen`: maximum delay (a power of 2)
// * `len`: current allpass comb-filter delay in samples
// * `g`: allpass coefficient (its sign is negated internally)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// rev2_test = src : fi.rev2(2048, 64, 0.6);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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);`)
//
// #### Usage
//
// ```
// _ : allpass_fcomb5(maxdel,N,aN) : _
// _ : allpass_fcomb1a(maxdel,N,aN) : _
// ```
//
// Where:
//
// * `maxdel`: maximum delay (a power of 2)
// * `N`: current (float) allpass comb-filter delay between 0 and maxdel
// * `aN`: coefficient of z^(-N) in the transfer function numerator
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// allpass_fcomb5_test = src : fi.allpass_fcomb5(2048, 64.5, 0.6);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// iir_test = src : fi.iir((0.5, 0.5), (0.3));
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
iir(bv,av) = ma.sub ~ fir(av) : fir(bv);

//-----------------------------`(fi.)fir`---------------------------------
// FIR filter (convolution of FIR filter coefficients with a signal).
// `fir` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : fir(bv) : _
// ```
//
// Where:
//
// * `bv`: `b0,b1,...,bn`, 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);
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// fir_test = src : fi.fir((0.2, 0.2, 0.2, 0.2, 0.2));
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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,ba.take(i+1,kv) * x@i); // take() defined in basics.lib
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// si = library("signals.lib");
// src = os.osc(440);
// convN_test = (src <: si.bus(3)) : fi.convN(3, (0.3, 0.2, 0.1, 0.05));
// ```

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 "LicenseRef-STK-4.3";
convN(N,kv) = sum(i,N, @(i)*ba.take(i+1,kv)); // take() defined in basics.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 "LicenseRef-STK-4.3";
conv(kv) = fir(kv);

//----------------`(fi.)tf1`, `(fi.)tf2` and `(fi.)tf3`----------------------
// tfN = N'th-order direct-form digital filter.
// `tf2` is a standard Faust function.
//
// #### 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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// tf1_test = src : fi.tf1(0.5, 0.25, -0.4);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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 "LicenseRef-STK-4.3";
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 "LicenseRef-STK-4.3";
tf3(b0,b1,b2,b3,a1,a2,a3) = iir((b0,b1,b2,b3),(a1,a2,a3));

//-----------------------------`(fi.)TF2`--------------------------------
// Tunable second-order direct-form digital filter, the "original"
// biquad of music.lib: y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] -
// a1*y[n-1] - a2*y[n-2]. Kept for comparison and for the compatibility
// libraries (maxmsp.lib, instruments.lib); new code should use `tf2`
// instead.
//
// #### Usage
//
// ```
// _ : TF2(b0,b1,b2,a1,a2) : _
// ```
//
// Where:
//
// * `b0`, `b1`, `b2`: the feedforward coefficients
// * `a1`, `a2`: the feedback coefficients
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// TF2_legacy_test = os.osc(440) : fi.TF2(0.2, 0.4, 0.2, -0.5, 0.3);
// ```
//------------------------------------------------------------
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// notchw_test = src : fi.notchw(200, 1000);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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
//
// #### References
//
// <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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// tf21_test = src : fi.tf21(0.1, 0.2, 0.1, -0.5, 0.06);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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 "LicenseRef-STK-4.3";
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 "LicenseRef-STK-4.3";
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 "LicenseRef-STK-4.3";
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.

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

//-------------------------------`(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)`.
//
// #### Test
// ```
// fi = library("filters.lib");
// si = library("signals.lib");
// av2sv_test = fi.av2sv((-0.4, 0.1)) : si.bus(2);
// ```
//
// #### References
//
//   <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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// si = library("signals.lib");
// bvav2nuv_test = fi.bvav2nuv((0.1, 0.2, 0.3), (-0.4, 0.1)) : si.bus(3);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// iir_lat2_test = src : fi.iir_lat2((0.1, 0.2, 0.3), (-0.4, 0.1));
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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.
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// si = library("signals.lib");
// src = os.osc(440);
// allpassnt_test = src : fi.allpassnt(2, (0.3, -0.2)) : si.bus(3);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// iir_kl_test = src : fi.iir_kl((0.1, 0.2, 0.3), (-0.4, 0.1));
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// si = library("signals.lib");
// src = os.osc(440);
// allpassnklt_test = src : fi.allpassnklt(2, (0.3, -0.2)) : si.bus(3);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// iir_lat1_test = src : fi.iir_lat1((0.1, 0.2, 0.3), (-0.4, 0.1));
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// si = library("signals.lib");
// src = os.osc(440);
// allpassn1mt_test = src : fi.allpassn1mt(2, (0.3, -0.2)) : si.bus(3);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// iir_nl_test = src : fi.iir_nl((0.1, 0.2, 0.3), (-0.4, 0.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 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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// si = library("signals.lib");
// src = os.osc(440);
// allpassnnlt_test = src : fi.allpassnnlt(2, (0.3, -0.2)) : si.bus(3);
// ```
//
// #### 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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// tf2np_test = src : fi.tf2np(0.6, 0.3, 0.2, -0.5, 0.2);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";

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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// wgr_test = fi.wgr(440, 0.995, src);
// ```
//
// #### 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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// nlf2_test = fi.nlf2(440, 0.995, src);
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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`: first allpass coefficient
// * `a2`: second allpass coefficient
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// apnl_test = fi.apnl(0.5, -0.5, src);
// ```
//
// #### References
//
// * "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 "LicenseRef-STK-4.3";
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); };
// };
// ```
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// dual_src = os.osc(440), os.osc(660);
// scatN_test = dual_src : fi.scatN(2, (1, 1), _);
// ```
//
// #### 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 "LicenseRef-STK-4.3";

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) : +~_;
// ```
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// scat_test = src : fi.scat(0.5, _);
// ```
//
// #### References
//
// * <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 "LicenseRef-STK-4.3";

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 `(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`.
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// allpassn_test = src : fi.allpassn(3, (0.3, 0.2, 0.1));
// ```
//
// #### 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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// allpassnn_test = src : fi.allpassnn(3, (0.3, 0.2, 0.1));
// ```
//----------------------------------------------
// 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 "LicenseRef-STK-4.3";

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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// allpassnkl_test = src : fi.allpassnkl(3, (0.3, 0.2, 0.1));
// ```
//----------------------------------------------
// Kelly-Lochbaum:
declare allpassnkl author "Julius O. Smith III";
declare allpassnkl copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare allpassnkl license "LicenseRef-STK-4.3";

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

//---------------`(fi.)allpassn1m`-----------------
// 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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// allpassn1m_test = src : fi.allpassn1m(3, (0.3, 0.2, 0.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 "LicenseRef-STK-4.3";

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.
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// ma = library("maths.lib");
// src = os.osc(440);
// tf2s_test = src : fi.tf2s(0, 0, 1, sqrt(2), 1, ma.PI*ma.SR/2);
// ```
//
// #### References
//
// * <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 "LicenseRef-STK-4.3";
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;
};

//----------------------------------------------
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// ma = library("maths.lib");
// src = os.osc(440);
// tf2snp_test = src : fi.tf2snp(0, 0, 1, sqrt(2), 1, ma.PI*ma.SR/2);
// ```
//----------------------------------------------

//-----------------------------`(fi.)tf1snp`-------------------------------
// First-order special case of tf2snp above.
//
// #### Usage
//
// ```
// _ : tf1snp(b1,b0,a0,w1) : _
// ```
//
// Where:
//
// * `b1`: analog numerator coefficient for `s`
// * `b0`: analog numerator constant coefficient
// * `a0`: analog denominator constant coefficient
// * `w1`: desired digital frequency (radians/second) corresponding to analog frequency `1 rad/sec`
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// ma = library("maths.lib");
// src = os.osc(440);
// tf1snp_test = src : fi.tf1snp(0, 1, 1, ma.PI*ma.SR/2);
// ```
//----------------------------------------------
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 "LicenseRef-STK-4.3";
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) : _
// ```
//
// Where:
//
// * `b3`: analog numerator coefficient of `s^3`
// * `b2`: analog numerator coefficient of `s^2`
// * `b1`: analog numerator coefficient of `s`
// * `b0`: analog numerator constant term
// * `a3`: analog denominator coefficient of `s^3`
// * `a2`: analog denominator coefficient of `s^2`
// * `a1`: analog denominator coefficient of `s`
// * `a0`: analog denominator constant term
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// tf3slf_test = src : fi.tf3slf(0, 0, 0, 1, 1, 2, 2, 1);
// ```
//----------------------------------------------
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 "LicenseRef-STK-4.3";
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`: analog numerator coefficient of `s`
// * `b0`: analog numerator constant term
// * `a0`: analog denominator constant term in `s + a0`
// * `w1`: desired digital frequency in radians per second corresponding to analog frequency `1 rad/s`
//
// 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.
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// ma = library("maths.lib");
// src = os.osc(440);
// tf1s_test = src : fi.tf1s(0, 1, 1, ma.PI*ma.SR/2);
// ```
//
// #### References
//
// * <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 "LicenseRef-STK-4.3";
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) : _
// ```
//
// Where:
//
// * `b2`, `b1`, `b0`: analog lowpass numerator coefficients
// * `a1`, `a0`: analog lowpass denominator coefficients
// * `w1`: half the desired passband width in radians/second
// * `wc`: desired center frequency in radians/second
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// ma = library("maths.lib");
// src = os.osc(440);
// tf2sb_test = src : fi.tf2sb(0, 0, 1, sqrt(2), 1, 2*ma.PI*200, 2*ma.PI*1000);
// ```
//----------------------------------------------
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 "LicenseRef-STK-4.3";
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) : _
// ```
//
// Where:
//
// * `b1`, `b0`: analog numerator coefficients
// * `a0`: analog denominator constant coefficient
// * `w1`: half the desired passband width in radians/second
// * `wc`: desired center frequency in radians/second
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// ma = library("maths.lib");
// src = os.osc(440);
// tf1sb_test = src : fi.tf1sb(0, 1, 1, 2*ma.PI*200, 2*ma.PI*1000);
// ```
//----------------------------------------------
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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// resonlp_test = src : fi.resonlp(1000, 2, 0.8);
// ```
//---------------------------------------------------------------------
// 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 "LicenseRef-STK-4.3";

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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// resonhp_test = fi.resonhp(1000, 2, 0.8, src);
// ```
//---------------------------------------------------------------------
// 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 "LicenseRef-STK-4.3";

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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// resonbp_test = src : fi.resonbp(1000, 2, 0.8);
// ```
//---------------------------------------------------------------------
// 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 "LicenseRef-STK-4.3";

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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lowpass_test = src : fi.lowpass(4, 2000);
// ```
//
// #### 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 "LicenseRef-STK-4.3";

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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// highpass_test = src : fi.highpass(4, 500);
// ```
//
// #### 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 highpass license "LicenseRef-STK-4.3";

highpass(N,fc) = lowpass0_highpass1(1,N,fc);


//-------------`(fi.)lowpass0_highpass1`--------------
// Generic Butterworth lowpass/highpass filter: the shared implementation
// behind `fi.lowpass` and `fi.highpass`, with a selector choosing between
// the two responses.
//
// #### Usage
//
// ```
// _ : lowpass0_highpass1(s,N,fc) : _
// ```
//
// Where:
//
// * `s`: response selector: 0 for lowpass, 1 for highpass (a constant numerical expression)
// * `N`: filter order (a constant numerical expression)
// * `fc`: -3dB cutoff frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lowpass0_highpass1_test = src : fi.lowpass0_highpass1(0, 2, 1000);
// ```
//------------------------------
declare lowpass0_highpass1 author "Julius O. Smith III";
declare lowpass0_highpass1 copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare lowpass0_highpass1 license "LicenseRef-STK-4.3";

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.)highpass_plus_lowpass`----------------
// Sum of the order-`N` Butterworth highpass and lowpass responses at the same
// cutoff: an allpass used as delay equalizer by `fi.filterbank` and the
// filter banks below. For odd orders 1, 3 and 5, canceling pole-zero pairs
// are removed.
//
// #### Usage
//
// ```
// _ : highpass_plus_lowpass(N,fc) : _
// ```
//
// Where:
//
// * `N`: filter order (a constant numerical expression)
// * `fc`: crossover frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// highpass_plus_lowpass_test = os.osc(440) : fi.highpass_plus_lowpass(3, 1000);
// ```
//------------------------------
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 "LicenseRef-STK-4.3";

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

//--------------------`(fi.)highpass_minus_lowpass`----------------
// Difference between the order-`N` Butterworth highpass and lowpass responses
// at the same cutoff: an allpass used as delay equalizer by the
// dc-inverted filter banks (`fi.mth_octave_filterbank_alt`). For odd orders
// 3 and 5, canceling pole-zero pairs are removed.
//
// #### Usage
//
// ```
// _ : highpass_minus_lowpass(N,fc) : _
// ```
//
// Where:
//
// * `N`: filter order (a constant numerical expression)
// * `fc`: crossover frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// highpass_minus_lowpass_test = os.osc(440) : fi.highpass_minus_lowpass(3, 1000);
// ```
//------------------------------
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 "LicenseRef-STK-4.3";
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);
};

//--------------------`(fi.)highpass_plus_lowpass_even`----------------
// Sum of the even-order Butterworth highpass and lowpass responses,
// `highpass(N,fc) + lowpass(N,fc)`. Takes two copies of the input signal
// (both normally fed the same signal); `fi.highpass_plus_lowpass` selects
// this case automatically for even `N`.
//
// #### Usage
//
// ```
// _,_ : highpass_plus_lowpass_even(N,fc) : _
// ```
//
// Where:
//
// * `N`: even filter order (a constant numerical expression)
// * `fc`: crossover frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// highpass_plus_lowpass_even_test = os.osc(440), os.osc(440) : fi.highpass_plus_lowpass_even(4, 1000);
// ```
//------------------------------
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 "LicenseRef-STK-4.3";
highpass_plus_lowpass_even(N,fc) = highpass(N,fc) + lowpass(N,fc);

//--------------------`(fi.)highpass_minus_lowpass_even`----------------
// Difference between the even-order Butterworth highpass and lowpass
// responses, `highpass(N,fc) - lowpass(N,fc)`. Takes two copies of the input
// signal (both normally fed the same signal); `fi.highpass_minus_lowpass`
// selects this case automatically for even `N`.
//
// #### Usage
//
// ```
// _,_ : highpass_minus_lowpass_even(N,fc) : _
// ```
//
// Where:
//
// * `N`: even filter order (a constant numerical expression)
// * `fc`: crossover frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// highpass_minus_lowpass_even_test = os.osc(440), os.osc(440) : fi.highpass_minus_lowpass_even(4, 1000);
// ```
//------------------------------
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_minus_lowpass_even license "LicenseRef-STK-4.3";
highpass_minus_lowpass_even(N,fc) = highpass(N,fc) - lowpass(N,fc);


//--------------------`(fi.)highpass_plus_lowpass_odd`-----------------
// Sum of the odd-order highpass and lowpass responses.
//
// #### Usage
//
// ```
// _,_ : highpass_plus_lowpass_odd(N,fc) : _
// ```
//
// Where:
//
// * `N`: odd Butterworth order
// * `fc`: cutoff frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// highpass_plus_lowpass_odd_test = os.osc(440), os.osc(440) : fi.highpass_plus_lowpass_odd(3, 1000);
// ```
// FIXME: Rewrite the following, as for orders 3 and 5 above,
//        to eliminate pole-zero cancellations:
//------------------------------
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 "LicenseRef-STK-4.3";
highpass_plus_lowpass_odd(N,fc) = highpass(N,fc) + lowpass(N,fc);

//--------------------`(fi.)highpass_minus_lowpass_odd`----------------
// Difference between the odd-order highpass and lowpass responses.
//
// #### Usage
//
// ```
// _,_ : highpass_minus_lowpass_odd(N,fc) : _
// ```
//
// Where:
//
// * `N`: odd Butterworth order
// * `fc`: cutoff frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// highpass_minus_lowpass_odd_test = os.osc(440), os.osc(440) : fi.highpass_minus_lowpass_odd(3, 1000);
// ```
// 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);
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_minus_lowpass_odd license "LicenseRef-STK-4.3";


//==========================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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lowpass3e_test = src : fi.lowpass3e(1000);
// ```
//
// #### 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 "LicenseRef-STK-4.3";

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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lowpass6e_test = src : fi.lowpass6e(1000);
// ```
//
// #### 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 "LicenseRef-STK-4.3";

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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// highpass3e_test = src : fi.highpass3e(1000);
// ```
//-------------------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// highpass6e_test = src : fi.highpass6e(1000);
// ```
//-------------------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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`.
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// bandpass_test = src : fi.bandpass(2, 500, 1500);
// ```
//-------------------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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`.
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// bandstop_test = src : fi.bandstop(2, 500, 1500);
// ```
//-------------------------------------------------------------------------
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 "LicenseRef-STK-4.3";
bandstop(Nh,fl,fu) = bandpass0_bandstop1(1,Nh,fl,fu);

//--------------------`(fi.)bandpass0_bandstop1`----------------
// Generic Butterworth bandpass/bandstop filter: the shared implementation
// behind `fi.bandpass` and `fi.bandstop`, with a selector choosing between
// the two responses.
//
// #### Usage
//
// ```
// _ : bandpass0_bandstop1(s,Nh,fl,fu) : _
// ```
//
// Where:
//
// * `s`: response selector: 0 for bandpass, 1 for bandstop (a constant numerical expression)
// * `Nh`: HALF the filter order, i.e. the number of second-order sections (a constant numerical expression)
// * `fl`: lower -3dB band edge frequency in Hz
// * `fu`: upper -3dB band edge frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// bandpass0_bandstop1_test = src : fi.bandpass0_bandstop1(0, 2, 500, 1500);
// ```
//-------------------------------------------------
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 "LicenseRef-STK-4.3";
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)`.
//
// #### Usage
//
// ```
// _ : bandpass6e(fl,fu) : _
// ```
//
// Where:
//
// * `fl`: lower -3dB band edge frequency in Hz
// * `fu`: upper -3dB band edge frequency in Hz
//--------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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)`.
//
// #### Usage
//
// ```
// _ : bandpass12e(fl,fu) : _
// ```
//
// Where:
//
// * `fl`: lower -3dB band edge frequency in Hz
// * `fu`: upper -3dB band edge frequency in Hz
//--------------------------------------------------------------
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 "LicenseRef-STK-4.3";
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 "LicenseRef-STK-4.3";
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
};

//-----------------------------`(fi.)pospass6e`--------------------------------
// Positive-pass filter like `pospass`, but built on the order-6 elliptic
// lowpass `lowpass6e` instead of a Butterworth lowpass: a steeper
// transition band for the same order.
//
// #### Usage
//
// ```
// _ : pospass6e(fc) : _,_
// ```
//
// Where:
//
// * `fc`: lower cutoff frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// pospass6e_test = os.osc(440) : fi.pospass6e(100);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";
pospass6e(fc) = fi._pospass0(lpf) with {
  lpf = fi.lowpass6e(ma.SR/4 - fc); // Elliptic lowpass, order 6
};

_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
};

//------------------------`(fi.)hilbert`---------------------------
// Approximate Hilbert transform: twice the imaginary part of the analytic
// signal computed by `fi.pospass`, i.e. a signal in quadrature (-90 degrees)
// with the in-phase (real) output of the same filter. This is the one-line
// construction given in the `fi.pospass` documentation.
//
// The output carries the group delay of the underlying filter, so pair it
// with the real part of `fi.pospass` (not with the raw input) when phase
// alignment matters. Accuracy improves with the order `N` and with the input
// frequency: rejecting the negative-frequency image of a component at `f` Hz
// requires the internal half-band lowpass to cut inside a transition band of
// width `2*f`, so low frequencies demand high orders (measured at 48 kHz:
// with `N = 8`, a 5 kHz sine yields an analytic pair with 0.5% envelope
// ripple and quadrature to 1e-4, while a 1 kHz sine needs `N` around 32 for
// comparable quality).
//
// #### Usage
//
// ```
// _ : hilbert(N,fc) : _
// ```
//
// Where:
//
// * `N`: order of the underlying Butterworth positive-pass filter (a constant numerical expression)
// * `fc`: lower cutoff frequency in Hz of the positive-pass filter
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// hilbert_test = os.osc(440) : fi.hilbert(4, 20);
// ```
//
// #### References
//
// * <https://ccrma.stanford.edu/~jos/mdft/Analytic_Signals_Hilbert_Transform.html>
// * <https://ccrma.stanford.edu/~jos/sasp/Hilbert_Transform.html>
//------------------------------------------------------------
declare hilbert author "Julius O. Smith III";
declare hilbert copyright "Copyright (C) 2003-2019 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare hilbert license "LicenseRef-STK-4.3";
hilbert(N,fc) = fi.pospass(N,fc) : !,*(2.0);


//=================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.)lowshelf`----------------------
// 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).
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lowshelf_test = src : fi.lowshelf(3, 6, 500);
// ```
//
// #### 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 "LicenseRef-STK-4.3";
lowshelf(N,L0,fx) = filterbank(N,(fx)) : _, *(ba.db2linear(L0)) :> _;

// Special cases and optimization:

//----------------------`(fi.)low_shelf`----------------------
// Default low shelf filter: 3rd-order Butterworth case of `fi.lowshelf`
// (`low_shelf = lowshelf(3)`), boosting or cutting by `L0` dB between dc and
// `fx`. `low_shelf` is a standard Faust function. See `fi.lowshelf` for the
// full documentation, including the shelf shape.
//
// #### Usage
//
// ```
// _ : low_shelf(L0,fx) : _
// ```
//
// Where:
//
// * `L0`: desired level (dB) between dc and `fx` (boost `L0>0` or cut `L0<0`)
// * `fx`: transition frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// low_shelf_test = src : fi.low_shelf(6, 500);
// ```
//-----------------------------------------------
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 "LicenseRef-STK-4.3";
low_shelf  = lowshelf(3); // default = 3rd order Butterworth


//-----------------------`(fi.)low_shelf1`-----------------------
// First-order low shelf filter, an optimized special case of `fi.lowshelf(1)`:
// boosts or cuts by `L0` dB between dc and `fx`.
//
// #### Usage
//
// ```
// _ : low_shelf1(L0,fx) : _
// ```
//
// Where:
//
// * `L0`: desired level (dB) between dc and `fx` (boost `L0>0` or cut `L0<0`)
// * `fx`: transition frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// low_shelf1_test = fi.low_shelf1(2, 500, src);
// ```
//-----------------------------------------------
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 "LicenseRef-STK-4.3";
low_shelf1(L0,fx,x) = x + (ba.db2linear(L0)-1)*lowpass(1,fx,x); // optimized


//----------------------`(fi.)low_shelf1_l`----------------------
// First-order low shelf filter like `fi.low_shelf1`, but taking a LINEAR gain
// `G0` instead of a level in dB (avoids computing `ba.db2linear` at run time
// when the gain is already linear).
//
// #### Usage
//
// ```
// _ : low_shelf1_l(G0,fx) : _
// ```
//
// Where:
//
// * `G0`: desired linear gain between dc and `fx`
// * `fx`: transition frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// low_shelf1_l_test = fi.low_shelf1_l(2, 500, src);
// ```
//-----------------------------------------------
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 "LicenseRef-STK-4.3";
low_shelf1_l(G0,fx,x) = x + (G0-1)*lowpass(1,fx,x); // optimized


//----------------------`(fi.)lowshelf_other_freq`----------------------
// Second corner frequency of a low shelf: given the shelf order `N`, level
// `L0` and transition frequency `fx` of `fi.lowshelf(N,L0,fx)`, returns the
// frequency at which the response returns to 0 dB (see the "Shelf Shape"
// section of `fi.lowshelf`). This is a frequency computation, not a filter.
//
// #### Usage
//
// ```
// lowshelf_other_freq(N,L0,fx) : _
// ```
//
// Where:
//
// * `N`: shelf order (odd)
// * `L0`: shelf level in dB
// * `fx`: transition frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// lowshelf_other_freq_test = fi.lowshelf_other_freq(3, 6, 500);
// ```
//-----------------------------------------------
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 "LicenseRef-STK-4.3";
lowshelf_other_freq(N, L0, fx) = ba.db2linear(ba.linear2db(fx) + L0/N); // convenience

//--------------`(fi.)highshelf`--------------
// 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.
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// highshelf_test = src : fi.highshelf(3, 6, 2000);
// ```
//
// #### 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 "LicenseRef-STK-4.3";
highshelf(N,Lpi,fx) = filterbank(N,(fx)) : *(ba.db2linear(Lpi)), _ :> _;

// Special cases and optimization:

//----------------------`(fi.)high_shelf`----------------------
// Default high shelf filter: 3rd-order Butterworth case of `fi.highshelf`
// (`high_shelf = highshelf(3)`), boosting or cutting by `Lpi` dB between `fx`
// and SR/2. `high_shelf` is a standard Faust function. See `fi.highshelf`
// for the full documentation.
//
// #### Usage
//
// ```
// _ : high_shelf(Lpi,fx) : _
// ```
//
// Where:
//
// * `Lpi`: desired level (dB) between `fx` and SR/2 (boost `Lpi>0` or cut `Lpi<0`)
// * `fx`: transition frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// high_shelf_test = src : fi.high_shelf(6, 2000);
// ```
//---------------------------------------------------
high_shelf = highshelf(3); // default = 3rd order Butterworth

//----------------------`(fi.)high_shelf1`----------------------
// First-order high shelf filter, an optimized special case of
// `fi.highshelf(1)`: boosts or cuts by `Lpi` dB between `fx` and SR/2.
//
// #### Usage
//
// ```
// _ : high_shelf1(Lpi,fx) : _
// ```
//
// Where:
//
// * `Lpi`: desired level (dB) between `fx` and SR/2 (boost `Lpi>0` or cut `Lpi<0`)
// * `fx`: transition frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// high_shelf1_test = fi.high_shelf1(6, 2000, src);
// ```
//---------------------------------------------------
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 "LicenseRef-STK-4.3";
high_shelf1(Lpi,fx,x) = x + (ba.db2linear(Lpi)-1)*highpass(1,fx,x); // optimized

//----------------------`(fi.)high_shelf1_l`----------------------
// First-order high shelf filter like `fi.high_shelf1`, but taking a LINEAR
// gain `Gpi` instead of a level in dB.
//
// #### Usage
//
// ```
// _ : high_shelf1_l(Gpi,fx) : _
// ```
//
// Where:
//
// * `Gpi`: desired linear gain between `fx` and SR/2
// * `fx`: transition frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// high_shelf1_l_test = fi.high_shelf1_l(2, 2000, src);
// ```
//---------------------------------------------------
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 "LicenseRef-STK-4.3";
high_shelf1_l(Gpi,fx,x) = x + (Gpi-1)*highpass(1,fx,x); //optimized

// shelf transitions between frequency fx and this one:

//----------------------`(fi.)highshelf_other_freq`----------------------
// Second corner frequency of a high shelf: given the shelf order `N`, level
// `Lpi` and transition frequency `fx` of `fi.highshelf(N,Lpi,fx)`, returns
// the frequency at which the response leaves 0 dB. This is a frequency
// computation, not a filter.
//
// #### Usage
//
// ```
// highshelf_other_freq(N,Lpi,fx) : _
// ```
//
// Where:
//
// * `N`: shelf order (odd)
// * `Lpi`: shelf level in dB
// * `fx`: transition frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// highshelf_other_freq_test = fi.highshelf_other_freq(3, 6, 2000);
// ```
//---------------------------------------------------
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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// peak_eq_test = src : fi.peak_eq(6, 1000, 200);
// ```
//
// #### 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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// peak_eq_cq_test = src : fi.peak_eq_cq(6, 1000, 4);
// ```
//
// #### 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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// ma = library("maths.lib");
// src = os.osc(440);
// peak_eq_rm_test = src : fi.peak_eq_rm(6, 1000, tan(ma.PI*200/ma.SR));
// ```
//
// #### References
//
// 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 "LicenseRef-STK-4.3";
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)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// spectral_tilt_test = src : fi.spectral_tilt(4, 200, 2000, -0.5);
// ```
//
// #### Example test program
//
// See `dm.spectral_tilt_demo` and the documentation for `no.pink_noise`.
//
// #### References
//
// 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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// levelfilter_test = fi.levelfilter(0.1, 200, src);
// ```
//
// #### References
//
// * <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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// levelfilterN_test = src : fi.levelfilterN(3, 200, 0.1);
// ```
//
// #### References
//
// * <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 "LicenseRef-STK-4.3";
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]`, `(fi.)mth_octave_filterbank`-------------------------
// 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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// mth_octave_filterbank_test = sig : fi.mth_octave_filterbank(3, 2, 8000, 2);
// ```
//------------------------------------------------------------
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 "LicenseRef-STK-4.3";

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 "LicenseRef-STK-4.3";

//------------------------`(fi.)mth_octave_filterbank_alt`-------------------------
//
// #### Usage
//
// ```
// _ : mth_octave_filterbank_alt(O,M,ftop,N) : par(i,N,_)
// ```
//
// Where:
//
// * `O`: filter order used to split each band
// * `M`: number of band-slices per octave
// * `ftop`: highest band-split crossover frequency
// * `N`: total number of bands
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// mth_octave_filterbank_alt_test = sig : fi.mth_octave_filterbank_alt(3, 2, 8000, 2);
// ```
//----------------------------------------------------------
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 "LicenseRef-STK-4.3";

//------------------------`(fi.)mth_octave_filterbank3`-------------------------
//
// #### Usage
//
// ```
// _ : mth_octave_filterbank3(M,ftop,N) : par(i,N,_)
// ```
//
// Where:
//
// * `M`: number of band-slices per octave
// * `ftop`: highest band-split crossover frequency
// * `N`: total number of bands
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// mth_octave_filterbank3_test = sig : fi.mth_octave_filterbank3(2, 8000, 2);
// ```
//----------------------------------------------------------
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 "LicenseRef-STK-4.3";

//------------------------`(fi.)mth_octave_filterbank5`-------------------------
//
// #### Usage
//
// ```
// _ : mth_octave_filterbank5(M,ftop,N) : par(i,N,_)
// ```
//
// Where:
//
// * `M`: number of band-slices per octave
// * `ftop`: highest band-split crossover frequency
// * `N`: total number of bands
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// mth_octave_filterbank5_test = sig : fi.mth_octave_filterbank5(2, 8000, 2);
// ```
//----------------------------------------------------------
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 "LicenseRef-STK-4.3";

//------------------------`(fi.)mth_octave_filterbank_default`-------------------------
// Default Mth-octave filter bank: the 5th-order Butterworth case of
// `fi.mth_octave_filterbank` (`mth_octave_filterbank_default =
// mth_octave_filterbank5`). See `fi.mth_octave_filterbank[n]` for the full
// documentation.
//
// #### Usage
//
// ```
// _ : mth_octave_filterbank_default(M,ftop,N) : par(i,N,_)
// ```
//
// Where:
//
// * `M`: number of band-slices per octave (a constant numerical expression)
// * `ftop`: highest band-split crossover frequency in Hz
// * `N`: total number of bands, including dc and Nyquist (a constant numerical expression)
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// mth_octave_filterbank_default_test = sig : fi.mth_octave_filterbank_default(2, 8000, 2);
// ```
//----------------------------------------------------------
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)) : _,_,_
// ```
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// filterbank_test = src : fi.filterbank(3, (500, 2000));
// ```
//---------------------------------------------------
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 "LicenseRef-STK-4.3";
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)) : _,_,_
// ```
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// filterbanki_test = src : fi.filterbanki(3, (500, 2000));
// ```
//---------------------------------------------------
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 "LicenseRef-STK-4.3";
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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// svf_lp_test = fi.svf.lp(1000, 0.707, sig);
// ```
//---------------------------------------------------
declare svf author "Oleg Nesterov";
declare svf copyright "Copyright (C) 2020 Oleg Nesterov <oleg@redhat.com>";
declare svf license "LicenseRef-STK-4.3";

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].
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// svf_morph_test = fi.svf_morph(1000, 0.707, 1, sig);
// ```
//
// #### 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;
// };
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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].
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// svf_notch_morph_test = fi.svf_notch_morph(1000, 0.707, 1, sig);
// ```
//
// #### 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;
// };
// ```
//
// #### References
//
// <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 "LicenseRef-STK-4.3";
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
// 
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// SVFTPT_SVF_test = fi.SVFTPT.SVF(1000, 0.707, sig);
// ```
//------------------------------------------------------------------------------
declare SVFTPT author "Dario Sanfilippo";
declare SVFTPT copyright
    "Copyright (C) 2024 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare SVFTPT license "MIT";
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.)dynamicSmoothing`------------------------------------------------
//
// 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:
//
// ```
// _ : dynamicSmoothing(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
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// dynamicSmoothing_test = fi.dynamicSmoothing(0.5, 500, sig);
// ```
//
// #### References
//
// * <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";
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.
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// sig = os.osc(440);
// oneEuro_test = sig : fi.oneEuro(1, 0.5, 5);
// ```
//
// #### 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.
//
// #### References
//
// 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`: lowpass cutoff in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// lowpassLR4_test = src : fi.lowpassLR4(1000);
// ```
//------------------------------------------------------------------------------
declare lowpassLR4 author "Dario Sanfilippo";
declare lowpassLR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare lowpassLR4 license "LicenseRef-STK-4.3";
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`: highpass cutoff in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// highpassLR4_test = src : fi.highpassLR4(1000);
// ```
//------------------------------------------------------------------------------
declare highpassLR4 author "Dario Sanfilippo";
declare highpassLR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare highpassLR4 license "LicenseRef-STK-4.3";
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`: crossover split cutoff in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// crossover2LR4_test = src : fi.crossover2LR4(1000);
// ```
//------------------------------------------------------------------------------
declare crossover2LR4 author "Dario Sanfilippo";
declare crossover2LR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare crossover2LR4 license "LicenseRef-STK-4.3";
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`: crossover lower split cutoff in Hz
// * `cf2`: crossover upper split cutoff in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// crossover3LR4_test = src : fi.crossover3LR4(500, 2000);
// ```
//------------------------------------------------------------------------------
declare crossover3LR4 author "Dario Sanfilippo";
declare crossover3LR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare crossover3LR4 license "LicenseRef-STK-4.3";
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`: crossover lower split cutoff in Hz
// * `cf2`: crossover mid split cutoff in Hz
// * `cf3`: crossover upper split cutoff in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// crossover4LR4_test = src : fi.crossover4LR4(300, 1000, 3000);
// ```
//------------------------------------------------------------------------------
declare crossover4LR4 author "Dario Sanfilippo";
declare crossover4LR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare crossover4LR4 license "LicenseRef-STK-4.3";
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`: lowest crossover cutoff frequency in Hz
// * `cf2`: second crossover cutoff frequency in Hz
// * `cf3`: third crossover cutoff frequency in Hz
// * `cf4`: fourth crossover cutoff frequency in Hz
// * `cf5`: fifth crossover cutoff frequency in Hz
// * `cf6`: sixth crossover cutoff frequency in Hz
// * `cf7`: highest crossover cutoff frequency in Hz
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// crossover8LR4_test = src : fi.crossover8LR4(100, 200, 400, 800, 1600, 3200, 6400);
// ```
//------------------------------------------------------------------------------
declare crossover8LR4 author "Dario Sanfilippo";
declare crossover8LR4 copyright
    "Copyright (C) 2022 Dario Sanfilippo <sanfilippo.dario@gmail.com>";
declare crossover8LR4 license "LicenseRef-STK-4.3";
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 : _
// ```
//
// #### Test
// ```
// fi = library("filters.lib");
// os = library("oscillators.lib");
// src = os.osc(440);
// itu_r_bs_1770_4_kfilter_test = src : fi.itu_r_bs_1770_4_kfilter;
// ```
//
// #### References
//
// * <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";

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`: averaging frame duration 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 "LicenseRef-STK-4.3";
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`: time, in seconds, for the system to decay by 1/e,
// or to reach 1-1/e of its final value.
//
// #### References
//
// * <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 "LicenseRef-STK-4.3";
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`: time, in seconds, for the system to decay by 1/1000,
// or to reach 1-1/1000 of its final value.
//
// #### References
//
// * <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 "LicenseRef-STK-4.3";
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`: 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.
//
// #### References
//
// 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 "LicenseRef-STK-4.3";
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";
kalman = kalmanEnv.kalmanFilter;


//============================Adaptive Filters=================================
//=============================================================================
//
// Least-mean-squares adaptive FIR filters: N weights updated by parallel
// gradient recursions, no automatic differentiation involved. The classic
// applications are system identification, echo/noise cancellation and
// adaptive notching.

//-----------`(fi.)lms`, `(fi.)nlms`-------------------------------------------
// Adaptive N-tap FIR filter trained by the (delayed-update) LMS rule:
// the filter output `y` is the dot product of the adaptive weights with the
// last `N` input samples, the error `e = d - y` is fed back, and each weight
// moves along its instantaneous gradient: `w_k <- w_k + mu * e * x@k`.
// `nlms` is the normalized variant: the step is divided by the energy of
// the `N` input samples currently in the filter, making the convergence
// rate independent of the input level (use it unless the input power is
// known and constant).
//
// The weight update uses the one-sample-delayed error (the standard
// "delayed LMS" that a causal sample loop imposes), which does not change
// the converged solution.
//
// #### Usage
//
// ```
// x, d : lms(N, mu) : _, _   // outputs: y (filter output), e (error)
// x, d : nlms(N, mu) : _, _
// ```
//
// Where:
//
// * `N`: number of adaptive weights (a constant numerical expression)
// * `mu`: adaptation step; for `lms`, stability requires
//   `mu < 2/(N*power(x))`, while for `nlms` values in (0, 2) are stable and
//   0.1-1.0 is typical
// * `x`: input signal (the reference fed to the adaptive FIR)
// * `d`: desired signal
//
// A converged identification of an unknown system H feeding
// `d = H(x)` leaves the weights equal to the first `N` samples of the
// impulse response of H, `y` equal to `d`, and `e` near zero.
//
// #### Test
// ```
// fi = library("filters.lib");
// no = library("noises.lib");
// lms_test = no.noise, (no.noise : @(3)*0.5 + no.noise@(1)*0.25) : fi.lms(8, 0.01);
// nlms_test = no.noise, (no.noise : @(3)*0.5 + no.noise@(1)*0.25) : fi.nlms(8, 0.5);
// ```
//
// #### References
//
// * S. Haykin, "Adaptive Filter Theory", Prentice Hall.
// * B. Widrow and S.D. Stearns, "Adaptive Signal Processing", Prentice Hall.
//-----------------------------------------------------------------------------
declare lms license "MIT";
declare nlms license "MIT";
lms(N, mu, x, d) = adaptFIR(N, mu, x, d);
nlms(N, mu, x, d) = adaptFIR(N, mu / (ma.EPSILON + sum(k, N, x@k * x@k)), x, d);

//-----------`(fi.)adaptFIR`---------------------------------------------------
// Generic adaptive N-tap FIR engine that `lms` and `nlms` are built on.
// The filter output `y` is the dot product of the adaptive weights with the
// last `N` input samples, the error `e = d - y` is the single feedback
// signal, and each weight integrates its (one-sample-delayed) gradient
// step: `w_k <- w_k + step * e * x@k`.
//
// Unlike `mu` in `lms`, the `step` parameter is a **signal**, which makes
// this the extension point for the whole family of adaptation rules:
// `nlms` is `adaptFIR` with the step divided by the instantaneous input
// energy, a leaky or gated LMS is `adaptFIR` with a step that decays or is
// forced to 0 while adaptation should be frozen, a scheduled step
// implements search-then-converge, and so on. Use `lms`/`nlms` directly
// unless you need a custom rule.
//
// #### Usage
//
// ```
// x, d : adaptFIR(N, step) : _, _   // outputs: y (filter output), e (error)
// ```
//
// Where:
//
// * `N`: number of adaptive weights (a constant numerical expression)
// * `step`: adaptation step, as a signal evaluated at every sample; the
//   `lms` stability bound `step < 2/(N*power(x))` applies to its
//   instantaneous value
// * `x`: input signal (the reference fed to the adaptive FIR)
// * `d`: desired signal
//
// #### Example
//
// ```
// // NLMS is one normalization away:
// nlms(N, mu, x, d) = adaptFIR(N, mu / (ma.EPSILON + sum(k, N, x@k * x@k)), x, d);
// // gated adaptation: freeze the weights below an input-level threshold
// gatedLMS(N, mu, th, x, d) = adaptFIR(N, mu * (abs(x) > th), x, d);
// ```
//
// #### Test
// ```
// fi = library("filters.lib");
// no = library("noises.lib");
// adaptFIR_test = no.noise, (no.noise : @(3)*0.5 + no.noise@(1)*0.25)
//     : fi.adaptFIR(8, 0.01);
// ```
//
// #### References
//
// * S. Haykin, "Adaptive Filter Theory", Prentice Hall.
//-----------------------------------------------------------------------------
declare adaptFIR license "MIT";
adaptFIR(N, step, x, d) = (loop ~ _) : (!, _, _)
with {
    loop(ep) = e, y, e
    with {
        // w_k[n] = sum_{m<n} step[m] * e[m] * x[m-k]  (delayed LMS update)
        w(k) = (step' * ep * x@(k+1)) : + ~ _;
        y = sum(k, N, w(k) * x@k);
        e = d - y;
    };
};


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

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


// Deprecated backward-compatibility alias, kept for one version and
// scheduled for removal: the canonical name now carries the internal
// underscore prefix (see CONTRIBUTING.md, "Variables and identifiers scoping").
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 "LicenseRef-STK-4.3";
declare pospass0 deprecated "internal core of pospass and pospass6e, renamed _pospass0";
pospass0 = _pospass0;