X-Git-Url: https://scm.cri.ensmp.fr/git/Faustine.git/blobdiff_plain/c7f552fd8888da2f0d8cfb228fe0f28d3df3a12c..b4b6f2ea75b9f0f3ca918f5b84016610bf7a4d4f:/interpretor/preprocessor/faust-0.9.47mr3/architecture/filter.lib diff --git a/interpretor/preprocessor/faust-0.9.47mr3/architecture/filter.lib b/interpretor/preprocessor/faust-0.9.47mr3/architecture/filter.lib new file mode 100644 index 0000000..db89294 --- /dev/null +++ b/interpretor/preprocessor/faust-0.9.47mr3/architecture/filter.lib @@ -0,0 +1,1581 @@ +// 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 (_):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); +};