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