+++ /dev/null
-// filter.lib - digital filters of various types useful in audio and beyond
-
-declare name "Faust Filter Library";
-declare author "Julius O. Smith (jos at ccrma.stanford.edu)";
-declare copyright "Julius O. Smith III";
-declare version "1.29";
-declare license "STK-4.3"; // Synthesis Tool Kit 4.3 (MIT style license)
-declare reference "https://ccrma.stanford.edu/~jos/filters/";
-
-import("music.lib"); // delay, frac and, from math.lib, SR and PI
-
-//---------------------- zero(z) --------------------------
-// z = location of zero along real axis in z-plane
-// Difference equation: y(n) = x(n) - z * x(n-1)
-// Reference: https://ccrma.stanford.edu/~jos/filters/One_Zero.html
-
-zero(z) = _ <: _,mem : _,*(z) : -;
-
-//------------------------ pole(p) ---------------------------
-// p = pole location = feedback coefficient
-// Could also be called a "leaky integrator".
-// Difference equation: y(n) = x(n) + p * y(n-1)
-// Reference: https://ccrma.stanford.edu/~jos/filters/One_Pole.html
-
-pole(p) = + ~ *(p);
-
-//---------------------- integrator --------------------------
-// pole(1) [implemented separately for block-diagram clarity]
-
-integrator = + ~ _ ;
-
-//----------------------- tau2pole ---------------------------
-// tau2pole(tau) returns a real pole giving exponential decay with
-// tau = time-constant in seconds
-//
-tau2pole(tau) = exp(-1.0/(tau*SR));
-
-//---------------------- smooth(s) --------------------------
-// Exponential smoothing by a unity-dc-gain one-pole lowpass
-//
-// USAGE: smooth(tau2pole(tau)), where
-// tau = desired smoothing time constant in seconds,
-// or
-// smooth(s), where s = smoothness between 0 and 1.
-// s=0 for no smoothing
-// s=0.999 is "very smooth"
-// s>1 is unstable, and s=1 yields the zero signal for all inputs.
-// The exponential time-constant is approximately
-// 1/(1-s) samples, when s is close to (but less than) 1.
-// Reference:
-// https://ccrma.stanford.edu/~jos/mdft/Convolution_Example_2_ADSR.html
-
-smooth(s) = *(1.0 - s) : + ~ *(s);
-
-//------------------- dcblockerat(fb) -----------------------
-// fb = "break frequency" in Hz, i.e., -3 dB gain frequency.
-// The amplitude response is substantially flat above fb,
-// and sloped at about +6 dB/octave below fb.
-// Derived from the analog transfer function
-// H(s) = s / (s + 2*PI*fb)
-// by the low-frequency-matching bilinear transform method
-// (i.e., the standard frequency-scaling constant 2*SR).
-// Reference:
-// https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html
-
-dcblockerat(fb) = *(b0) : zero(1) : pole(p)
-with {
- wn = PI*fb/SR;
- b0 = 1.0 / (1 + wn);
- p = (1 - wn) * b0;
-};
-
-//---------------------- dcblocker --------------------------
-// Default dc blocker has -3dB point near 35 Hz (at 44.1 kHz)
-// and high-frequency gain near 1.0025 (due to no scaling)
-//
-dcblocker = zero(1) : pole(0.995);
-
-//------------ notchw(width,freq), notch(freq) --------------
-// width = "notch width" in Hz (approximate)
-// freq = "notch frequency" in Hz
-// Reference:
-// https://ccrma.stanford.edu/~jos/pasp/Phasing_2nd_Order_Allpass_Filters.html
-
-notchw(width,freq) = tf2(b0,b1,b2,a1,a2)
-with {
- fb = 0.5*width; // First design a dcblockerat(width/2)
- wn = PI*fb/SR;
- b0db = 1.0 / (1 + wn);
- p = (1 - wn) * b0db; // This is our pole radius.
- // Now place unit-circle zeros at desired angles:
- tn = 2*PI*freq/SR;
- a2 = p * p;
- a2p1 = 1+a2;
- a1 = -a2p1*cos(tn);
- b1 = a1;
- b0 = 0.5*a2p1;
- b2 = b0;
-};
-
-//========================= Comb Filters ===============================
-
-//----------------------- ff_comb, ff_fcomb ----------------------------
-// Feed-Forward Comb Filter
-//
-// USAGE:
-// _ : ff_comb(maxdel,intdel,b0,bM) : _
-// _ : ff_fcomb(maxdel,del,b0,bM) : _
-// where
-// maxdel = maximum delay (a power of 2)
-// intdel = current (integer) comb-filter delay between 0 and maxdel
-// del = current (float) comb-filter delay between 0 and maxdel
-// b0 = gain applied to delay-line input
-// bM = gain applied to delay-line output and then summed with input
-//
-// Note that ff_comb requires integer delays (uses delay() internally)
-// while ff_fcomb takes floating-point delays (uses fdelay() internally).
-//
-// REFERENCE:
-// https://ccrma.stanford.edu/~jos/pasp/Feedforward_Comb_Filters.html
-
-ff_comb (maxdel,M,b0,bM) = _ <: *(b0), bM * delay(maxdel,M) : + ;
-ff_fcomb(maxdel,M,b0,bM) = _ <: *(b0), bM * fdelay(maxdel,M) : + ;
-
-// Typical special case:
-ffcombfilter(maxdel,del,g) = ff_comb(maxdel,del,1,g);
-
-//----------------------- fb_comb, fb_fcomb, rev1 -----------------------
-// Feed-Back Comb Filter
-//
-// USAGE:
-// _ : fb_comb(maxdel,intdel,b0,aN) : _
-// _ : fb_fcomb(maxdel,del,b0,aN) : _
-// _ : rev1(maxdel,del,-aN) : _
-// where
-// maxdel = maximum delay (a power of 2)
-// intdel = current (integer) comb-filter delay between 0 and maxdel
-// del = current (float) comb-filter delay between 0 and maxdel
-// b0 = gain applied to delay-line input and forwarded to output
-// aN = minus the gain applied to delay-line output before
-// summing with the input and feeding to the delay line
-//
-// Reference:
-// https://ccrma.stanford.edu/~jos/pasp/Feedback_Comb_Filters.html
-
-fb_comb (maxdel,N,b0,aN) = (+ <: delay(maxdel,N),_) ~ *(-aN) : !,*(b0);
-fb_fcomb(maxdel,N,b0,aN) = (+ <: fdelay(maxdel,N),_) ~ *(-aN) : !,*(b0);
-
-// The "rev1 section" dates back to the 1960s in computer-music reverberation.
-// See the jcrev and brassrev in effect.lib for usage examples.
-rev1(maxdel,N,g) = fb_comb (maxdel,N,1,-g);
-
-// Typical special case:
-fbcombfilter(maxdel,intdel,g) = (+ : delay(maxdel,intdel)) ~ *(g);
-ffbcombfilter(maxdel,del,g) = (+ : fdelay(maxdel,del)) ~ *(g);
-
-//------------------- allpass_comb, allpass_fcomb, rev2 -----------------
-// Schroeder Allpass Comb Filter
-//
-// USAGE:
-// _ : allpass_comb (maxdel,intdel,aN) : _
-// _ : allpass_fcomb(maxdel,del,aN) : _
-// _ : rev2(maxdel,del,-aN) : _
-// where
-// maxdel = maximum delay (a power of 2)
-// intdel = current (integer) comb-filter delay between 0 and maxdel
-// del = current (float) comb-filter delay between 0 and maxdel
-// aN = minus the feedback gain
-//
-// Note that allpass_comb(maxlen,len,aN) =
-// ff_comb(maxlen,len,aN,1) :
-// fb_comb(maxlen,len-1,1,aN);
-// which is a direct-form-1 implementation, requiring two delay lines.
-// The implementation here is direct-form-2 requiring only one delay line.
-//
-// REFERENCES:
-// https://ccrma.stanford.edu/~jos/pasp/Allpass_Two_Combs.html
-// https://ccrma.stanford.edu/~jos/pasp/Schroeder_Allpass_Sections.html
-// https://ccrma.stanford.edu/~jos/filters/Four_Direct_Forms.html
-
-allpass_comb(maxdel,N,aN) = (+ <:
- delay(maxdel,N-1),*(aN)) ~ *(-aN)
- : mem,_ : + ;
-
-// The "rev2 section" dates back to the 1960s in computer-music reverberation:
-rev2(maxlen,len,g) = allpass_comb(maxlen,len,-g);
-
-//================ Direct-Form Digital Filter Sections ================
-
-// Specified by transfer-function polynomials B(z)/A(z) as in matlab
-
-//---------------------------- iir (tfN) -------------------------------
-// Nth-order Infinite-Impulse-Response (IIR) digital filter,
-// implemented in terms of the Transfer-Function (TF) coefficients.
-// Such filter structures are termed "direct form".
-//
-// USAGE:
-// _ : iir(bcoeffs,acoeffs) : _
-// where
-// order = filter order (int) = max(#poles,#zeros)
-// bcoeffs = (b0,b1,...,b_order) = TF numerator coefficients
-// acoeffs = (a1,...,a_order) = TF denominator coeffs (a0=1)
-//
-// REFERENCE:
-// https://ccrma.stanford.edu/~jos/filters/Four_Direct_Forms.html
-
-iir(bv,av) = sub ~ fir(av) : fir(bv);
-
-//----------------------------- sub ---------------------------------
-sub(x,y) = y-x; // move to math.lib?
-
-//----------------------------- fir ---------------------------------
-// FIR filter (convolution of FIR filter coefficients with a signal)
-//
-// USAGE:
-// _ : fir(bv) : _
-// where bv = b0,b1,...,bn is a parallel bank of coefficient signals.
-// NOTE: bv is processed using pattern-matching at compile time,
-// so it must have this normal form (parallel signals).
-// EXAMPLE: Smoothing white noise with a five-point moving average:
-// bv = .2,.2,.2,.2,.2;
-// process = noise : fir(bv);
-// EQUIVALENT (note double parens):
-// process = noise : fir((.2,.2,.2,.2,.2));
-
-fir(bv) = conv(bv);
-
-//--------------------------- conv, convN -------------------------------
-// Convolution of input signal with given coefficients
-//
-// USAGE:
-// _ : conv((k1,k2,k3,...,kN)) : _; // Argument = one signal bank
-// _ : convN(N,(k1,k2,k3,...)) : _; // Useful when N < count((k1,...))
-
-convN(N,kv,x) = sum(i,N,take(i+1,kv) * x@i); // take() defined in math.lib
-
-conv(kv,x) = sum(i,count(kv),take(i+1,kv) * x@i); // count() from math.lib
-
-// Named special cases:
-//----------------------------- tf1, tf2 ---------------------------------
-// tfN = N'th-order direct-form digital filter
-tf1(b0,b1,a1) = _ <: *(b0), (mem : *(b1)) :> + ~ *(0-a1);
-tf2(b0,b1,b2,a1,a2) = iir((b0,b1,b2),(a1,a2)); // cf. TF2 in music.lib)
-
-//===================== Ladder/Lattice Digital Filters ======================
-// Ladder and lattice digital filters generally have superior numerical
-// properties relative to direct-form digital filters. They can be derived
-// from digital waveguide filters, which gives them a physical interpretation.
-
-// REFERENCES:
-// F. Itakura and S. Saito: "Digital Filtering Techniques for Speech Analysis and Synthesis",
-// 7th Int. Cong. Acoustics, Budapest, 25 C 1, 1971.
-// J. D. Markel and A. H. Gray: Linear Prediction of Speech, New York: Springer Verlag, 1976.
-// https://ccrma.stanford.edu/~jos/pasp/Conventional_Ladder_Filters.html
-
-//------------------------------ block, crossn,crossn1 ----------------------------------
-// signal block/crossing utilities
-// (move to math.lib?)
-
-// block - terminate n signals (goes with bus(n) in math.lib)
-block(n) = par(i,n,!);
-
-// crossnn - cross two bus(n)s:
-crossnn(n) = bus(n),bus(n) <: block(n),bus(n),bus(n),block(n);
-
-// crossn1 - cross bus(n) and bus(1):
-crossn1(n) = bus(n),(bus(1)<:bus(n)) <: block(n),bus(n),bus(n),block(n):bus(1),block(n-1),bus(n);
-
-//------------------------------- av2sv -----------------------------------
-// Compute reflection coefficients sv from transfer-function denominator av
-//
-// USAGE:
-// sv = av2sv(av)
-// where
-// av = parallel signal bank a1,...,aN
-// sv = parallel signal bank s1,...,sN
-// where si = ith reflection coefficient, and
-// ai = coefficient of z^(-i) in the filter
-// transfer-function denominator A(z).
-//
-// REFERENCE:
-// https://ccrma.stanford.edu/~jos/filters/Step_Down_Procedure.html
-// (where reflection coefficients are denoted by k rather than s).
-
-av2sv(av) = par(i,M,s(i+1)) with {
- M = count(av);
- s(m) = sr(M-m+1); // m=1..M
- sr(m) = Ari(m,M-m+1); // s_{M-1-m}
- Ari(m,i) = take(i+1,Ar(m-1));
- //step-down recursion for lattice/ladder digital filters:
- Ar(0) = (1,av); // Ar(m) is order M-m (i.e. "reverse-indexed")
- Ar(m) = 1,par(i,M-m, (Ari(m,i+1) - sr(m)*Ari(m,M-m-i))/(1-sr(m)*sr(m)));
-};
-
-//---------------------------- bvav2nuv --------------------------------
-// Compute lattice tap coefficients from transfer-function coefficients
-//
-// USAGE:
-// nuv = bvav2nuv(bv,av)
-// where
-// av = parallel signal bank a1,...,aN
-// bv = parallel signal bank b0,b1,...,aN
-// nuv = parallel signal bank nu1,...,nuN
-// where nui is the i'th tap coefficient,
-// bi is the coefficient of z^(-i) in the filter numerator,
-// ai is the coefficient of z^(-i) in the filter denominator
-
-bvav2nuv(bv,av) = par(m,M+1,nu(m)) with {
- M = count(av);
- nu(m) = take(m+1,Pr(M-m)); // m=0..M
- // lattice/ladder tap parameters:
- Pr(0) = bv; // Pr(m) is order M-m, 'r' means "reversed"
- Pr(m) = par(i,M-m+1, (Pri(m,i) - nu(M-m+1)*Ari(m,M-m-i+1)));
- Pri(m,i) = take(i+1,Pr(m-1));
- Ari(m,i) = take(i+1,Ar(m-1));
- //step-down recursion for lattice/ladder digital filters:
- Ar(0) = (1,av); // Ar(m) is order M-m (recursion index must start at constant)
- Ar(m) = 1,par(i,M-m, (Ari(m,i+1) - sr(m)*Ari(m,M-m-i))/(1-sr(m)*sr(m)));
- sr(m) = Ari(m,M-m+1); // s_{M-1-m}
-};
-
-//--------------------------- iir_lat2, allpassnt -----------------------
-
-iir_lat2(bv,av) = allpassnt(M,sv) : sum(i,M+1,*(take(M-i+1,tg)))
-with {
- M = count(av);
- sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
- tg = bvav2nuv(bv,av); // tg = vector of tap gains
-};
-
-// two-multiply lattice allpass (nested order-1 direct-form-ii allpasses):
-allpassnt(0,sv) = _;
-allpassnt(n,sv) =
-//0: x <: ((+ <: (allpassnt(n-1,sv)),*(s))~(*(-s))) : _',_ :+
- _ : ((+ <: (allpassnt(n-1,sv),*(s)))~*(-s)) : fsec(n)
-with {
- fsec(1) = crossnn(1) : _, (_<:mem,_) : +,_;
- fsec(n) = crossn1(n) : _, (_<:mem,_),par(i,n-1,_) : +, par(i,n,_);
- innertaps(n) = par(i,n,_);
- s = take(n,sv); // reflection coefficient s = sin(theta)
-};
-
-//------------------------------- iir_kl, allpassnklt -------------------------
-iir_kl(bv,av) = allpassnklt(M,sv) : sum(i,M+1,*(tghr(i)))
-with {
- M = count(av);
- sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
- tg = bvav2nuv(bv,av); // tg = vector of tap gains for 2mul case
- tgr(i) = take(M+1-i,tg);
- tghr(n) = tgr(n)/pi(n);
- pi(0) = 1;
- pi(n) = pi(n-1)*(1+take(M-n+1,sv)); // all sign parameters '+'
-};
-
-// Kelly-Lochbaum ladder allpass with tap lines:
-allpassnklt(0,sv) = _;
-allpassnklt(n,sv) = _ <: *(s),(*(1+s) : (+
- : allpassnklt(n-1,sv))~(*(-s))) : fsec(n)
-with {
- fsec(1) = _, (_<:mem*(1-s),_) : sumandtaps(n);
- fsec(n) = _, (_<:mem*(1-s),_), par(i,n-1,_) : sumandtaps(n);
- s = take(n,sv);
- sumandtaps(n) = +,par(i,n,_);
-};
-
-
-//------------------------------- iir_lat1, allpassn1mt -------------------------
-iir_lat1(bv,av) = allpassn1mt(M,sv) : sum(i,M+1,*(tghr(i+1)))
-with {
- M = count(av);
- sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
- tg = bvav2nuv(bv,av); // tg = vector of tap gains
- tgr(i) = take(M+2-i,tg); // i=1..M+1 (for "takability")
- tghr(n)=tgr(n)/pi(n);
- pi(1) = 1;
- pi(n) = pi(n-1)*(1+take(M-n+2,sv)); // all sign parameters '+'
-};
-
-// one-multiply lattice allpass with tap lines:
-allpassn1mt(0,sv) = _;
-allpassn1mt(n,sv)= _ <: _,_ : ((+:*(s) <: _,_),_ : _,+ : crossnn(1)
- : allpassn1mt(n-1,sv),_)~(*(-1)) : fsec(n)
-with {
-//0: fsec(n) = _',_ : +
- fsec(1) = crossnn(1) : _, (_<:mem,_) : +,_;
- fsec(n) = crossn1(n) : _, (_<:mem,_),par(i,n-1,_) : +, par(i,n,_);
- innertaps(n) = par(i,n,_);
- s = take(n,sv); // reflection coefficient s = sin(theta)
-};
-
-//------------------------------- iir_nl, allpassnnlt -------------------------
-// Normalized ladder filter
-//
-// REFERENCES:
-// J. D. Markel and A. H. Gray, Linear Prediction of Speech, New York: Springer Verlag, 1976.
-// https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html
-
-iir_nl(bv,av) = allpassnnlt(M,sv) : sum(i,M+1,*(tghr(i)))
-with {
- M = count(av);
- sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
- tg = bvav2nuv(bv,av); // tg = vector of tap gains for 2mul case
- tgr(i) = take(M+1-i,tg);
- tghr(n) = tgr(n)/pi(n);
- pi(0) = 1;
- s(n) = take(M-n+1,sv);
- c(n) = sqrt(max(0,1-s(n)*s(n))); // compiler crashes on sqrt(-)
- pi(n) = pi(n-1)*c(n);
-};
-
-// Normalized ladder allpass with tap lines:
-allpassnnlt(0,sv) = _;
-allpassnnlt(n,scl*(sv)) = allpassnnlt(n,par(i,count(sv),scl*(sv(i))));
-allpassnnlt(n,sv) = _ <: *(s),(*(c) : (+
- : allpassnnlt(n-1,sv))~(*(-s))) : fsec(n)
-with {
- fsec(1) = _, (_<:mem*(c),_) : sumandtaps(n);
- fsec(n) = _, (_<:mem*(c),_), par(i,n-1,_) : sumandtaps(n);
- s = take(n,sv);
- c = sqrt(max(0,1-s*s));
- sumandtaps(n) = +,par(i,n,_);
-};
-
-//========================= Useful special cases ============================
-
-//-------------------------------- tf2np ------------------------------------
-// tf2np - biquad based on a stable second-order Normalized Ladder Filter
-// (more robust to modulation than tf2 and protected against instability)
-tf2np(b0,b1,b2,a1,a2) = allpassnnlt(M,sv) : sum(i,M+1,*(tghr(i)))
-with {
- smax = 0.9999; // maximum reflection-coefficient magnitude allowed
- s2 = max(-smax, min(smax,a2)); // Project both reflection-coefficients
- s1 = max(-smax, min(smax,a1/(1+a2))); // into the defined stability-region.
- sv = (s1,s2); // vector of sin(theta) reflection coefficients
- M = 2;
- nu(2) = b2;
- nu(1) = b1 - b2*a1;
- nu(0) = (b0-b2*a2) - nu(1)*s1;
- tg = (nu(0),nu(1),nu(2));
- tgr(i) = take(M+1-i,tg); // vector of tap gains for 2mul case
- tghr(n) = tgr(n)/pi(n); // apply pi parameters for NLF case
- pi(0) = 1;
- s(n) = take(M-n+1,sv);
- c(n) = sqrt(1-s(n)*s(n));
- pi(n) = pi(n-1)*c(n);
-};
-
-//----------------------------- wgr ---------------------------------
-// Second-order transformer-normalized digital waveguide resonator
-// USAGE:
-// _ : wgr(f,r) : _
-// where
-// f = resonance frequency (Hz)
-// r = loss factor for exponential decay
-// (set to 1 to make a numerically stable oscillator)
-//
-// REFERENCES:
-// https://ccrma.stanford.edu/~jos/pasp/Power_Normalized_Waveguide_Filters.html
-// https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html
-//
-wgr(f,r,x) = (*(G),_<:_,((+:*(C))<:_,_),_:+,_,_:+(x),-) ~ cross : _,*(0-gi)
-with {
- C = cos(2*PI*f/SR);
- gi = sqrt(max(0,(1+C)/(1-C))); // compensate amplitude (only needed when
- G = r*(1-1' + gi')/gi; // frequency changes substantially)
- cross = _,_ <: !,_,_,!;
-};
-
-//----------------------------- nlf2 --------------------------------
-// Second order normalized digital waveguide resonator
-// USAGE:
-// _ : nlf2(f,r) : _
-// where
-// f = resonance frequency (Hz)
-// r = loss factor for exponential decay
-// (set to 1 to make a sinusoidal oscillator)
-//
-// REFERENCE:
-// https://ccrma.stanford.edu/~jos/pasp/Power_Normalized_Waveguide_Filters.html
-//
-nlf2(f,r,x) = ((_<:_,_),(_<:_,_) : (*(s),*(c),*(c),*(0-s)) :>
- (*(r),+(x))) ~ cross
-with {
- th = 2*PI*f/SR;
- c = cos(th);
- s = sin(th);
- cross = _,_ <: !,_,_,!;
-};
-
-//===================== Ladder/Lattice Allpass Filters ======================
-// An allpass filter has gain 1 at every frequency, but variable phase.
-// Ladder/lattice allpass filters are specified by reflection coefficients.
-// They are defined here as nested allpass filters, hence the names allpassn*.
-//
-// REFERENCES
-// 1. https://ccrma.stanford.edu/~jos/pasp/Conventional_Ladder_Filters.html
-// https://ccrma.stanford.edu/~jos/pasp/Nested_Allpass_Filters.html
-// 2. Linear Prediction of Speech, Markel and Gray, Springer Verlag, 1976
-//
-// QUICK GUIDE
-// allpassn - two-multiply lattice - each section is two multiply-adds
-// allpassnn - normalized form - four multiplies and two adds per section,
-// but coefficients can be time varying and nonlinear without
-// "parametric amplification" (modulation of signal energy).
-// allpassnkl - Kelly-Lochbaum form - four multiplies and two adds per
-// section, but all signals have an immediate physical
-// interpretation as traveling pressure waves, etc.
-// allpassn1m - One-multiply form - one multiply and three adds per section.
-// Normally the most efficient in special-purpose hardware.
-//
-// TYPICAL USAGE
-// _ : allpassn(N,sv) : _
-// where
-// N = allpass order (number of ladder or lattice sections)
-// sv = (s1,s2,...,sN) = reflection coefficients (between -1 and 1).
-// For allpassnn only, sv is replaced by tv, where sv(i) = sin(tv(i)),
-// where tv(i) may range between -PI and PI.
-//
-// two-multiply:
-allpassn(0,sv) = _;
-allpassn(n,sv) = _ <: ((+ <: (allpassn(n-1,sv)),*(s))~(*(-s))) : _',_ :+
-with { s = take(n,sv); };
-
-// power-normalized (reflection coefficients s = sin(t)):
-allpassnn(0,tv) = _;
-allpassnn(n,tv) = _ <: *(s), (*(c) : (+
- : allpassnn(n-1,tv))~(*(-s))) : _, mem*c : +
-with { c=cos(take(n,tv)); s=sin(take(n,tv)); };
-
-// power-normalized with sparse delays dv(n)>1:
-allpassnns(0,tv,dmax,dv) = _;
-allpassnns(n,tv,dmax,dv) = _ <: *(s), (*(c) : (+ : dl
- : allpassnns(n-1,tv,dmax,dv))~(*(-s))) : _, mem*c : +
-with { c=cos(take(n,tv)); s=sin(take(n,tv));
- dl=delay(dmax,(take(n,dv)-1)); };
-
-// Kelly-Lochbaum:
-allpassnkl(0,sv) = _;
-allpassnkl(n,sv) = _ <: *(s),(*(1+s) : (+
- : allpassnkl(n-1,sv))~(*(-s))) : _, mem*(1-s) : +
-with { s = take(n,sv); };
-
-// one-multiply:
-allpassn1m(0,sv) = _;
-allpassn1m(n,sv)= _ <: _,_ : ((+:*(s) <: _,_),_ : _,+ : cross
- : allpassn1m(n-1,sv),_)~(*(-1)) : _',_ : +
-with {s = take(n,sv); cross = _,_ <: !,_,_,!; };
-
-//===== Digital Filter Sections Specified as Analog Filter Sections =====
-//
-//------------------------- tf2s, tf2snp --------------------------------
-// Second-order direct-form digital filter,
-// specified by ANALOG transfer-function polynomials B(s)/A(s),
-// and a frequency-scaling parameter. Digitization via the
-// bilinear transform is built in.
-//
-// USAGE: tf2s(b2,b1,b0,a1,a0,w1), where
-//
-// b2 s^2 + b1 s + b0
-// H(s) = --------------------
-// s^2 + a1 s + a0
-//
-// and w1 is the desired digital frequency (in radians/second)
-// corresponding to analog frequency 1 rad/sec (i.e., s = j).
-//
-// EXAMPLE: A second-order ANALOG Butterworth lowpass filter,
-// normalized to have cutoff frequency at 1 rad/sec,
-// has transfer function
-//
-// 1
-// H(s) = -----------------
-// s^2 + a1 s + 1
-//
-// where a1 = sqrt(2). Therefore, a DIGITAL Butterworth lowpass
-// cutting off at SR/4 is specified as tf2s(0,0,1,sqrt(2),1,PI*SR/2);
-//
-// METHOD: Bilinear transform scaled for exact mapping of w1.
-// REFERENCE:
-// https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html
-//
-tf2s(b2,b1,b0,a1,a0,w1) = tf2(b0d,b1d,b2d,a1d,a2d)
-with {
- c = 1/tan(w1*0.5/SR); // bilinear-transform scale-factor
- csq = c*c;
- d = a0 + a1 * c + csq;
- b0d = (b0 + b1 * c + b2 * csq)/d;
- b1d = 2 * (b0 - b2 * csq)/d;
- b2d = (b0 - b1 * c + b2 * csq)/d;
- a1d = 2 * (a0 - csq)/d;
- a2d = (a0 - a1*c + csq)/d;
-};
-
-// tf2snp = tf2s but using a protected normalized ladder filter for tf2:
-tf2snp(b2,b1,b0,a1,a0,w1) = tf2np(b0d,b1d,b2d,a1d,a2d)
-with {
- c = 1/tan(w1*0.5/SR); // bilinear-transform scale-factor
- csq = c*c;
- d = a0 + a1 * c + csq;
- b0d = (b0 + b1 * c + b2 * csq)/d;
- b1d = 2 * (b0 - b2 * csq)/d;
- b2d = (b0 - b1 * c + b2 * csq)/d;
- a1d = 2 * (a0 - csq)/d;
- a2d = (a0 - a1*c + csq)/d;
-};
-
-//----------------------------- tf1s --------------------------------
-// First-order direct-form digital filter,
-// specified by ANALOG transfer-function polynomials B(s)/A(s),
-// and a frequency-scaling parameter.
-//
-// USAGE: tf1s(b1,b0,a0,w1), where
-//
-// b1 s + b0
-// H(s) = ----------
-// s + a0
-//
-// and w1 is the desired digital frequency (in radians/second)
-// corresponding to analog frequency 1 rad/sec (i.e., s = j).
-//
-// EXAMPLE: A first-order ANALOG Butterworth lowpass filter,
-// normalized to have cutoff frequency at 1 rad/sec,
-// has transfer function
-//
-// 1
-// H(s) = -------
-// s + 1
-//
-// so b0 = a0 = 1 and b1 = 0. Therefore, a DIGITAL first-order
-// Butterworth lowpass with gain -3dB at SR/4 is specified as
-//
-// tf1s(0,1,1,PI*SR/2); // digital half-band order 1 Butterworth
-//
-// METHOD: Bilinear transform scaled for exact mapping of w1.
-// REFERENCE:
-// https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html
-//
-tf1s(b1,b0,a0,w1) = tf1(b0d,b1d,a1d)
-with {
- c = 1/tan((w1)*0.5/SR); // bilinear-transform scale-factor
- d = a0 + c;
- b1d = (b0 - b1*c) / d;
- b0d = (b0 + b1*c) / d;
- a1d = (a0 - c) / d;
-};
-
-//----------------------------- tf2sb --------------------------------
-// Bandpass mapping of tf2s: In addition to a frequency-scaling parameter
-// w1 (set to HALF the desired passband width in rad/sec),
-// there is a desired center-frequency parameter wc (also in rad/s).
-// Thus, tf2sb implements a fourth-order digital bandpass filter section
-// specified by the coefficients of a second-order analog lowpass prototpe
-// section. Such sections can be combined in series for higher orders.
-// The order of mappings is (1) frequency scaling (to set lowpass cutoff w1),
-// (2) bandpass mapping to wc, then (3) the bilinear transform, with the
-// usual scale parameter 2*SR. Algebra carried out in maxima and pasted here.
-//
-tf2sb(b2,b1,b0,a1,a0,w1,wc) =
- iir((b0d/a0d,b1d/a0d,b2d/a0d,b3d/a0d,b4d/a0d),(a1d/a0d,a2d/a0d,a3d/a0d,a4d/a0d)) with {
- T = 1.0/float(SR);
- b0d = (4*b0*w1^2+8*b2*wc^2)*T^2+8*b1*w1*T+16*b2;
- b1d = 4*b2*wc^4*T^4+4*b1*wc^2*w1*T^3-16*b1*w1*T-64*b2;
- b2d = 6*b2*wc^4*T^4+(-8*b0*w1^2-16*b2*wc^2)*T^2+96*b2;
- b3d = 4*b2*wc^4*T^4-4*b1*wc^2*w1*T^3+16*b1*w1*T-64*b2;
- b4d = (b2*wc^4*T^4-2*b1*wc^2*w1*T^3+(4*b0*w1^2+8*b2*wc^2)*T^2-8*b1*w1*T +16*b2)
- + b2*wc^4*T^4+2*b1*wc^2*w1*T^3;
- a0d = wc^4*T^4+2*a1*wc^2*w1*T^3+(4*a0*w1^2+8*wc^2)*T^2+8*a1*w1*T+16;
- a1d = 4*wc^4*T^4+4*a1*wc^2*w1*T^3-16*a1*w1*T-64;
- a2d = 6*wc^4*T^4+(-8*a0*w1^2-16*wc^2)*T^2+96;
- a3d = 4*wc^4*T^4-4*a1*wc^2*w1*T^3+16*a1*w1*T-64;
- a4d = wc^4*T^4-2*a1*wc^2*w1*T^3+(4*a0*w1^2+8*wc^2)*T^2-8*a1*w1*T+16;
-};
-
-//----------------------------- tf1sb --------------------------------
-// First-to-second-order lowpass-to-bandpass section mapping,
-// analogous to tf2sb above.
-//
-tf1sb(b1,b0,a0,w1,wc) = tf2(b0d/a0d,b1d/a0d,b2d/a0d,a1d/a0d,a2d/a0d) with {
- T = 1.0/float(SR);
- a0d = wc^2*T^2+2*a0*w1*T+4;
- b0d = b1*wc^2*T^2 +2*b0*w1*T+4*b1;
- b1d = 2*b1*wc^2*T^2-8*b1;
- b2d = b1*wc^2*T^2-2*b0*w1*T+4*b1;
- a1d = 2*wc^2*T^2-8;
- a2d = wc^2*T^2-2*a0*w1*T+4;
-};
-
-//====================== Simple Resonator Filters ======================
-
-// resonlp = 2nd-order lowpass with corner resonance:
-resonlp(fc,Q,gain) = tf2s(b2,b1,b0,a1,a0,wc)
-with {
- wc = 2*PI*fc;
- a1 = 2/Q;
- a0 = 1;
- b2 = 0;
- b1 = 0;
- b0 = gain;
-};
-
-// resonhp = 2nd-order highpass with corner resonance:
-resonhp(fc,Q,gain,x) = gain*x-resonlp(fc,Q,gain,x);
-
-// resonbp = 2nd-order bandpass
-resonbp(fc,Q,gain) = tf2s(b2,b1,b0,a1,a0,wc)
-with {
- wc = 2*PI*fc;
- a1 = 2/Q;
- a0 = 1;
- b2 = 0;
- b1 = gain;
- b0 = 0;
-};
-
-//================ Butterworth Lowpass/Highpass Filters ======================
-// Nth-order Butterworth lowpass or highpass filters
-//
-// USAGE:
-// _ : lowpass(N,fc) : _
-// _ : highpass(N,fc) : _
-// where
-// N = filter order (number of poles) [nonnegative integer]
-// fc = desired cut-off frequency (-3dB frequency) in Hz
-// REFERENCE:
-// https://ccrma.stanford.edu/~jos/filters/Butterworth_Lowpass_Design.html
-// 'butter' function in Octave ("[z,p,g] = butter(N,1,'s');")
-// ACKNOWLEDGMENT
-// Generalized recursive formulation initiated by Yann Orlarey.
-
-lowpass(N,fc) = lowpass0_highpass1(0,N,fc);
-highpass(N,fc) = lowpass0_highpass1(1,N,fc);
-lowpass0_highpass1(s,N,fc) = lphpr(s,N,N,fc)
-with {
- lphpr(s,0,N,fc) = _;
- lphpr(s,1,N,fc) = tf1s(s,1-s,1,2*PI*fc);
- lphpr(s,O,N,fc) = lphpr(s,(O-2),N,fc) : tf2s(s,0,1-s,a1s,1,w1) with {
- parity = N % 2;
- S = (O-parity)/2; // current section number
- a1s = -2*cos(-PI + (1-parity)*PI/(2*N) + (S-1+parity)*PI/N);
- w1 = 2*PI*fc;
- };
-};
-
-//========== Special Filter-Bank Delay-Equalizing Allpass Filters ===========
-//
-// These special allpass filters are needed by filterbank et al. below.
-// They are equivalent to (lowpass(N,fc) +|- highpass(N,fc))/2, but with
-// canceling pole-zero pairs removed (which occurs for odd N).
-
-//-------------------- lowpass_plus|minus_highpass ------------------
-
-highpass_plus_lowpass(1,fc) = _;
-highpass_plus_lowpass(3,fc) = tf2s(1,-1,1,1,1,w1) with { w1 = 2*PI*fc; };
-highpass_minus_lowpass(3,fc) = tf1s(-1,1,1,w1) with { w1 = 2*PI*fc; };
-highpass_plus_lowpass(5,fc) = tf2s(1,-a11,1,a11,1,w1)
-with {
- a11 = 1.618033988749895;
- w1 = 2*PI*fc;
-};
-highpass_minus_lowpass(5,fc) = tf1s(1,-1,1,w1) : tf2s(1,-a12,1,a12,1,w1)
-with {
- a12 = 0.618033988749895;
- w1 = 2*PI*fc;
-};
-
-// Catch-all definitions for generality - even order is done:
-
-highpass_plus_lowpass(N,fc) = switch_odd_even(N%2,N,fc) with {
- switch_odd_even(0,N,fc) = highpass_plus_lowpass_even(N,fc);
- switch_odd_even(1,N,fc) = highpass_plus_lowpass_odd(N,fc);
-};
-
-highpass_minus_lowpass(N,fc) = switch_odd_even(N%2,N,fc) with {
- switch_odd_even(0,N,fc) = highpass_minus_lowpass_even(N,fc);
- switch_odd_even(1,N,fc) = highpass_minus_lowpass_odd(N,fc);
-};
-
-highpass_plus_lowpass_even(N,fc) = highpass(N,fc) + lowpass(N,fc);
-highpass_minus_lowpass_even(N,fc) = highpass(N,fc) - lowpass(N,fc);
-
-// FIXME: Rewrite the following, as for orders 3 and 5 above,
-// to eliminate pole-zero cancellations:
-highpass_plus_lowpass_odd(N,fc) = highpass(N,fc) + lowpass(N,fc);
-highpass_minus_lowpass_odd(N,fc) = highpass(N,fc) - lowpass(N,fc);
-
-//===================== Elliptic (Cauer) Lowpass Filters ===================
-// USAGE:
-// _ : lowpass3e(fc) : _
-// _ : lowpass6e(fc) : _
-// where fc = -3dB frequency in Hz
-//
-// REFERENCES:
-// http://en.wikipedia.org/wiki/Elliptic_filter
-// functions 'ncauer' and 'ellip' in Octave
-
-//----------------------------- lowpass3e -----------------------------
-// Third-order Elliptic (Cauer) lowpass filter
-// DESIGN: For spectral band-slice level display (see octave_analyzer3e):
-// [z,p,g] = ncauer(Rp,Rs,3); % analog zeros, poles, and gain, where
-// Rp = 60 % dB ripple in stopband
-// Rs = 0.2 % dB ripple in passband
-//
-lowpass3e(fc) = tf2s(b21,b11,b01,a11,a01,w1) : tf1s(0,1,a02,w1)
-with {
- a11 = 0.802636764161030; // format long; poly(p(1:2)) % in octave
- a01 = 1.412270893774204;
- a02 = 0.822445908998816; // poly(p(3)) % in octave
- b21 = 0.019809144837789; // poly(z)
- b11 = 0;
- b01 = 1.161516418982696;
- w1 = 2*PI*fc;
-};
-
-//----------------------------- lowpass6e -----------------------------
-// Sixth-order Elliptic/Cauer lowpass filter
-// DESIGN: For spectral band-slice level display (see octave_analyzer6e):
-// [z,p,g] = ncauer(Rp,Rs,6); % analog zeros, poles, and gain, where
-// Rp = 80 % dB ripple in stopband
-// Rs = 0.2 % dB ripple in passband
-//
-lowpass6e(fc) =
- tf2s(b21,b11,b01,a11,a01,w1) :
- tf2s(b22,b12,b02,a12,a02,w1) :
- tf2s(b23,b13,b03,a13,a03,w1)
-with {
- b21 = 0.000099999997055;
- a21 = 1;
- b11 = 0;
- a11 = 0.782413046821645;
- b01 = 0.000433227200555;
- a01 = 0.245291508706160;
- b22 = 1;
- a22 = 1;
- b12 = 0;
- a12 = 0.512478641889141;
- b02 = 7.621731298870603;
- a02 = 0.689621364484675;
- b23 = 1;
- a23 = 1;
- b13 = 0;
- a13 = 0.168404871113589;
- b03 = 53.536152954556727;
- a03 = 1.069358407707312;
- w1 = 2*PI*fc;
-};
-
-//===================== Elliptic Highpass Filters =====================
-// USAGE:
-// _ : highpass3e(fc) : _
-// _ : highpass6e(fc) : _
-// where fc = -3dB frequency in Hz
-
-//----------------------------- highpass3e -----------------------------
-// Third-order Elliptic (Cauer) highpass filter
-// DESIGN: Inversion of lowpass3e wrt unit circle in s plane (s <- 1/s)
-//
-highpass3e(fc) = tf2s(b01/a01,b11/a01,b21/a01,a11/a01,1/a01,w1) :
- tf1s(1/a02,0,1/a02,w1)
-with {
- a11 = 0.802636764161030;
- a01 = 1.412270893774204;
- a02 = 0.822445908998816;
- b21 = 0.019809144837789;
- b11 = 0;
- b01 = 1.161516418982696;
- w1 = 2*PI*fc;
-};
-
-//----------------------------- highpass6e -----------------------------
-// Sixth-order Elliptic/Cauer highpass filter
-// METHOD: Inversion of lowpass3e wrt unit circle in s plane (s <- 1/s)
-//
-highpass6e(fc) =
- tf2s(b01/a01,b11/a01,b21/a01,a11/a01,1/a01,w1) :
- tf2s(b02/a02,b12/a02,b22/a02,a12/a02,1/a02,w1) :
- tf2s(b03/a03,b13/a03,b23/a03,a13/a03,1/a03,w1)
-with {
- b21 = 0.000099999997055;
- a21 = 1;
- b11 = 0;
- a11 = 0.782413046821645;
- b01 = 0.000433227200555;
- a01 = 0.245291508706160;
- b22 = 1;
- a22 = 1;
- b12 = 0;
- a12 = 0.512478641889141;
- b02 = 7.621731298870603;
- a02 = 0.689621364484675;
- b23 = 1;
- a23 = 1;
- b13 = 0;
- a13 = 0.168404871113589;
- b03 = 53.536152954556727;
- a03 = 1.069358407707312;
- w1 = 2*PI*fc;
-};
-
-//================== Butterworth Bandpass/Bandstop Filters =====================
-// Order 2*Nh Butterworth bandpass filter made using the transformation
-// s <- s + wc^2/s on lowpass(Nh), where wc is the desired bandpass center
-// frequency. The lowpass(Nh) cutoff w1 is half the desired bandpass width.
-// A notch-like "bandstop" filter is similarly made from highpass(Nh).
-//
-// USAGE:
-// _ : bandpass(Nh,fl,fu) : _
-// _ : bandstop(Nh,fl,fu) : _
-// where
-// Nh = HALF the desired bandpass/bandstop order (which is therefore even)
-// fl = lower -3dB frequency in Hz
-// fu = upper -3dB frequency in Hz
-// Thus, the passband (stopband) width is fu-fl,
-// and its center frequency is (fl+fu)/2.
-//
-// REFERENCE:
-// http://cnx.org/content/m16913/latest/
-//
-bandpass(Nh,fl,fu) = bandpass0_bandstop1(0,Nh,fl,fu);
-bandstop(Nh,fl,fu) = bandpass0_bandstop1(1,Nh,fl,fu);
-bandpass0_bandstop1(s,Nh,fl,fu) = bpbsr(s,Nh,Nh,fl,fu)
-with {
- wl = 2*PI*fl; // digital (z-plane) lower passband edge
- wu = 2*PI*fu; // digital (z-plane) upper passband edge
-
- c = 2.0*SR; // bilinear transform scaling used in tf2sb, tf1sb
- wla = c*tan(wl/c); // analog (s-splane) lower cutoff
- wua = c*tan(wu/c); // analog (s-splane) upper cutoff
-
- wc = sqrt(wla*wua); // s-plane center frequency
- w1 = wua - wc^2/wua; // s-plane lowpass prototype cutoff
-
- bpbsr(s,0,Nh,fl,fu) = _;
- bpbsr(s,1,Nh,fl,fu) = tf1sb(s,1-s,1,w1,wc);
- bpbsr(s,O,Nh,fl,fu) = bpbsr(s,O-2,Nh,fl,fu) : tf2sb(s,0,(1-s),a1s,1,w1,wc)
- with {
- parity = Nh % 2;
- S = (O-parity)/2; // current section number
- a1s = -2*cos(-PI + (1-parity)*PI/(2*Nh) + (S-1+parity)*PI/Nh);
- };
-};
-
-//======================= Elliptic Bandpass Filters ============================
-
-//----------------------------- bandpass6e -----------------------------
-// Order 12 elliptic bandpass filter analogous to bandpass(6) above.
-//
-bandpass6e(fl,fu) = tf2sb(b21,b11,b01,a11,a01,w1,wc) : tf1sb(0,1,a02,w1,wc)
-with {
- a11 = 0.802636764161030; // In octave: format long; poly(p(1:2))
- a01 = 1.412270893774204;
- a02 = 0.822445908998816; // poly(p(3))
- b21 = 0.019809144837789; // poly(z)
- b11 = 0;
- b01 = 1.161516418982696;
-
- wl = 2*PI*fl; // digital (z-plane) lower passband edge
- wu = 2*PI*fu; // digital (z-plane) upper passband edge
-
- c = 2.0*SR; // bilinear transform scaling used in tf2sb, tf1sb
- wla = c*tan(wl/c); // analog (s-splane) lower cutoff
- wua = c*tan(wu/c); // analog (s-splane) upper cutoff
-
- wc = sqrt(wla*wua); // s-plane center frequency
- w1 = wua - wc^2/wua; // s-plane lowpass cutoff
-};
-
-//----------------------------- bandpass12e -----------------------------
-
-bandpass12e(fl,fu) =
- tf2sb(b21,b11,b01,a11,a01,w1,wc) :
- tf2sb(b22,b12,b02,a12,a02,w1,wc) :
- tf2sb(b23,b13,b03,a13,a03,w1,wc)
-with { // octave script output:
- b21 = 0.000099999997055;
- a21 = 1;
- b11 = 0;
- a11 = 0.782413046821645;
- b01 = 0.000433227200555;
- a01 = 0.245291508706160;
- b22 = 1;
- a22 = 1;
- b12 = 0;
- a12 = 0.512478641889141;
- b02 = 7.621731298870603;
- a02 = 0.689621364484675;
- b23 = 1;
- a23 = 1;
- b13 = 0;
- a13 = 0.168404871113589;
- b03 = 53.536152954556727;
- a03 = 1.069358407707312;
-
- wl = 2*PI*fl; // digital (z-plane) lower passband edge
- wu = 2*PI*fu; // digital (z-plane) upper passband edge
-
- c = 2.0*SR; // bilinear transform scaling used in tf2sb, tf1sb
- wla = c*tan(wl/c); // analog (s-splane) lower cutoff
- wua = c*tan(wu/c); // analog (s-splane) upper cutoff
-
- wc = sqrt(wla*wua); // s-plane center frequency
- w1 = wua - wc^2/wua; // s-plane lowpass cutoff
-};
-
-//================= Parametric Equalizers (Shelf, Peaking) ==================
-// REFERENCES
-// - http://en.wikipedia.org/wiki/Equalization
-// - Digital Audio Signal Processing, Udo Zolzer, Wiley, 1999, p. 124
-// - http://www.harmony-central.com/Computer/Programming/Audio-EQ-Cookbook.txt
-// http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
-// - https://ccrma.stanford.edu/~jos/filters/Low_High_Shelving_Filters.html
-// - https://ccrma.stanford.edu/~jos/filters/Peaking_Equalizers.html
-// - maxmsp.lib in the Faust distribution
-// - bandfilter.dsp in the faust2pd distribution
-
-//----------------------------- low_shelf ------------------------------------
-// First-order "low shelf" filter (gain boost|cut between dc and some frequency)
-// USAGE: lowshelf(L0,fx), where
-// L0 = desired boost (dB) between dc and fx
-// fx = desired transition frequency (Hz) from boost to unity gain
-// The gain at SR/2 is constrained to be 1.
-//
-low_shelf = low_shelf3; // default
-low_shelf1(L0,fx,x) = x + (db2linear(L0)-1)*lowpass(1,fx,x);
-low_shelf1_l(G0,fx,x) = x + (G0-1)*lowpass(1,fx,x);
-low_shelf3(L0,fx,x) = x + (db2linear(L0)-1)*lowpass(3,fx,x);
-low_shelf5(L0,fx,x) = x + (db2linear(L0)-1)*lowpass(5,fx,x);
-
-//----------------------------- high_shelf -----------------------------------
-// First-order "high shelf" filter (gain boost|cut above some frequency)
-//
-// USAGE: high_shelf(Lpi,fx), where
-// Lpi = desired boost or cut (dB) between fx and SR/2
-// fx = desired transition frequency in Hz
-// The gain at dc is constrained to be 1
-//
-high_shelf=high_shelf3; //default
-high_shelf1(Lpi,fx,x) = x + (db2linear(Lpi)-1)*highpass(1,fx,x);
-high_shelf1_l(Gpi,fx,x) = x + (Gpi-1)*highpass(1,fx,x);
-high_shelf3(Lpi,fx,x) = x + (db2linear(Lpi)-1)*highpass(3,fx,x);
-high_shelf5(Lpi,fx,x) = x + (db2linear(Lpi)-1)*highpass(5,fx,x);
-
-//-------------------------------- peak_eq -----------------------------------
-// Second order "peaking equalizer" section
-// (gain boost or cut near some frequency)
-// Also called a "parametric equalizer" section
-// USAGE: _ : peak_eq(Lfx,fx,B) : _;
-// where
-// Lfx = level (dB) at fx
-// fx = peak frequency (Hz)
-// B = bandwidth (B) of peak in Hz
-//
-peak_eq(Lfx,fx,B) = tf2s(1,b1s,1,a1s,1,wx) with {
- T = float(1.0/SR);
- Bw = B*T/sin(wx*T); // prewarp s-bandwidth for more accuracy in z-plane
- a1 = PI*Bw;
- b1 = g*a1;
- g = db2linear(abs(Lfx));
- b1s = select2(Lfx>0,a1,b1); // When Lfx>0, pole dominates bandwidth
- a1s = select2(Lfx>0,b1,a1); // When Lfx<0, zero dominates
- wx = 2*PI*fx;
-};
-
-//------------------------------- peak_eq_cq ---------------------------------
-// Constant-Q second order peaking equalizer section
-// USAGE: _ : peak_eq_cq(Lfx,fx,Q) : _;
-// where
-// Lfx = level (dB) at fx
-// fx = boost or cut frequency (Hz)
-// Q = "Quality factor" = fx/B where B = bandwidth of peak in Hz
-//
-peak_eq_cq(Lfx,fx,Q) = peak_eq(Lfx,fx,fx/Q);
-
-//------------------------------- peak_eq_rm ---------------------------------
-// Regalia-Mitra second order peaking equalizer section
-// USAGE: _ : peak_eq_rm(Lfx,fx,tanPiBT) : _;
-// where
-// Lfx = level (dB) at fx
-// fx = boost or cut frequency (Hz)
-// tanPiBT = tan(PI*B/SR), where B = -3dB bandwidth (Hz) when 10^(Lfx/20) = 0
-// ~ PI*B/SR for narrow bandwidths B
-//
-// REFERENCE:
-// P.A. Regalia, S.K. Mitra, and P.P. Vaidyanathan,
-// "The Digital All-Pass Filter: A Versatile Signal Processing Building Block"
-// Proceedings of the IEEE, 76(1):19-37, Jan. 1988. (See pp. 29-30.)
-//
-peak_eq_rm(Lfx,fx,tanPiBT) = _ <: _,A,_ : +,- : *(0.5),*(K/2.0) : + with {
- A = tf2(k2, k1*(1+k2), 1, k1*(1+k2), k2) <: _,_; // allpass
- k1 = 0.0 - cos(2.0*PI*fx/SR);
- k2 = (1.0 - tanPiBT)/(1.0 + tanPiBT);
- K = db2linear(Lfx);
-};
-
-//-------------------------- parametric_eq_demo ------------------------------
-// USAGE: _ : parametric_eq_demo: _ ;
-parametric_eq_demo = // input_signal :
- low_shelf(LL,FL) :
- peak_eq(LP,FP,BP) :
- high_shelf(LH,FH)
-// Recommended:
-// : mth_octave_spectral_level_demo(2) // half-octave spectrum analyzer
-with {
- eq_group(x) = hgroup("[0] PARAMETRIC EQ SECTIONS
- [tooltip: See Faust's filter.lib for info and pointers]",x);
-
- ls_group(x) = eq_group(vgroup("[1] Low Shelf",x));
- LL = ls_group(hslider("[0] Low Boost|Cut [unit:dB] [style:knob]
- [tooltip: Amount of low-frequency boost or cut in decibels]",
- 0,-40,40,0.1));
- FL = ls_group(hslider("[1] Transition Frequency [unit:Hz] [style:knob]
- [tooltip: Transition-frequency from boost (cut) to unity gain]",
- 200,1,5000,1));
-
- pq_group(x) = eq_group(vgroup("[2] Peaking Equalizer
- [tooltip: Parametric Equalizer sections from filter.lib]",x));
- LP = pq_group(hslider("[0] Peak Boost|Cut [unit:dB] [style:knob]
- [tooltip: Amount of local boost or cut in decibels]",
- 0,-40,40,0.1));
- FP = pq_group(hslider("[1] Peak Frequency [unit:PK] [style:knob]
- [tooltip: Peak Frequency in Piano Key (PK) units (A-440= 49 PK)]",
- 49,1,100,1)) : smooth(0.999) : pianokey2hz
- with { pianokey2hz(x) = 440.0*pow(2.0, (x-49.0)/12); };
-
- Q = pq_group(hslider("[2] Peak Q [style:knob]
- [tooltip: Quality factor (Q) of the peak = center-frequency/bandwidth]",
- 40,1,50,0.1));
-
- BP = FP/Q;
-
- hs_group(x) = eq_group(vgroup("[3] High Shelf
- [tooltip: A high shelf provides a boost or cut
- above some frequency]",x));
- LH = hs_group(hslider("[0] High Boost|Cut [unit:dB] [style:knob]
- [tooltip: Amount of high-frequency boost or cut in decibels]",
- 0,-40,40,.1));
- FH = hs_group(hslider("[1] Transition Frequency [unit:Hz] [style:knob]
- [tooltip: Transition-frequency from boost (cut) to unity gain]",
- 8000,20,10000,1));
-};
-
-//========================= Lagrange Interpolation ========================
-// Reference:
-// https://ccrma.stanford.edu/~jos/pasp/Lagrange_Interpolation.html
-//
-//------------------ fdelay1, fdelay2, fdelay3, fdelay4 ---------------
-// Delay lines interpolated using Lagrange interpolation
-// USAGE: _ : fdelayN(maxdelay, delay, inputsignal) : _
-// (exactly like fdelay in music.lib)
-// where N=1,2,3, or 4 is the order of the Lagrange interpolation polynomial.
-//
-// NOTE: requested delay should not be less than (N-1)/2.
-//
-// NOTE: While the implementations below appear to use multiple delay lines,
-// they in fact use only one thanks to optimization by the Faust compiler.
-
-// first-order case (linear interpolation) - equivalent to fdelay in music.lib
-// delay d in [0,1]
-fdelay1(n,d,x) = delay(n,id,x)*(1 - fd) + delay(n,id+1,x)*fd
-with {
- id = int(d);
- fd = frac(d);
-};
-
-// second-order (quadratic) case, delay in [0.5,1.5]
-// delay d should be at least 0.5
-fdelay2(n,d,x) = delay(n,id,x)*(1-fd)*(2-fd)/2
- + delay(n,id+1,x)*(2-fd)*fd
- + delay(n,id+2,x)*(fd-1)*fd/2
-with {
- o = 0.49999; // offset to make life easy for interpolator
- dmo = d - o; // assumed nonnegative
- id = int(dmo);
- fd = o + frac(dmo);
-};
-
-// third-order (cubic) case, delay in [1,2]
-// delay d should be at least 1
-fdelay3(n,d,x) = delay(n,id,x) * (0-fdm1*fdm2*fdm3)/6
- + delay(n,id+1,x) * fd*fdm2*fdm3/2
- + delay(n,id+2,x) * (0-fd*fdm1*fdm3)/2
- + delay(n,id+3,x) * fd*fdm1*fdm2/6
-with {
- id = int(d-1);
- fd = 1+frac(d);
- fdm1 = fd-1;
- fdm2 = fd-2;
- fdm3 = fd-3;
-};
-
-// fourth-order (quartic) case, delay in [1.5,2.5]
-// delay d should be at least 1.5
-fdelay4(n,d,x) = delay(n,id,x) * fdm1*fdm2*fdm3*fdm4/24
- + delay(n,id+1,x) * (0-fd*fdm2*fdm3*fdm4)/6
- + delay(n,id+2,x) * fd*fdm1*fdm3*fdm4/4
- + delay(n,id+3,x) * (0-fd*fdm1*fdm2*fdm4)/6
- + delay(n,id+4,x) * fd*fdm1*fdm2*fdm3/24
-with {
-//v1: o = 1;
- o = 1.49999;
- dmo = d - o; // assumed nonnegative
- id = int(dmo);
- fd = o + frac(dmo);
- fdm1 = fd-1;
- fdm2 = fd-2;
- fdm3 = fd-3;
- fdm4 = fd-4;
-};
-
-// fifth-order case, delay in [2,3]
-// delay d should be at least 2
-fdelay5(n,d,x) =
- delay(n,id,x) * -fdm1*fdm2*fdm3*fdm4*fdm5/120
- + delay(n,id+1,x) * fd* fdm2*fdm3*fdm4*fdm5/24
- + delay(n,id+2,x) * -fd*fdm1* fdm3*fdm4*fdm5/12
- + delay(n,id+3,x) * fd*fdm1*fdm2* fdm4*fdm5/12
- + delay(n,id+4,x) * -fd*fdm1*fdm2*fdm3* fdm5/24
- + delay(n,id+5,x) * fd*fdm1*fdm2*fdm3*fdm4 /120
-with {
-//v1: o = 1;
- o = 1.99999;
- dmo = d - o; // assumed nonnegative
- id = int(dmo);
- fd = o + frac(dmo);
- fdm1 = fd-1;
- fdm2 = fd-2;
- fdm3 = fd-3;
- fdm4 = fd-4;
- fdm5 = fd-5;
-};
-
-//====================== Thiran Allpass Interpolation =====================
-// Reference:
-// https://ccrma.stanford.edu/~jos/pasp/Thiran_Allpass_Interpolators.html
-//
-//---------------- fdelay1a, fdelay2a, fdelay3a, fdelay4a -------------
-// Delay lines interpolated using Thiran allpass interpolation
-// USAGE: fdelayNa(maxdelay, delay, inputsignal)
-// (exactly like fdelay in music.lib)
-// where N=1,2,3, or 4 is the order of the Thiran interpolation filter,
-// and the delay argument is at least N - 1/2.
-//
-// (Move the following and similar notes above to filter-lib-doc.txt?)
-//
-// NOTE: The interpolated delay should not be less than N - 1/2.
-// (The allpass delay ranges from N - 1/2 to N + 1/2.)
-// This constraint can be alleviated by altering the code,
-// but be aware that allpass filters approach zero delay
-// by means of pole-zero cancellations.
-// The delay range [N-1/2,N+1/2] is not optimal. What is?
-//
-// NOTE: Delay arguments too small will produce an UNSTABLE allpass!
-//
-// NOTE: Because allpass interpolation is recursive, it is not as robust
-// as Lagrange interpolation under time-varying conditions.
-// (You may hear clicks when changing the delay rapidly.)
-//
-// first-order allpass interpolation, delay d in [0.5,1.5]
-fdelay1a(n,d,x) = delay(n,id,x) : tf1(eta,1,eta)
-with {
- o = 0.49999; // offset to make life easy for allpass
- dmo = d - o; // assumed nonnegative
- id = int(dmo);
- fd = o + frac(dmo);
- eta = (1-fd)/(1+fd); // allpass coefficient
-};
-
-// second-order allpass delay in [1.5,2.5]
-fdelay2a(n,d,x) = delay(n,id,x) : tf2(a2,a1,1,a1,a2)
-with {
- o = 1.49999;
- dmo = d - o; // delay range is [order-1/2, order+1/2]
- id = int(dmo);
- fd = o + frac(dmo);
- a1o2 = (2-fd)/(1+fd); // share some terms (the compiler does this anyway)
- a1 = 2*a1o2;
- a2 = a1o2*(1-fd)/(2+fd);
-};
-
-// third-order allpass delay in [2.5,3.5]
-// delay d should be at least 2.5
-fdelay3a(n,d,x) = delay(n,id,x) : iir((a3,a2,a1,1),(a1,a2,a3))
-with {
- o = 2.49999;
- dmo = d - o;
- id = int(dmo);
- fd = o + frac(dmo);
- a1o3 = (3-fd)/(1+fd);
- a2o3 = a1o3*(2-fd)/(2+fd);
- a1 = 3*a1o3;
- a2 = 3*a2o3;
- a3 = a2o3*(1-fd)/(3+fd);
-};
-
-// fourth-order allpass delay in [3.5,4.5]
-// delay d should be at least 3.5
-fdelay4a(n,d,x) = delay(n,id,x) : tf4(a4,a3,a2,a1,1,a1,a2,a3,a4)
-with {
- o = 3.49999;
- dmo = d - o;
- id = int(dmo);
- fd = o + frac(dmo);
- a1o4 = (4-fd)/(1+fd);
- a2o6 = a1o4*(3-fd)/(2+fd);
- a3o4 = a2o6*(2-fd)/(3+fd);
- a1 = 4*a1o4;
- a2 = 6*a2o6;
- a3 = 4*a3o4;
- a4 = a3o4*(1-fd)/(4+fd);
-};
-
-//================ Mth-Octave Filter-Banks and Spectrum-Analyzers ============
-// Mth-octave filter-banks and spectrum-analyzers split the input signal into a
-// bank of parallel signals, one for each spectral band. The parameters are
-//
-// M = number of band-slices per octave (>1)
-// N = total number of bands (>2)
-// ftop = upper bandlimit of the Mth-octave bands (<SR/2)
-//
-// In addition to the Mth-octave output signals, there is a highpass signal
-// containing frequencies from ftop to SR/2, and a "dc band" lowpass signal
-// containing frequencies from 0 (dc) up to the start of the Mth-octave bands.
-// Thus, the N output signals are
-//
-// highpass(ftop), MthOctaveBands(M,N-2,ftop), dcBand(ftop*2^(-M*(N-1)))
-//
-// A FILTER-BANK is defined here as a signal bandsplitter having the
-// property that summing its output signals gives an allpass-filtered
-// version of the filter-bank input signal. A more conventional term for
-// this is an "allpass-complementary filter bank". If the allpass filter
-// is a pure delay (and possible scaling), the filter bank is said to be
-// a "perfect-reconstruction filter bank" (see Vaidyanathan-1993 cited
-// below for details). A "graphic equalizer", in which band signals
-// are scaled by gains and summed, should be based on a filter bank.
-//
-// A SPECTRUM-ANALYZER is defined here as any band-split whose bands span
-// the relevant spectrum, but whose band-signals do not
-// necessarily sum to the original signal, either exactly or to within an
-// allpass filtering. Spectrum analyzer outputs are normally at least nearly
-// "power complementary", i.e., the power spectra of the individual bands
-// sum to the original power spectrum (to within some negligible tolerance).
-//
-// The filter-banks below are implemented as Butterworth or Elliptic
-// spectrum-analyzers followed by delay equalizers that make them
-// allpass-complementary.
-//
-// INCREASING CHANNEL ISOLATION
-// Go to higher filter orders - see Regalia et al. or Vaidyanathan (cited
-// below) regarding the construction of more aggressive recursive
-// filter-banks using elliptic or Chebyshev prototype filters.
-//
-// REFERENCES
-// - "Tree-structured complementary filter banks using all-pass sections",
-// Regalia et al., IEEE Trans. Circuits & Systems, CAS-34:1470-1484, Dec. 1987
-// - "Multirate Systems and Filter Banks", P. Vaidyanathan, Prentice-Hall, 1993
-// - Elementary filter theory: https://ccrma.stanford.edu/~jos/filters/
-//
-//------------------------- mth_octave_analyzer ----------------------------
-//
-// USAGE
-// _ : mth_octave_analyzer(O,M,ftop,N) : par(i,N,_); // Oth-order Butterworth
-// _ : mth_octave_analyzer6e(M,ftop,N) : par(i,N,_); // 6th-order elliptic
-//
-// where
-// O = order of filter used to split each frequency band into two
-// M = number of band-slices per octave
-// ftop = highest band-split crossover frequency (e.g., 20 kHz)
-// N = total number of bands (including dc and Nyquist)
-//
-// ACKNOWLEDGMENT
-// Recursive band-splitting formulation improved by Yann Orlarey.
-
-mth_octave_analyzer6e(M,ftop,N) = _ <: bsplit(N-1) with {
- fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
- lp(n) = lowpass6e(fc(n)); // 6th-order elliptic - see other choices above
- hp(n) = highpass6e(fc(n)); // (search for lowpass* and highpass*)
- bsplit(0) = _;
- bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
-};
-
-// Butterworth analyzers may be cascaded with allpass
-// delay-equalizers to make (allpass-complementary) filter banks:
-
-mth_octave_analyzer(O,M,ftop,N) = _ <: bsplit(N-1) with {
- fc(n) = ftop * 2^(float(n-N+1)/float(M));
- lp(n) = lowpass(O,fc(n)); // Order O Butterworth
- hp(n) = highpass(O,fc(n));
- bsplit(0) = _;
- bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
-};
-
-mth_octave_analyzer3(M,ftop,N) = mth_octave_analyzer(3,M,ftop,N);
-mth_octave_analyzer5(M,ftop,N) = mth_octave_analyzer(5,M,ftop,N);
-mth_octave_analyzer_default = mth_octave_analyzer6e; // default analyzer
-
-//------------------------ mth_octave_filterbank -------------------------
-// Allpass-complementary filter banks based on Butterworth band-splitting.
-// For Butterworth band-splits, the needed delay equalizer is easily found.
-
-mth_octave_filterbank(O,M,ftop,N) =
- mth_octave_analyzer(O,M,ftop,N) :
- delayeq(N) with {
- fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
- ap(n) = highpass_plus_lowpass(O,fc(n)); // delay-equalizing allpass
- delayeq(N) = par(i,N-2,apchain(i+1)), _, _;
- apchain(i) = seq(j,N-1-i,ap(j+1));
-};
-
-// dc-inverted version. This reduces the delay-equalizer order for odd O.
-// Negating the input signal makes the dc band noninverting
-// and all higher bands sign-inverted (if preferred).
-mth_octave_filterbank_alt(O,M,ftop,N) =
- mth_octave_analyzer(O,M,ftop,N) : delayeqi(O,N) with {
- fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
- ap(n) = highpass_minus_lowpass(O,fc(n)); // half the order of 'plus' case
- delayeqi(N) = par(i,N-2,apchain(i+1)), _, *(-1.0);
- apchain(i) = seq(j,N-1-i,ap(j+1));
-};
-
-// Note that even-order cases require complex coefficients.
-// See Vaidyanathan 1993 and papers cited there for more info.
-
-mth_octave_filterbank3(M,ftop,N) = mth_octave_filterbank_alt(3,M,ftop,N);
-mth_octave_filterbank5(M,ftop,N) = mth_octave_filterbank(5,M,ftop,N);
-
-mth_octave_filterbank_default = mth_octave_filterbank5;
-
-//======================= Mth-Octave Spectral Level =========================
-// Spectral Level: Display (in bar graphs) the average signal level in each
-// spectral band.
-//
-//------------------------ mth_octave_spectral_level -------------------------
-// USAGE: _ : mth_octave_spectral_level(M,ftop,NBands,tau,dB_offset);
-// where
-// M = bands per octave
-// ftop = lower edge frequency of top band
-// NBands = number of passbands (including highpass and dc bands),
-// tau = spectral display averaging-time (time constant) in seconds,
-// dB_offset = constant dB offset in all band level meters.
-//
-mth_octave_spectral_level6e(M,ftop,N,tau,dB_offset) = _<:
- _,mth_octave_analyzer6e(M,ftop,N) :
- _,(display:>_):attach with {
- display = par(i,N,dbmeter(i));
- dbmeter(i) = abs : smooth(tau2pole(tau)) : linear2db : +(dB_offset) :
- meter(N-i-1);
- meter(i) = speclevel_group(vbargraph("[%2i] [unit:dB]
- [tooltip: Spectral Band Level in dB]", -50, 10));
- // Can M be included in the label string somehow?
- speclevel_group(x) = hgroup("[0] CONSTANT-Q SPECTRUM ANALYZER (6E)
- [tooltip: See Faust's filter.lib for documentation and references]", x);
-};
-
-mth_octave_spectral_level_default = mth_octave_spectral_level6e;
-spectral_level = mth_octave_spectral_level(2,10000,20); // simplest case
-
-//---------------------- mth_octave_spectral_level_demo ----------------------
-// Demonstrate mth_octave_spectral_level in a standalone GUI.
-//
-// USAGE: _ : mth_octave_spectral_level_demo(BandsPerOctave);
-
-mth_octave_spectral_level_demo(M) =
- mth_octave_spectral_level_default(M,ftop,N,tau,dB_offset)
-with {
- // Span nearly 10 octaves so that lowest band-edge is at
- // ftop*2^(-Noct+2) = 40 Hz when ftop=10 kHz:
- N = int(10*M); // without 'int()', segmentation fault observed for M=1.67
- ftop = 10000;
- ctl_group(x) = hgroup("[1] SPECTRUM ANALYZER CONTROLS", x);
- tau = ctl_group(hslider("[0] Level Averaging Time [unit:sec]
- [tooltip: band-level averaging time in seconds]",
- 0.1,0,1,0.01));
- dB_offset = ctl_group(hslider("[1] Level dB Offset [unit:dB]
- [tooltip: Level offset in decibels]",
- 50,0,100,1));
-};
-
-spectral_level_demo = mth_octave_spectral_level_demo(1.5); // 2/3 octave
-
-//---------------- (third|half)_octave_(analyzer|filterbank) -----------------
-
-// Named special cases of mth_octave_* with defaults filled in:
-
-third_octave_analyzer(N) = mth_octave_analyzer_default(3,10000,N);
-third_octave_filterbank(N) = mth_octave_filterbank_default(3,10000,N);
-// Third-Octave Filter-Banks have been used in audio for over a century.
-// See, e.g.,
-// Acoustics [the book], by L. L. Beranek
-// Amer. Inst. Physics for the Acoustical Soc. America,
-// http://asa.aip.org/publications.html, 1986 (1st ed.1954)
-
-// Third-octave bands across the audio spectrum are too wide for current
-// typical computer screens, so half-octave bands are the default:
-half_octave_analyzer(N) = mth_octave_analyzer_default(2,10000,N);
-half_octave_filterbank(N) = mth_octave_filterbank_default(2,10000,N);
-
-octave_filterbank(N) = mth_octave_filterbank_default(1,10000,N);
-octave_analyzer(N) = mth_octave_analyzer_default(1,10000,N);
-
-//=========================== Filter-Bank Demos ==============================
-// Graphic Equalizer: Each filter-bank output signal routes through a fader.
-//
-// USAGE: _ : mth_octave_filterbank_demo(M) : _
-// where
-// M = number of bands per octave
-
-mth_octave_filterbank_demo(M) = bp1(bp,mthoctavefilterbankdemo) with {
- bp1 = component("effect.lib").bypass1;
- mofb_group(x) = vgroup("CONSTANT-Q FILTER BANK (Butterworth dyadic tree)
- [tooltip: See Faust's filter.lib for documentation and references]", x);
- bypass_group(x) = mofb_group(hgroup("[0]", x));
- slider_group(x) = mofb_group(hgroup("[1]", x));
- N = 10*M; // total number of bands (highpass band, octave-bands, dc band)
- ftop = 10000;
- mthoctavefilterbankdemo = chan;
- chan = mth_octave_filterbank_default(M,ftop,N) :
- sum(i,N,(*(db2linear(fader(N-i)))));
- fader(i) = slider_group(vslider("[%2i] [unit:dB]
- [tooltip: Bandpass filter gain in dB]", -10, -70, 10, 0.1)) :
- smooth(0.999);
- bp = bypass_group(checkbox("[0] Bypass
- [tooltip: When this is checked, the filter-bank has no effect]"));
-};
-
-filterbank_demo = mth_octave_filterbank_demo(1); // octave-bands = default
-
-//=========== Arbritary-Crossover Filter-Banks and Spectrum Analyzers ========
-// These are similar to the Mth-octave filter-banks above, except that the
-// band-split frequencies are passed explicitly as arguments.
-//
-// USAGE:
-// _ : filterbank (O,freqs) : par(i,N,_); // Butterworth band-splits
-// _ : filterbanki(O,freqs) : par(i,N,_); // Inverted-dc version
-// _ : analyzer (O,freqs) : par(i,N,_); // No delay equalizer
-//
-// where
-// O = band-split filter order (ODD integer required for filterbank[i])
-// freqs = (fc1,fc2,...,fcNs) [in numerically ascending order], where
-// Ns=N-1 is the number of octave band-splits
-// (total number of bands N=Ns+1).
-//
-// If frequencies are listed explicitly as arguments, enclose them in parens:
-//
-// _ : filterbank(3,(fc1,fc2)) : _,_,_
-//
-// ACKNOWLEDGMENT
-// Technique for processing a variable number of signal arguments due
-// to Yann Orlarey (as is the entire Faust framework!)
-//
-//------------------------------ analyzer --------------------------------------
-analyzer(O,lfreqs) = _ <: bsplit(nb) with
-{
- nb = count(lfreqs);
- fc(n) = take(n, lfreqs);
- lp(n) = lowpass(O,fc(n));
- hp(n) = highpass(O,fc(n));
- bsplit(0) = _;
- bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
-};
-
-//----------------------------- filterbank -------------------------------------
-filterbank(O,lfreqs) = analyzer(O,lfreqs) : delayeq with
-{
- nb = count(lfreqs);
- fc(n) = take(n, lfreqs);
- ap(n) = highpass_plus_lowpass(O,fc(n));
- delayeq = par(i,nb-1,apchain(nb-1-i)),_,_;
- apchain(0) = _;
- apchain(i) = ap(i) : apchain(i-1);
-};
-
-//----------------------------- filterbanki ------------------------------------
-filterbanki(O,lfreqs) = _ <: bsplit(nb) with
-{
- fc(n) = take(n, lfreqs);
- lp(n) = lowpass(O,fc(n));
- hp(n) = highpass(O,fc(n));
- ap(n) = highpass_minus_lowpass(O,fc(n));
- bsplit(0) = *(-1.0);
- bsplit(i) = (hp(i) : delayeq(i-1)), (lp(i) <: bsplit(i-1));
- delayeq(0) = _; // moving the *(-1) here inverts all outputs BUT dc
- delayeq(i) = ap(i) : delayeq(i-1);
-};
projects / Faustine.git / blobdiff