1 // filter.lib - digital filters of various types useful in audio and beyond
3 declare name "Faust Filter Library";
4 declare author "Julius O. Smith (jos at ccrma.stanford.edu)";
5 declare copyright "Julius O. Smith III";
6 declare version "1.29";
7 declare license "STK-4.3"; // Synthesis Tool Kit 4.3 (MIT style license)
8 declare reference "https://ccrma.stanford.edu/~jos/filters/";
10 import("music.lib"); // delay, frac and, from math.lib, SR and PI
12 //---------------------- zero(z) --------------------------
13 // z = location of zero along real axis in z-plane
14 // Difference equation: y(n) = x(n) - z * x(n-1)
15 // Reference: https://ccrma.stanford.edu/~jos/filters/One_Zero.html
17 zero(z) = _ <: _,mem : _,*(z) : -;
19 //------------------------ pole(p) ---------------------------
20 // p = pole location = feedback coefficient
21 // Could also be called a "leaky integrator".
22 // Difference equation: y(n) = x(n) + p * y(n-1)
23 // Reference: https://ccrma.stanford.edu/~jos/filters/One_Pole.html
25 pole(p) = + ~ *(p);
27 //---------------------- integrator --------------------------
28 // pole(1) [implemented separately for block-diagram clarity]
30 integrator = + ~ _ ;
32 //----------------------- tau2pole ---------------------------
33 // tau2pole(tau) returns a real pole giving exponential decay with
34 // tau = time-constant in seconds
35 //
36 tau2pole(tau) = exp(-1.0/(tau*SR));
38 //---------------------- smooth(s) --------------------------
39 // Exponential smoothing by a unity-dc-gain one-pole lowpass
40 //
41 // USAGE: smooth(tau2pole(tau)), where
42 // tau = desired smoothing time constant in seconds,
43 // or
44 // smooth(s), where s = smoothness between 0 and 1.
45 // s=0 for no smoothing
46 // s=0.999 is "very smooth"
47 // s>1 is unstable, and s=1 yields the zero signal for all inputs.
48 // The exponential time-constant is approximately
49 // 1/(1-s) samples, when s is close to (but less than) 1.
50 // Reference:
53 smooth(s) = *(1.0 - s) : + ~ *(s);
55 //------------------- dcblockerat(fb) -----------------------
56 // fb = "break frequency" in Hz, i.e., -3 dB gain frequency.
57 // The amplitude response is substantially flat above fb,
58 // and sloped at about +6 dB/octave below fb.
59 // Derived from the analog transfer function
60 // H(s) = s / (s + 2*PI*fb)
61 // by the low-frequency-matching bilinear transform method
62 // (i.e., the standard frequency-scaling constant 2*SR).
63 // Reference:
64 // https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html
66 dcblockerat(fb) = *(b0) : zero(1) : pole(p)
67 with {
68 wn = PI*fb/SR;
69 b0 = 1.0 / (1 + wn);
70 p = (1 - wn) * b0;
71 };
73 //---------------------- dcblocker --------------------------
74 // Default dc blocker has -3dB point near 35 Hz (at 44.1 kHz)
75 // and high-frequency gain near 1.0025 (due to no scaling)
76 //
77 dcblocker = zero(1) : pole(0.995);
79 //------------ notchw(width,freq), notch(freq) --------------
80 // width = "notch width" in Hz (approximate)
81 // freq = "notch frequency" in Hz
82 // Reference:
83 // https://ccrma.stanford.edu/~jos/pasp/Phasing_2nd_Order_Allpass_Filters.html
85 notchw(width,freq) = tf2(b0,b1,b2,a1,a2)
86 with {
87 fb = 0.5*width; // First design a dcblockerat(width/2)
88 wn = PI*fb/SR;
89 b0db = 1.0 / (1 + wn);
90 p = (1 - wn) * b0db; // This is our pole radius.
91 // Now place unit-circle zeros at desired angles:
92 tn = 2*PI*freq/SR;
93 a2 = p * p;
94 a2p1 = 1+a2;
95 a1 = -a2p1*cos(tn);
96 b1 = a1;
97 b0 = 0.5*a2p1;
98 b2 = b0;
99 };
101 //========================= Comb Filters ===============================
103 //----------------------- ff_comb, ff_fcomb ----------------------------
104 // Feed-Forward Comb Filter
105 //
106 // USAGE:
107 // _ : ff_comb(maxdel,intdel,b0,bM) : _
108 // _ : ff_fcomb(maxdel,del,b0,bM) : _
109 // where
110 // maxdel = maximum delay (a power of 2)
111 // intdel = current (integer) comb-filter delay between 0 and maxdel
112 // del = current (float) comb-filter delay between 0 and maxdel
113 // b0 = gain applied to delay-line input
114 // bM = gain applied to delay-line output and then summed with input
115 //
116 // Note that ff_comb requires integer delays (uses delay() internally)
117 // while ff_fcomb takes floating-point delays (uses fdelay() internally).
118 //
119 // REFERENCE:
120 // https://ccrma.stanford.edu/~jos/pasp/Feedforward_Comb_Filters.html
122 ff_comb (maxdel,M,b0,bM) = _ <: *(b0), bM * delay(maxdel,M) : + ;
123 ff_fcomb(maxdel,M,b0,bM) = _ <: *(b0), bM * fdelay(maxdel,M) : + ;
125 // Typical special case:
126 ffcombfilter(maxdel,del,g) = ff_comb(maxdel,del,1,g);
128 //----------------------- fb_comb, fb_fcomb, rev1 -----------------------
129 // Feed-Back Comb Filter
130 //
131 // USAGE:
132 // _ : fb_comb(maxdel,intdel,b0,aN) : _
133 // _ : fb_fcomb(maxdel,del,b0,aN) : _
134 // _ : rev1(maxdel,del,-aN) : _
135 // where
136 // maxdel = maximum delay (a power of 2)
137 // intdel = current (integer) comb-filter delay between 0 and maxdel
138 // del = current (float) comb-filter delay between 0 and maxdel
139 // b0 = gain applied to delay-line input and forwarded to output
140 // aN = minus the gain applied to delay-line output before
141 // summing with the input and feeding to the delay line
142 //
143 // Reference:
144 // https://ccrma.stanford.edu/~jos/pasp/Feedback_Comb_Filters.html
146 fb_comb (maxdel,N,b0,aN) = (+ <: delay(maxdel,N),_) ~ *(-aN) : !,*(b0);
147 fb_fcomb(maxdel,N,b0,aN) = (+ <: fdelay(maxdel,N),_) ~ *(-aN) : !,*(b0);
149 // The "rev1 section" dates back to the 1960s in computer-music reverberation.
150 // See the jcrev and brassrev in effect.lib for usage examples.
151 rev1(maxdel,N,g) = fb_comb (maxdel,N,1,-g);
153 // Typical special case:
154 fbcombfilter(maxdel,intdel,g) = (+ : delay(maxdel,intdel)) ~ *(g);
155 ffbcombfilter(maxdel,del,g) = (+ : fdelay(maxdel,del)) ~ *(g);
157 //------------------- allpass_comb, allpass_fcomb, rev2 -----------------
158 // Schroeder Allpass Comb Filter
159 //
160 // USAGE:
161 // _ : allpass_comb (maxdel,intdel,aN) : _
162 // _ : allpass_fcomb(maxdel,del,aN) : _
163 // _ : rev2(maxdel,del,-aN) : _
164 // where
165 // maxdel = maximum delay (a power of 2)
166 // intdel = current (integer) comb-filter delay between 0 and maxdel
167 // del = current (float) comb-filter delay between 0 and maxdel
168 // aN = minus the feedback gain
169 //
170 // Note that allpass_comb(maxlen,len,aN) =
171 // ff_comb(maxlen,len,aN,1) :
172 // fb_comb(maxlen,len-1,1,aN);
173 // which is a direct-form-1 implementation, requiring two delay lines.
174 // The implementation here is direct-form-2 requiring only one delay line.
175 //
176 // REFERENCES:
177 // https://ccrma.stanford.edu/~jos/pasp/Allpass_Two_Combs.html
178 // https://ccrma.stanford.edu/~jos/pasp/Schroeder_Allpass_Sections.html
179 // https://ccrma.stanford.edu/~jos/filters/Four_Direct_Forms.html
181 allpass_comb(maxdel,N,aN) = (+ <:
182 delay(maxdel,N-1),*(aN)) ~ *(-aN)
183 : mem,_ : + ;
185 // The "rev2 section" dates back to the 1960s in computer-music reverberation:
186 rev2(maxlen,len,g) = allpass_comb(maxlen,len,-g);
188 //================ Direct-Form Digital Filter Sections ================
190 // Specified by transfer-function polynomials B(z)/A(z) as in matlab
192 //---------------------------- iir (tfN) -------------------------------
193 // Nth-order Infinite-Impulse-Response (IIR) digital filter,
194 // implemented in terms of the Transfer-Function (TF) coefficients.
195 // Such filter structures are termed "direct form".
196 //
197 // USAGE:
198 // _ : iir(bcoeffs,acoeffs) : _
199 // where
200 // order = filter order (int) = max(#poles,#zeros)
201 // bcoeffs = (b0,b1,...,b_order) = TF numerator coefficients
202 // acoeffs = (a1,...,a_order) = TF denominator coeffs (a0=1)
203 //
204 // REFERENCE:
205 // https://ccrma.stanford.edu/~jos/filters/Four_Direct_Forms.html
207 iir(bv,av) = sub ~ fir(av) : fir(bv);
209 //----------------------------- sub ---------------------------------
210 sub(x,y) = y-x; // move to math.lib?
212 //----------------------------- fir ---------------------------------
213 // FIR filter (convolution of FIR filter coefficients with a signal)
214 //
215 // USAGE:
216 // _ : fir(bv) : _
217 // where bv = b0,b1,...,bn is a parallel bank of coefficient signals.
218 // NOTE: bv is processed using pattern-matching at compile time,
219 // so it must have this normal form (parallel signals).
220 // EXAMPLE: Smoothing white noise with a five-point moving average:
221 // bv = .2,.2,.2,.2,.2;
222 // process = noise : fir(bv);
223 // EQUIVALENT (note double parens):
224 // process = noise : fir((.2,.2,.2,.2,.2));
226 fir(bv) = conv(bv);
228 //--------------------------- conv, convN -------------------------------
229 // Convolution of input signal with given coefficients
230 //
231 // USAGE:
232 // _ : conv((k1,k2,k3,...,kN)) : _; // Argument = one signal bank
233 // _ : convN(N,(k1,k2,k3,...)) : _; // Useful when N < count((k1,...))
235 convN(N,kv,x) = sum(i,N,take(i+1,kv) * x@i); // take() defined in math.lib
237 conv(kv,x) = sum(i,count(kv),take(i+1,kv) * x@i); // count() from math.lib
239 // Named special cases:
240 //----------------------------- tf1, tf2 ---------------------------------
241 // tfN = N'th-order direct-form digital filter
242 tf1(b0,b1,a1) = _ <: *(b0), (mem : *(b1)) :> + ~ *(0-a1);
243 tf2(b0,b1,b2,a1,a2) = iir((b0,b1,b2),(a1,a2)); // cf. TF2 in music.lib)
245 //===================== Ladder/Lattice Digital Filters ======================
246 // Ladder and lattice digital filters generally have superior numerical
247 // properties relative to direct-form digital filters. They can be derived
248 // from digital waveguide filters, which gives them a physical interpretation.
250 // REFERENCES:
251 // F. Itakura and S. Saito: "Digital Filtering Techniques for Speech Analysis and Synthesis",
252 // 7th Int. Cong. Acoustics, Budapest, 25 C 1, 1971.
253 // J. D. Markel and A. H. Gray: Linear Prediction of Speech, New York: Springer Verlag, 1976.
256 //------------------------------ block, crossn,crossn1 ----------------------------------
257 // signal block/crossing utilities
258 // (move to math.lib?)
260 // block - terminate n signals (goes with bus(n) in math.lib)
261 block(n) = par(i,n,!);
263 // crossnn - cross two bus(n)s:
264 crossnn(n) = bus(n),bus(n) <: block(n),bus(n),bus(n),block(n);
266 // crossn1 - cross bus(n) and bus(1):
267 crossn1(n) = bus(n),(bus(1)<:bus(n)) <: block(n),bus(n),bus(n),block(n):bus(1),block(n-1),bus(n);
269 //------------------------------- av2sv -----------------------------------
270 // Compute reflection coefficients sv from transfer-function denominator av
271 //
272 // USAGE:
273 // sv = av2sv(av)
274 // where
275 // av = parallel signal bank a1,...,aN
276 // sv = parallel signal bank s1,...,sN
277 // where si = ith reflection coefficient, and
278 // ai = coefficient of z^(-i) in the filter
279 // transfer-function denominator A(z).
280 //
281 // REFERENCE:
282 // https://ccrma.stanford.edu/~jos/filters/Step_Down_Procedure.html
283 // (where reflection coefficients are denoted by k rather than s).
285 av2sv(av) = par(i,M,s(i+1)) with {
286 M = count(av);
287 s(m) = sr(M-m+1); // m=1..M
288 sr(m) = Ari(m,M-m+1); // s_{M-1-m}
289 Ari(m,i) = take(i+1,Ar(m-1));
290 //step-down recursion for lattice/ladder digital filters:
291 Ar(0) = (1,av); // Ar(m) is order M-m (i.e. "reverse-indexed")
292 Ar(m) = 1,par(i,M-m, (Ari(m,i+1) - sr(m)*Ari(m,M-m-i))/(1-sr(m)*sr(m)));
293 };
295 //---------------------------- bvav2nuv --------------------------------
296 // Compute lattice tap coefficients from transfer-function coefficients
297 //
298 // USAGE:
299 // nuv = bvav2nuv(bv,av)
300 // where
301 // av = parallel signal bank a1,...,aN
302 // bv = parallel signal bank b0,b1,...,aN
303 // nuv = parallel signal bank nu1,...,nuN
304 // where nui is the i'th tap coefficient,
305 // bi is the coefficient of z^(-i) in the filter numerator,
306 // ai is the coefficient of z^(-i) in the filter denominator
308 bvav2nuv(bv,av) = par(m,M+1,nu(m)) with {
309 M = count(av);
310 nu(m) = take(m+1,Pr(M-m)); // m=0..M
312 Pr(0) = bv; // Pr(m) is order M-m, 'r' means "reversed"
313 Pr(m) = par(i,M-m+1, (Pri(m,i) - nu(M-m+1)*Ari(m,M-m-i+1)));
314 Pri(m,i) = take(i+1,Pr(m-1));
315 Ari(m,i) = take(i+1,Ar(m-1));
316 //step-down recursion for lattice/ladder digital filters:
317 Ar(0) = (1,av); // Ar(m) is order M-m (recursion index must start at constant)
318 Ar(m) = 1,par(i,M-m, (Ari(m,i+1) - sr(m)*Ari(m,M-m-i))/(1-sr(m)*sr(m)));
319 sr(m) = Ari(m,M-m+1); // s_{M-1-m}
320 };
322 //--------------------------- iir_lat2, allpassnt -----------------------
324 iir_lat2(bv,av) = allpassnt(M,sv) : sum(i,M+1,*(take(M-i+1,tg)))
325 with {
326 M = count(av);
327 sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
328 tg = bvav2nuv(bv,av); // tg = vector of tap gains
329 };
331 // two-multiply lattice allpass (nested order-1 direct-form-ii allpasses):
332 allpassnt(0,sv) = _;
333 allpassnt(n,sv) =
334 //0: x <: ((+ <: (allpassnt(n-1,sv)),*(s))~(*(-s))) : _',_ :+
335 _ : ((+ <: (allpassnt(n-1,sv),*(s)))~*(-s)) : fsec(n)
336 with {
337 fsec(1) = crossnn(1) : _, (_<:mem,_) : +,_;
338 fsec(n) = crossn1(n) : _, (_<:mem,_),par(i,n-1,_) : +, par(i,n,_);
339 innertaps(n) = par(i,n,_);
340 s = take(n,sv); // reflection coefficient s = sin(theta)
341 };
343 //------------------------------- iir_kl, allpassnklt -------------------------
344 iir_kl(bv,av) = allpassnklt(M,sv) : sum(i,M+1,*(tghr(i)))
345 with {
346 M = count(av);
347 sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
348 tg = bvav2nuv(bv,av); // tg = vector of tap gains for 2mul case
349 tgr(i) = take(M+1-i,tg);
350 tghr(n) = tgr(n)/pi(n);
351 pi(0) = 1;
352 pi(n) = pi(n-1)*(1+take(M-n+1,sv)); // all sign parameters '+'
353 };
355 // Kelly-Lochbaum ladder allpass with tap lines:
356 allpassnklt(0,sv) = _;
357 allpassnklt(n,sv) = _ <: *(s),(*(1+s) : (+
358 : allpassnklt(n-1,sv))~(*(-s))) : fsec(n)
359 with {
360 fsec(1) = _, (_<:mem*(1-s),_) : sumandtaps(n);
361 fsec(n) = _, (_<:mem*(1-s),_), par(i,n-1,_) : sumandtaps(n);
362 s = take(n,sv);
363 sumandtaps(n) = +,par(i,n,_);
364 };
367 //------------------------------- iir_lat1, allpassn1mt -------------------------
368 iir_lat1(bv,av) = allpassn1mt(M,sv) : sum(i,M+1,*(tghr(i+1)))
369 with {
370 M = count(av);
371 sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
372 tg = bvav2nuv(bv,av); // tg = vector of tap gains
373 tgr(i) = take(M+2-i,tg); // i=1..M+1 (for "takability")
374 tghr(n)=tgr(n)/pi(n);
375 pi(1) = 1;
376 pi(n) = pi(n-1)*(1+take(M-n+2,sv)); // all sign parameters '+'
377 };
379 // one-multiply lattice allpass with tap lines:
380 allpassn1mt(0,sv) = _;
381 allpassn1mt(n,sv)= _ <: _,_ : ((+:*(s) <: _,_),_ : _,+ : crossnn(1)
382 : allpassn1mt(n-1,sv),_)~(*(-1)) : fsec(n)
383 with {
384 //0: fsec(n) = _',_ : +
385 fsec(1) = crossnn(1) : _, (_<:mem,_) : +,_;
386 fsec(n) = crossn1(n) : _, (_<:mem,_),par(i,n-1,_) : +, par(i,n,_);
387 innertaps(n) = par(i,n,_);
388 s = take(n,sv); // reflection coefficient s = sin(theta)
389 };
391 //------------------------------- iir_nl, allpassnnlt -------------------------
393 //
394 // REFERENCES:
395 // J. D. Markel and A. H. Gray, Linear Prediction of Speech, New York: Springer Verlag, 1976.
396 // https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html
398 iir_nl(bv,av) = allpassnnlt(M,sv) : sum(i,M+1,*(tghr(i)))
399 with {
400 M = count(av);
401 sv = av2sv(av); // sv = vector of sin(theta) reflection coefficients
402 tg = bvav2nuv(bv,av); // tg = vector of tap gains for 2mul case
403 tgr(i) = take(M+1-i,tg);
404 tghr(n) = tgr(n)/pi(n);
405 pi(0) = 1;
406 s(n) = take(M-n+1,sv);
407 c(n) = sqrt(max(0,1-s(n)*s(n))); // compiler crashes on sqrt(-)
408 pi(n) = pi(n-1)*c(n);
409 };
411 // Normalized ladder allpass with tap lines:
412 allpassnnlt(0,sv) = _;
413 allpassnnlt(n,scl*(sv)) = allpassnnlt(n,par(i,count(sv),scl*(sv(i))));
414 allpassnnlt(n,sv) = _ <: *(s),(*(c) : (+
415 : allpassnnlt(n-1,sv))~(*(-s))) : fsec(n)
416 with {
417 fsec(1) = _, (_<:mem*(c),_) : sumandtaps(n);
418 fsec(n) = _, (_<:mem*(c),_), par(i,n-1,_) : sumandtaps(n);
419 s = take(n,sv);
420 c = sqrt(max(0,1-s*s));
421 sumandtaps(n) = +,par(i,n,_);
422 };
424 //========================= Useful special cases ============================
426 //-------------------------------- tf2np ------------------------------------
427 // tf2np - biquad based on a stable second-order Normalized Ladder Filter
428 // (more robust to modulation than tf2 and protected against instability)
429 tf2np(b0,b1,b2,a1,a2) = allpassnnlt(M,sv) : sum(i,M+1,*(tghr(i)))
430 with {
431 smax = 0.9999; // maximum reflection-coefficient magnitude allowed
432 s2 = max(-smax, min(smax,a2)); // Project both reflection-coefficients
433 s1 = max(-smax, min(smax,a1/(1+a2))); // into the defined stability-region.
434 sv = (s1,s2); // vector of sin(theta) reflection coefficients
435 M = 2;
436 nu(2) = b2;
437 nu(1) = b1 - b2*a1;
438 nu(0) = (b0-b2*a2) - nu(1)*s1;
439 tg = (nu(0),nu(1),nu(2));
440 tgr(i) = take(M+1-i,tg); // vector of tap gains for 2mul case
441 tghr(n) = tgr(n)/pi(n); // apply pi parameters for NLF case
442 pi(0) = 1;
443 s(n) = take(M-n+1,sv);
444 c(n) = sqrt(1-s(n)*s(n));
445 pi(n) = pi(n-1)*c(n);
446 };
448 //----------------------------- wgr ---------------------------------
449 // Second-order transformer-normalized digital waveguide resonator
450 // USAGE:
451 // _ : wgr(f,r) : _
452 // where
453 // f = resonance frequency (Hz)
454 // r = loss factor for exponential decay
455 // (set to 1 to make a numerically stable oscillator)
456 //
457 // REFERENCES:
458 // https://ccrma.stanford.edu/~jos/pasp/Power_Normalized_Waveguide_Filters.html
459 // https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html
460 //
461 wgr(f,r,x) = (*(G),_<:_,((+:*(C))<:_,_),_:+,_,_:+(x),-) ~ cross : _,*(0-gi)
462 with {
463 C = cos(2*PI*f/SR);
464 gi = sqrt(max(0,(1+C)/(1-C))); // compensate amplitude (only needed when
465 G = r*(1-1' + gi')/gi; // frequency changes substantially)
466 cross = _,_ <: !,_,_,!;
467 };
469 //----------------------------- nlf2 --------------------------------
470 // Second order normalized digital waveguide resonator
471 // USAGE:
472 // _ : nlf2(f,r) : _
473 // where
474 // f = resonance frequency (Hz)
475 // r = loss factor for exponential decay
476 // (set to 1 to make a sinusoidal oscillator)
477 //
478 // REFERENCE:
479 // https://ccrma.stanford.edu/~jos/pasp/Power_Normalized_Waveguide_Filters.html
480 //
481 nlf2(f,r,x) = ((_<:_,_),(_<:_,_) : (*(s),*(c),*(c),*(0-s)) :>
482 (*(r),+(x))) ~ cross
483 with {
484 th = 2*PI*f/SR;
485 c = cos(th);
486 s = sin(th);
487 cross = _,_ <: !,_,_,!;
488 };
490 //===================== Ladder/Lattice Allpass Filters ======================
491 // An allpass filter has gain 1 at every frequency, but variable phase.
492 // Ladder/lattice allpass filters are specified by reflection coefficients.
493 // They are defined here as nested allpass filters, hence the names allpassn*.
494 //
495 // REFERENCES
497 // https://ccrma.stanford.edu/~jos/pasp/Nested_Allpass_Filters.html
498 // 2. Linear Prediction of Speech, Markel and Gray, Springer Verlag, 1976
499 //
500 // QUICK GUIDE
501 // allpassn - two-multiply lattice - each section is two multiply-adds
502 // allpassnn - normalized form - four multiplies and two adds per section,
503 // but coefficients can be time varying and nonlinear without
504 // "parametric amplification" (modulation of signal energy).
505 // allpassnkl - Kelly-Lochbaum form - four multiplies and two adds per
506 // section, but all signals have an immediate physical
507 // interpretation as traveling pressure waves, etc.
508 // allpassn1m - One-multiply form - one multiply and three adds per section.
509 // Normally the most efficient in special-purpose hardware.
510 //
511 // TYPICAL USAGE
512 // _ : allpassn(N,sv) : _
513 // where
514 // N = allpass order (number of ladder or lattice sections)
515 // sv = (s1,s2,...,sN) = reflection coefficients (between -1 and 1).
516 // For allpassnn only, sv is replaced by tv, where sv(i) = sin(tv(i)),
517 // where tv(i) may range between -PI and PI.
518 //
519 // two-multiply:
520 allpassn(0,sv) = _;
521 allpassn(n,sv) = _ <: ((+ <: (allpassn(n-1,sv)),*(s))~(*(-s))) : _',_ :+
522 with { s = take(n,sv); };
524 // power-normalized (reflection coefficients s = sin(t)):
525 allpassnn(0,tv) = _;
526 allpassnn(n,tv) = _ <: *(s), (*(c) : (+
527 : allpassnn(n-1,tv))~(*(-s))) : _, mem*c : +
528 with { c=cos(take(n,tv)); s=sin(take(n,tv)); };
530 // power-normalized with sparse delays dv(n)>1:
531 allpassnns(0,tv,dmax,dv) = _;
532 allpassnns(n,tv,dmax,dv) = _ <: *(s), (*(c) : (+ : dl
533 : allpassnns(n-1,tv,dmax,dv))~(*(-s))) : _, mem*c : +
534 with { c=cos(take(n,tv)); s=sin(take(n,tv));
535 dl=delay(dmax,(take(n,dv)-1)); };
537 // Kelly-Lochbaum:
538 allpassnkl(0,sv) = _;
539 allpassnkl(n,sv) = _ <: *(s),(*(1+s) : (+
540 : allpassnkl(n-1,sv))~(*(-s))) : _, mem*(1-s) : +
541 with { s = take(n,sv); };
543 // one-multiply:
544 allpassn1m(0,sv) = _;
545 allpassn1m(n,sv)= _ <: _,_ : ((+:*(s) <: _,_),_ : _,+ : cross
546 : allpassn1m(n-1,sv),_)~(*(-1)) : _',_ : +
547 with {s = take(n,sv); cross = _,_ <: !,_,_,!; };
549 //===== Digital Filter Sections Specified as Analog Filter Sections =====
550 //
551 //------------------------- tf2s, tf2snp --------------------------------
552 // Second-order direct-form digital filter,
553 // specified by ANALOG transfer-function polynomials B(s)/A(s),
554 // and a frequency-scaling parameter. Digitization via the
555 // bilinear transform is built in.
556 //
557 // USAGE: tf2s(b2,b1,b0,a1,a0,w1), where
558 //
559 // b2 s^2 + b1 s + b0
560 // H(s) = --------------------
561 // s^2 + a1 s + a0
562 //
563 // and w1 is the desired digital frequency (in radians/second)
564 // corresponding to analog frequency 1 rad/sec (i.e., s = j).
565 //
566 // EXAMPLE: A second-order ANALOG Butterworth lowpass filter,
567 // normalized to have cutoff frequency at 1 rad/sec,
568 // has transfer function
569 //
570 // 1
571 // H(s) = -----------------
572 // s^2 + a1 s + 1
573 //
574 // where a1 = sqrt(2). Therefore, a DIGITAL Butterworth lowpass
575 // cutting off at SR/4 is specified as tf2s(0,0,1,sqrt(2),1,PI*SR/2);
576 //
577 // METHOD: Bilinear transform scaled for exact mapping of w1.
578 // REFERENCE:
579 // https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html
580 //
581 tf2s(b2,b1,b0,a1,a0,w1) = tf2(b0d,b1d,b2d,a1d,a2d)
582 with {
583 c = 1/tan(w1*0.5/SR); // bilinear-transform scale-factor
584 csq = c*c;
585 d = a0 + a1 * c + csq;
586 b0d = (b0 + b1 * c + b2 * csq)/d;
587 b1d = 2 * (b0 - b2 * csq)/d;
588 b2d = (b0 - b1 * c + b2 * csq)/d;
589 a1d = 2 * (a0 - csq)/d;
590 a2d = (a0 - a1*c + csq)/d;
591 };
593 // tf2snp = tf2s but using a protected normalized ladder filter for tf2:
594 tf2snp(b2,b1,b0,a1,a0,w1) = tf2np(b0d,b1d,b2d,a1d,a2d)
595 with {
596 c = 1/tan(w1*0.5/SR); // bilinear-transform scale-factor
597 csq = c*c;
598 d = a0 + a1 * c + csq;
599 b0d = (b0 + b1 * c + b2 * csq)/d;
600 b1d = 2 * (b0 - b2 * csq)/d;
601 b2d = (b0 - b1 * c + b2 * csq)/d;
602 a1d = 2 * (a0 - csq)/d;
603 a2d = (a0 - a1*c + csq)/d;
604 };
606 //----------------------------- tf1s --------------------------------
607 // First-order direct-form digital filter,
608 // specified by ANALOG transfer-function polynomials B(s)/A(s),
609 // and a frequency-scaling parameter.
610 //
611 // USAGE: tf1s(b1,b0,a0,w1), where
612 //
613 // b1 s + b0
614 // H(s) = ----------
615 // s + a0
616 //
617 // and w1 is the desired digital frequency (in radians/second)
618 // corresponding to analog frequency 1 rad/sec (i.e., s = j).
619 //
620 // EXAMPLE: A first-order ANALOG Butterworth lowpass filter,
621 // normalized to have cutoff frequency at 1 rad/sec,
622 // has transfer function
623 //
624 // 1
625 // H(s) = -------
626 // s + 1
627 //
628 // so b0 = a0 = 1 and b1 = 0. Therefore, a DIGITAL first-order
629 // Butterworth lowpass with gain -3dB at SR/4 is specified as
630 //
631 // tf1s(0,1,1,PI*SR/2); // digital half-band order 1 Butterworth
632 //
633 // METHOD: Bilinear transform scaled for exact mapping of w1.
634 // REFERENCE:
635 // https://ccrma.stanford.edu/~jos/pasp/Bilinear_Transformation.html
636 //
637 tf1s(b1,b0,a0,w1) = tf1(b0d,b1d,a1d)
638 with {
639 c = 1/tan((w1)*0.5/SR); // bilinear-transform scale-factor
640 d = a0 + c;
641 b1d = (b0 - b1*c) / d;
642 b0d = (b0 + b1*c) / d;
643 a1d = (a0 - c) / d;
644 };
646 //----------------------------- tf2sb --------------------------------
647 // Bandpass mapping of tf2s: In addition to a frequency-scaling parameter
648 // w1 (set to HALF the desired passband width in rad/sec),
649 // there is a desired center-frequency parameter wc (also in rad/s).
650 // Thus, tf2sb implements a fourth-order digital bandpass filter section
651 // specified by the coefficients of a second-order analog lowpass prototpe
652 // section. Such sections can be combined in series for higher orders.
653 // The order of mappings is (1) frequency scaling (to set lowpass cutoff w1),
654 // (2) bandpass mapping to wc, then (3) the bilinear transform, with the
655 // usual scale parameter 2*SR. Algebra carried out in maxima and pasted here.
656 //
657 tf2sb(b2,b1,b0,a1,a0,w1,wc) =
658 iir((b0d/a0d,b1d/a0d,b2d/a0d,b3d/a0d,b4d/a0d),(a1d/a0d,a2d/a0d,a3d/a0d,a4d/a0d)) with {
659 T = 1.0/float(SR);
660 b0d = (4*b0*w1^2+8*b2*wc^2)*T^2+8*b1*w1*T+16*b2;
661 b1d = 4*b2*wc^4*T^4+4*b1*wc^2*w1*T^3-16*b1*w1*T-64*b2;
662 b2d = 6*b2*wc^4*T^4+(-8*b0*w1^2-16*b2*wc^2)*T^2+96*b2;
663 b3d = 4*b2*wc^4*T^4-4*b1*wc^2*w1*T^3+16*b1*w1*T-64*b2;
664 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)
665 + b2*wc^4*T^4+2*b1*wc^2*w1*T^3;
666 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;
667 a1d = 4*wc^4*T^4+4*a1*wc^2*w1*T^3-16*a1*w1*T-64;
668 a2d = 6*wc^4*T^4+(-8*a0*w1^2-16*wc^2)*T^2+96;
669 a3d = 4*wc^4*T^4-4*a1*wc^2*w1*T^3+16*a1*w1*T-64;
670 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;
671 };
673 //----------------------------- tf1sb --------------------------------
674 // First-to-second-order lowpass-to-bandpass section mapping,
675 // analogous to tf2sb above.
676 //
677 tf1sb(b1,b0,a0,w1,wc) = tf2(b0d/a0d,b1d/a0d,b2d/a0d,a1d/a0d,a2d/a0d) with {
678 T = 1.0/float(SR);
679 a0d = wc^2*T^2+2*a0*w1*T+4;
680 b0d = b1*wc^2*T^2 +2*b0*w1*T+4*b1;
681 b1d = 2*b1*wc^2*T^2-8*b1;
682 b2d = b1*wc^2*T^2-2*b0*w1*T+4*b1;
683 a1d = 2*wc^2*T^2-8;
684 a2d = wc^2*T^2-2*a0*w1*T+4;
685 };
687 //====================== Simple Resonator Filters ======================
689 // resonlp = 2nd-order lowpass with corner resonance:
690 resonlp(fc,Q,gain) = tf2s(b2,b1,b0,a1,a0,wc)
691 with {
692 wc = 2*PI*fc;
693 a1 = 2/Q;
694 a0 = 1;
695 b2 = 0;
696 b1 = 0;
697 b0 = gain;
698 };
700 // resonhp = 2nd-order highpass with corner resonance:
701 resonhp(fc,Q,gain,x) = gain*x-resonlp(fc,Q,gain,x);
703 // resonbp = 2nd-order bandpass
704 resonbp(fc,Q,gain) = tf2s(b2,b1,b0,a1,a0,wc)
705 with {
706 wc = 2*PI*fc;
707 a1 = 2/Q;
708 a0 = 1;
709 b2 = 0;
710 b1 = gain;
711 b0 = 0;
712 };
714 //================ Butterworth Lowpass/Highpass Filters ======================
715 // Nth-order Butterworth lowpass or highpass filters
716 //
717 // USAGE:
718 // _ : lowpass(N,fc) : _
719 // _ : highpass(N,fc) : _
720 // where
721 // N = filter order (number of poles) [nonnegative integer]
722 // fc = desired cut-off frequency (-3dB frequency) in Hz
723 // REFERENCE:
724 // https://ccrma.stanford.edu/~jos/filters/Butterworth_Lowpass_Design.html
725 // 'butter' function in Octave ("[z,p,g] = butter(N,1,'s');")
726 // ACKNOWLEDGMENT
727 // Generalized recursive formulation initiated by Yann Orlarey.
729 lowpass(N,fc) = lowpass0_highpass1(0,N,fc);
730 highpass(N,fc) = lowpass0_highpass1(1,N,fc);
731 lowpass0_highpass1(s,N,fc) = lphpr(s,N,N,fc)
732 with {
733 lphpr(s,0,N,fc) = _;
734 lphpr(s,1,N,fc) = tf1s(s,1-s,1,2*PI*fc);
735 lphpr(s,O,N,fc) = lphpr(s,(O-2),N,fc) : tf2s(s,0,1-s,a1s,1,w1) with {
736 parity = N % 2;
737 S = (O-parity)/2; // current section number
738 a1s = -2*cos(-PI + (1-parity)*PI/(2*N) + (S-1+parity)*PI/N);
739 w1 = 2*PI*fc;
740 };
741 };
743 //========== Special Filter-Bank Delay-Equalizing Allpass Filters ===========
744 //
745 // These special allpass filters are needed by filterbank et al. below.
746 // They are equivalent to (lowpass(N,fc) +|- highpass(N,fc))/2, but with
747 // canceling pole-zero pairs removed (which occurs for odd N).
749 //-------------------- lowpass_plus|minus_highpass ------------------
751 highpass_plus_lowpass(1,fc) = _;
752 highpass_plus_lowpass(3,fc) = tf2s(1,-1,1,1,1,w1) with { w1 = 2*PI*fc; };
753 highpass_minus_lowpass(3,fc) = tf1s(-1,1,1,w1) with { w1 = 2*PI*fc; };
754 highpass_plus_lowpass(5,fc) = tf2s(1,-a11,1,a11,1,w1)
755 with {
756 a11 = 1.618033988749895;
757 w1 = 2*PI*fc;
758 };
759 highpass_minus_lowpass(5,fc) = tf1s(1,-1,1,w1) : tf2s(1,-a12,1,a12,1,w1)
760 with {
761 a12 = 0.618033988749895;
762 w1 = 2*PI*fc;
763 };
765 // Catch-all definitions for generality - even order is done:
767 highpass_plus_lowpass(N,fc) = switch_odd_even(N%2,N,fc) with {
768 switch_odd_even(0,N,fc) = highpass_plus_lowpass_even(N,fc);
769 switch_odd_even(1,N,fc) = highpass_plus_lowpass_odd(N,fc);
770 };
772 highpass_minus_lowpass(N,fc) = switch_odd_even(N%2,N,fc) with {
773 switch_odd_even(0,N,fc) = highpass_minus_lowpass_even(N,fc);
774 switch_odd_even(1,N,fc) = highpass_minus_lowpass_odd(N,fc);
775 };
777 highpass_plus_lowpass_even(N,fc) = highpass(N,fc) + lowpass(N,fc);
778 highpass_minus_lowpass_even(N,fc) = highpass(N,fc) - lowpass(N,fc);
780 // FIXME: Rewrite the following, as for orders 3 and 5 above,
781 // to eliminate pole-zero cancellations:
782 highpass_plus_lowpass_odd(N,fc) = highpass(N,fc) + lowpass(N,fc);
783 highpass_minus_lowpass_odd(N,fc) = highpass(N,fc) - lowpass(N,fc);
785 //===================== Elliptic (Cauer) Lowpass Filters ===================
786 // USAGE:
787 // _ : lowpass3e(fc) : _
788 // _ : lowpass6e(fc) : _
789 // where fc = -3dB frequency in Hz
790 //
791 // REFERENCES:
792 // http://en.wikipedia.org/wiki/Elliptic_filter
793 // functions 'ncauer' and 'ellip' in Octave
795 //----------------------------- lowpass3e -----------------------------
796 // Third-order Elliptic (Cauer) lowpass filter
797 // DESIGN: For spectral band-slice level display (see octave_analyzer3e):
798 // [z,p,g] = ncauer(Rp,Rs,3); % analog zeros, poles, and gain, where
799 // Rp = 60 % dB ripple in stopband
800 // Rs = 0.2 % dB ripple in passband
801 //
802 lowpass3e(fc) = tf2s(b21,b11,b01,a11,a01,w1) : tf1s(0,1,a02,w1)
803 with {
804 a11 = 0.802636764161030; // format long; poly(p(1:2)) % in octave
805 a01 = 1.412270893774204;
806 a02 = 0.822445908998816; // poly(p(3)) % in octave
807 b21 = 0.019809144837789; // poly(z)
808 b11 = 0;
809 b01 = 1.161516418982696;
810 w1 = 2*PI*fc;
811 };
813 //----------------------------- lowpass6e -----------------------------
814 // Sixth-order Elliptic/Cauer lowpass filter
815 // DESIGN: For spectral band-slice level display (see octave_analyzer6e):
816 // [z,p,g] = ncauer(Rp,Rs,6); % analog zeros, poles, and gain, where
817 // Rp = 80 % dB ripple in stopband
818 // Rs = 0.2 % dB ripple in passband
819 //
820 lowpass6e(fc) =
821 tf2s(b21,b11,b01,a11,a01,w1) :
822 tf2s(b22,b12,b02,a12,a02,w1) :
823 tf2s(b23,b13,b03,a13,a03,w1)
824 with {
825 b21 = 0.000099999997055;
826 a21 = 1;
827 b11 = 0;
828 a11 = 0.782413046821645;
829 b01 = 0.000433227200555;
830 a01 = 0.245291508706160;
831 b22 = 1;
832 a22 = 1;
833 b12 = 0;
834 a12 = 0.512478641889141;
835 b02 = 7.621731298870603;
836 a02 = 0.689621364484675;
837 b23 = 1;
838 a23 = 1;
839 b13 = 0;
840 a13 = 0.168404871113589;
841 b03 = 53.536152954556727;
842 a03 = 1.069358407707312;
843 w1 = 2*PI*fc;
844 };
846 //===================== Elliptic Highpass Filters =====================
847 // USAGE:
848 // _ : highpass3e(fc) : _
849 // _ : highpass6e(fc) : _
850 // where fc = -3dB frequency in Hz
852 //----------------------------- highpass3e -----------------------------
853 // Third-order Elliptic (Cauer) highpass filter
854 // DESIGN: Inversion of lowpass3e wrt unit circle in s plane (s <- 1/s)
855 //
856 highpass3e(fc) = tf2s(b01/a01,b11/a01,b21/a01,a11/a01,1/a01,w1) :
857 tf1s(1/a02,0,1/a02,w1)
858 with {
859 a11 = 0.802636764161030;
860 a01 = 1.412270893774204;
861 a02 = 0.822445908998816;
862 b21 = 0.019809144837789;
863 b11 = 0;
864 b01 = 1.161516418982696;
865 w1 = 2*PI*fc;
866 };
868 //----------------------------- highpass6e -----------------------------
869 // Sixth-order Elliptic/Cauer highpass filter
870 // METHOD: Inversion of lowpass3e wrt unit circle in s plane (s <- 1/s)
871 //
872 highpass6e(fc) =
873 tf2s(b01/a01,b11/a01,b21/a01,a11/a01,1/a01,w1) :
874 tf2s(b02/a02,b12/a02,b22/a02,a12/a02,1/a02,w1) :
875 tf2s(b03/a03,b13/a03,b23/a03,a13/a03,1/a03,w1)
876 with {
877 b21 = 0.000099999997055;
878 a21 = 1;
879 b11 = 0;
880 a11 = 0.782413046821645;
881 b01 = 0.000433227200555;
882 a01 = 0.245291508706160;
883 b22 = 1;
884 a22 = 1;
885 b12 = 0;
886 a12 = 0.512478641889141;
887 b02 = 7.621731298870603;
888 a02 = 0.689621364484675;
889 b23 = 1;
890 a23 = 1;
891 b13 = 0;
892 a13 = 0.168404871113589;
893 b03 = 53.536152954556727;
894 a03 = 1.069358407707312;
895 w1 = 2*PI*fc;
896 };
898 //================== Butterworth Bandpass/Bandstop Filters =====================
899 // Order 2*Nh Butterworth bandpass filter made using the transformation
900 // s <- s + wc^2/s on lowpass(Nh), where wc is the desired bandpass center
901 // frequency. The lowpass(Nh) cutoff w1 is half the desired bandpass width.
902 // A notch-like "bandstop" filter is similarly made from highpass(Nh).
903 //
904 // USAGE:
905 // _ : bandpass(Nh,fl,fu) : _
906 // _ : bandstop(Nh,fl,fu) : _
907 // where
908 // Nh = HALF the desired bandpass/bandstop order (which is therefore even)
909 // fl = lower -3dB frequency in Hz
910 // fu = upper -3dB frequency in Hz
911 // Thus, the passband (stopband) width is fu-fl,
912 // and its center frequency is (fl+fu)/2.
913 //
914 // REFERENCE:
915 // http://cnx.org/content/m16913/latest/
916 //
917 bandpass(Nh,fl,fu) = bandpass0_bandstop1(0,Nh,fl,fu);
918 bandstop(Nh,fl,fu) = bandpass0_bandstop1(1,Nh,fl,fu);
919 bandpass0_bandstop1(s,Nh,fl,fu) = bpbsr(s,Nh,Nh,fl,fu)
920 with {
921 wl = 2*PI*fl; // digital (z-plane) lower passband edge
922 wu = 2*PI*fu; // digital (z-plane) upper passband edge
924 c = 2.0*SR; // bilinear transform scaling used in tf2sb, tf1sb
925 wla = c*tan(wl/c); // analog (s-splane) lower cutoff
926 wua = c*tan(wu/c); // analog (s-splane) upper cutoff
928 wc = sqrt(wla*wua); // s-plane center frequency
929 w1 = wua - wc^2/wua; // s-plane lowpass prototype cutoff
931 bpbsr(s,0,Nh,fl,fu) = _;
932 bpbsr(s,1,Nh,fl,fu) = tf1sb(s,1-s,1,w1,wc);
933 bpbsr(s,O,Nh,fl,fu) = bpbsr(s,O-2,Nh,fl,fu) : tf2sb(s,0,(1-s),a1s,1,w1,wc)
934 with {
935 parity = Nh % 2;
936 S = (O-parity)/2; // current section number
937 a1s = -2*cos(-PI + (1-parity)*PI/(2*Nh) + (S-1+parity)*PI/Nh);
938 };
939 };
941 //======================= Elliptic Bandpass Filters ============================
943 //----------------------------- bandpass6e -----------------------------
944 // Order 12 elliptic bandpass filter analogous to bandpass(6) above.
945 //
946 bandpass6e(fl,fu) = tf2sb(b21,b11,b01,a11,a01,w1,wc) : tf1sb(0,1,a02,w1,wc)
947 with {
948 a11 = 0.802636764161030; // In octave: format long; poly(p(1:2))
949 a01 = 1.412270893774204;
950 a02 = 0.822445908998816; // poly(p(3))
951 b21 = 0.019809144837789; // poly(z)
952 b11 = 0;
953 b01 = 1.161516418982696;
955 wl = 2*PI*fl; // digital (z-plane) lower passband edge
956 wu = 2*PI*fu; // digital (z-plane) upper passband edge
958 c = 2.0*SR; // bilinear transform scaling used in tf2sb, tf1sb
959 wla = c*tan(wl/c); // analog (s-splane) lower cutoff
960 wua = c*tan(wu/c); // analog (s-splane) upper cutoff
962 wc = sqrt(wla*wua); // s-plane center frequency
963 w1 = wua - wc^2/wua; // s-plane lowpass cutoff
964 };
966 //----------------------------- bandpass12e -----------------------------
968 bandpass12e(fl,fu) =
969 tf2sb(b21,b11,b01,a11,a01,w1,wc) :
970 tf2sb(b22,b12,b02,a12,a02,w1,wc) :
971 tf2sb(b23,b13,b03,a13,a03,w1,wc)
972 with { // octave script output:
973 b21 = 0.000099999997055;
974 a21 = 1;
975 b11 = 0;
976 a11 = 0.782413046821645;
977 b01 = 0.000433227200555;
978 a01 = 0.245291508706160;
979 b22 = 1;
980 a22 = 1;
981 b12 = 0;
982 a12 = 0.512478641889141;
983 b02 = 7.621731298870603;
984 a02 = 0.689621364484675;
985 b23 = 1;
986 a23 = 1;
987 b13 = 0;
988 a13 = 0.168404871113589;
989 b03 = 53.536152954556727;
990 a03 = 1.069358407707312;
992 wl = 2*PI*fl; // digital (z-plane) lower passband edge
993 wu = 2*PI*fu; // digital (z-plane) upper passband edge
995 c = 2.0*SR; // bilinear transform scaling used in tf2sb, tf1sb
996 wla = c*tan(wl/c); // analog (s-splane) lower cutoff
997 wua = c*tan(wu/c); // analog (s-splane) upper cutoff
999 wc = sqrt(wla*wua); // s-plane center frequency
1000 w1 = wua - wc^2/wua; // s-plane lowpass cutoff
1001 };
1003 //================= Parametric Equalizers (Shelf, Peaking) ==================
1004 // REFERENCES
1005 // - http://en.wikipedia.org/wiki/Equalization
1006 // - Digital Audio Signal Processing, Udo Zolzer, Wiley, 1999, p. 124
1007 // - http://www.harmony-central.com/Computer/Programming/Audio-EQ-Cookbook.txt
1008 // http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
1009 // - https://ccrma.stanford.edu/~jos/filters/Low_High_Shelving_Filters.html
1010 // - https://ccrma.stanford.edu/~jos/filters/Peaking_Equalizers.html
1011 // - maxmsp.lib in the Faust distribution
1012 // - bandfilter.dsp in the faust2pd distribution
1014 //----------------------------- low_shelf ------------------------------------
1015 // First-order "low shelf" filter (gain boost|cut between dc and some frequency)
1016 // USAGE: lowshelf(L0,fx), where
1017 // L0 = desired boost (dB) between dc and fx
1018 // fx = desired transition frequency (Hz) from boost to unity gain
1019 // The gain at SR/2 is constrained to be 1.
1020 //
1021 low_shelf = low_shelf3; // default
1022 low_shelf1(L0,fx,x) = x + (db2linear(L0)-1)*lowpass(1,fx,x);
1023 low_shelf1_l(G0,fx,x) = x + (G0-1)*lowpass(1,fx,x);
1024 low_shelf3(L0,fx,x) = x + (db2linear(L0)-1)*lowpass(3,fx,x);
1025 low_shelf5(L0,fx,x) = x + (db2linear(L0)-1)*lowpass(5,fx,x);
1027 //----------------------------- high_shelf -----------------------------------
1028 // First-order "high shelf" filter (gain boost|cut above some frequency)
1029 //
1030 // USAGE: high_shelf(Lpi,fx), where
1031 // Lpi = desired boost or cut (dB) between fx and SR/2
1032 // fx = desired transition frequency in Hz
1033 // The gain at dc is constrained to be 1
1034 //
1035 high_shelf=high_shelf3; //default
1036 high_shelf1(Lpi,fx,x) = x + (db2linear(Lpi)-1)*highpass(1,fx,x);
1037 high_shelf1_l(Gpi,fx,x) = x + (Gpi-1)*highpass(1,fx,x);
1038 high_shelf3(Lpi,fx,x) = x + (db2linear(Lpi)-1)*highpass(3,fx,x);
1039 high_shelf5(Lpi,fx,x) = x + (db2linear(Lpi)-1)*highpass(5,fx,x);
1041 //-------------------------------- peak_eq -----------------------------------
1042 // Second order "peaking equalizer" section
1043 // (gain boost or cut near some frequency)
1044 // Also called a "parametric equalizer" section
1045 // USAGE: _ : peak_eq(Lfx,fx,B) : _;
1046 // where
1047 // Lfx = level (dB) at fx
1048 // fx = peak frequency (Hz)
1049 // B = bandwidth (B) of peak in Hz
1050 //
1051 peak_eq(Lfx,fx,B) = tf2s(1,b1s,1,a1s,1,wx) with {
1052 T = float(1.0/SR);
1053 Bw = B*T/sin(wx*T); // prewarp s-bandwidth for more accuracy in z-plane
1054 a1 = PI*Bw;
1055 b1 = g*a1;
1056 g = db2linear(abs(Lfx));
1057 b1s = select2(Lfx>0,a1,b1); // When Lfx>0, pole dominates bandwidth
1058 a1s = select2(Lfx>0,b1,a1); // When Lfx<0, zero dominates
1059 wx = 2*PI*fx;
1060 };
1062 //------------------------------- peak_eq_cq ---------------------------------
1063 // Constant-Q second order peaking equalizer section
1064 // USAGE: _ : peak_eq_cq(Lfx,fx,Q) : _;
1065 // where
1066 // Lfx = level (dB) at fx
1067 // fx = boost or cut frequency (Hz)
1068 // Q = "Quality factor" = fx/B where B = bandwidth of peak in Hz
1069 //
1070 peak_eq_cq(Lfx,fx,Q) = peak_eq(Lfx,fx,fx/Q);
1072 //------------------------------- peak_eq_rm ---------------------------------
1073 // Regalia-Mitra second order peaking equalizer section
1074 // USAGE: _ : peak_eq_rm(Lfx,fx,tanPiBT) : _;
1075 // where
1076 // Lfx = level (dB) at fx
1077 // fx = boost or cut frequency (Hz)
1078 // tanPiBT = tan(PI*B/SR), where B = -3dB bandwidth (Hz) when 10^(Lfx/20) = 0
1079 // ~ PI*B/SR for narrow bandwidths B
1080 //
1081 // REFERENCE:
1082 // P.A. Regalia, S.K. Mitra, and P.P. Vaidyanathan,
1083 // "The Digital All-Pass Filter: A Versatile Signal Processing Building Block"
1084 // Proceedings of the IEEE, 76(1):19-37, Jan. 1988. (See pp. 29-30.)
1085 //
1086 peak_eq_rm(Lfx,fx,tanPiBT) = _ <: _,A,_ : +,- : *(0.5),*(K/2.0) : + with {
1087 A = tf2(k2, k1*(1+k2), 1, k1*(1+k2), k2) <: _,_; // allpass
1088 k1 = 0.0 - cos(2.0*PI*fx/SR);
1089 k2 = (1.0 - tanPiBT)/(1.0 + tanPiBT);
1090 K = db2linear(Lfx);
1091 };
1093 //-------------------------- parametric_eq_demo ------------------------------
1094 // USAGE: _ : parametric_eq_demo: _ ;
1095 parametric_eq_demo = // input_signal :
1096 low_shelf(LL,FL) :
1097 peak_eq(LP,FP,BP) :
1098 high_shelf(LH,FH)
1099 // Recommended:
1100 // : mth_octave_spectral_level_demo(2) // half-octave spectrum analyzer
1101 with {
1102 eq_group(x) = hgroup(" PARAMETRIC EQ SECTIONS
1103 [tooltip: See Faust's filter.lib for info and pointers]",x);
1105 ls_group(x) = eq_group(vgroup(" Low Shelf",x));
1106 LL = ls_group(hslider(" Low Boost|Cut [unit:dB] [style:knob]
1107 [tooltip: Amount of low-frequency boost or cut in decibels]",
1108 0,-40,40,0.1));
1109 FL = ls_group(hslider(" Transition Frequency [unit:Hz] [style:knob]
1110 [tooltip: Transition-frequency from boost (cut) to unity gain]",
1111 200,1,5000,1));
1113 pq_group(x) = eq_group(vgroup(" Peaking Equalizer
1114 [tooltip: Parametric Equalizer sections from filter.lib]",x));
1115 LP = pq_group(hslider(" Peak Boost|Cut [unit:dB] [style:knob]
1116 [tooltip: Amount of local boost or cut in decibels]",
1117 0,-40,40,0.1));
1118 FP = pq_group(hslider(" Peak Frequency [unit:PK] [style:knob]
1119 [tooltip: Peak Frequency in Piano Key (PK) units (A-440= 49 PK)]",
1120 49,1,100,1)) : smooth(0.999) : pianokey2hz
1121 with { pianokey2hz(x) = 440.0*pow(2.0, (x-49.0)/12); };
1123 Q = pq_group(hslider(" Peak Q [style:knob]
1124 [tooltip: Quality factor (Q) of the peak = center-frequency/bandwidth]",
1125 40,1,50,0.1));
1127 BP = FP/Q;
1129 hs_group(x) = eq_group(vgroup(" High Shelf
1130 [tooltip: A high shelf provides a boost or cut
1131 above some frequency]",x));
1132 LH = hs_group(hslider(" High Boost|Cut [unit:dB] [style:knob]
1133 [tooltip: Amount of high-frequency boost or cut in decibels]",
1134 0,-40,40,.1));
1135 FH = hs_group(hslider(" Transition Frequency [unit:Hz] [style:knob]
1136 [tooltip: Transition-frequency from boost (cut) to unity gain]",
1137 8000,20,10000,1));
1138 };
1140 //========================= Lagrange Interpolation ========================
1141 // Reference:
1142 // https://ccrma.stanford.edu/~jos/pasp/Lagrange_Interpolation.html
1143 //
1144 //------------------ fdelay1, fdelay2, fdelay3, fdelay4 ---------------
1145 // Delay lines interpolated using Lagrange interpolation
1146 // USAGE: _ : fdelayN(maxdelay, delay, inputsignal) : _
1147 // (exactly like fdelay in music.lib)
1148 // where N=1,2,3, or 4 is the order of the Lagrange interpolation polynomial.
1149 //
1150 // NOTE: requested delay should not be less than (N-1)/2.
1151 //
1152 // NOTE: While the implementations below appear to use multiple delay lines,
1153 // they in fact use only one thanks to optimization by the Faust compiler.
1155 // first-order case (linear interpolation) - equivalent to fdelay in music.lib
1156 // delay d in [0,1]
1157 fdelay1(n,d,x) = delay(n,id,x)*(1 - fd) + delay(n,id+1,x)*fd
1158 with {
1159 id = int(d);
1160 fd = frac(d);
1161 };
1163 // second-order (quadratic) case, delay in [0.5,1.5]
1164 // delay d should be at least 0.5
1165 fdelay2(n,d,x) = delay(n,id,x)*(1-fd)*(2-fd)/2
1166 + delay(n,id+1,x)*(2-fd)*fd
1167 + delay(n,id+2,x)*(fd-1)*fd/2
1168 with {
1169 o = 0.49999; // offset to make life easy for interpolator
1170 dmo = d - o; // assumed nonnegative
1171 id = int(dmo);
1172 fd = o + frac(dmo);
1173 };
1175 // third-order (cubic) case, delay in [1,2]
1176 // delay d should be at least 1
1177 fdelay3(n,d,x) = delay(n,id,x) * (0-fdm1*fdm2*fdm3)/6
1178 + delay(n,id+1,x) * fd*fdm2*fdm3/2
1179 + delay(n,id+2,x) * (0-fd*fdm1*fdm3)/2
1180 + delay(n,id+3,x) * fd*fdm1*fdm2/6
1181 with {
1182 id = int(d-1);
1183 fd = 1+frac(d);
1184 fdm1 = fd-1;
1185 fdm2 = fd-2;
1186 fdm3 = fd-3;
1187 };
1189 // fourth-order (quartic) case, delay in [1.5,2.5]
1190 // delay d should be at least 1.5
1191 fdelay4(n,d,x) = delay(n,id,x) * fdm1*fdm2*fdm3*fdm4/24
1192 + delay(n,id+1,x) * (0-fd*fdm2*fdm3*fdm4)/6
1193 + delay(n,id+2,x) * fd*fdm1*fdm3*fdm4/4
1194 + delay(n,id+3,x) * (0-fd*fdm1*fdm2*fdm4)/6
1195 + delay(n,id+4,x) * fd*fdm1*fdm2*fdm3/24
1196 with {
1197 //v1: o = 1;
1198 o = 1.49999;
1199 dmo = d - o; // assumed nonnegative
1200 id = int(dmo);
1201 fd = o + frac(dmo);
1202 fdm1 = fd-1;
1203 fdm2 = fd-2;
1204 fdm3 = fd-3;
1205 fdm4 = fd-4;
1206 };
1208 // fifth-order case, delay in [2,3]
1209 // delay d should be at least 2
1210 fdelay5(n,d,x) =
1211 delay(n,id,x) * -fdm1*fdm2*fdm3*fdm4*fdm5/120
1212 + delay(n,id+1,x) * fd* fdm2*fdm3*fdm4*fdm5/24
1213 + delay(n,id+2,x) * -fd*fdm1* fdm3*fdm4*fdm5/12
1214 + delay(n,id+3,x) * fd*fdm1*fdm2* fdm4*fdm5/12
1215 + delay(n,id+4,x) * -fd*fdm1*fdm2*fdm3* fdm5/24
1216 + delay(n,id+5,x) * fd*fdm1*fdm2*fdm3*fdm4 /120
1217 with {
1218 //v1: o = 1;
1219 o = 1.99999;
1220 dmo = d - o; // assumed nonnegative
1221 id = int(dmo);
1222 fd = o + frac(dmo);
1223 fdm1 = fd-1;
1224 fdm2 = fd-2;
1225 fdm3 = fd-3;
1226 fdm4 = fd-4;
1227 fdm5 = fd-5;
1228 };
1230 //====================== Thiran Allpass Interpolation =====================
1231 // Reference:
1232 // https://ccrma.stanford.edu/~jos/pasp/Thiran_Allpass_Interpolators.html
1233 //
1234 //---------------- fdelay1a, fdelay2a, fdelay3a, fdelay4a -------------
1235 // Delay lines interpolated using Thiran allpass interpolation
1236 // USAGE: fdelayNa(maxdelay, delay, inputsignal)
1237 // (exactly like fdelay in music.lib)
1238 // where N=1,2,3, or 4 is the order of the Thiran interpolation filter,
1239 // and the delay argument is at least N - 1/2.
1240 //
1241 // (Move the following and similar notes above to filter-lib-doc.txt?)
1242 //
1243 // NOTE: The interpolated delay should not be less than N - 1/2.
1244 // (The allpass delay ranges from N - 1/2 to N + 1/2.)
1245 // This constraint can be alleviated by altering the code,
1246 // but be aware that allpass filters approach zero delay
1247 // by means of pole-zero cancellations.
1248 // The delay range [N-1/2,N+1/2] is not optimal. What is?
1249 //
1250 // NOTE: Delay arguments too small will produce an UNSTABLE allpass!
1251 //
1252 // NOTE: Because allpass interpolation is recursive, it is not as robust
1253 // as Lagrange interpolation under time-varying conditions.
1254 // (You may hear clicks when changing the delay rapidly.)
1255 //
1256 // first-order allpass interpolation, delay d in [0.5,1.5]
1257 fdelay1a(n,d,x) = delay(n,id,x) : tf1(eta,1,eta)
1258 with {
1259 o = 0.49999; // offset to make life easy for allpass
1260 dmo = d - o; // assumed nonnegative
1261 id = int(dmo);
1262 fd = o + frac(dmo);
1263 eta = (1-fd)/(1+fd); // allpass coefficient
1264 };
1266 // second-order allpass delay in [1.5,2.5]
1267 fdelay2a(n,d,x) = delay(n,id,x) : tf2(a2,a1,1,a1,a2)
1268 with {
1269 o = 1.49999;
1270 dmo = d - o; // delay range is [order-1/2, order+1/2]
1271 id = int(dmo);
1272 fd = o + frac(dmo);
1273 a1o2 = (2-fd)/(1+fd); // share some terms (the compiler does this anyway)
1274 a1 = 2*a1o2;
1275 a2 = a1o2*(1-fd)/(2+fd);
1276 };
1278 // third-order allpass delay in [2.5,3.5]
1279 // delay d should be at least 2.5
1280 fdelay3a(n,d,x) = delay(n,id,x) : iir((a3,a2,a1,1),(a1,a2,a3))
1281 with {
1282 o = 2.49999;
1283 dmo = d - o;
1284 id = int(dmo);
1285 fd = o + frac(dmo);
1286 a1o3 = (3-fd)/(1+fd);
1287 a2o3 = a1o3*(2-fd)/(2+fd);
1288 a1 = 3*a1o3;
1289 a2 = 3*a2o3;
1290 a3 = a2o3*(1-fd)/(3+fd);
1291 };
1293 // fourth-order allpass delay in [3.5,4.5]
1294 // delay d should be at least 3.5
1295 fdelay4a(n,d,x) = delay(n,id,x) : tf4(a4,a3,a2,a1,1,a1,a2,a3,a4)
1296 with {
1297 o = 3.49999;
1298 dmo = d - o;
1299 id = int(dmo);
1300 fd = o + frac(dmo);
1301 a1o4 = (4-fd)/(1+fd);
1302 a2o6 = a1o4*(3-fd)/(2+fd);
1303 a3o4 = a2o6*(2-fd)/(3+fd);
1304 a1 = 4*a1o4;
1305 a2 = 6*a2o6;
1306 a3 = 4*a3o4;
1307 a4 = a3o4*(1-fd)/(4+fd);
1308 };
1310 //================ Mth-Octave Filter-Banks and Spectrum-Analyzers ============
1311 // Mth-octave filter-banks and spectrum-analyzers split the input signal into a
1312 // bank of parallel signals, one for each spectral band. The parameters are
1313 //
1314 // M = number of band-slices per octave (>1)
1315 // N = total number of bands (>2)
1316 // ftop = upper bandlimit of the Mth-octave bands (<SR/2)
1317 //
1318 // In addition to the Mth-octave output signals, there is a highpass signal
1319 // containing frequencies from ftop to SR/2, and a "dc band" lowpass signal
1320 // containing frequencies from 0 (dc) up to the start of the Mth-octave bands.
1321 // Thus, the N output signals are
1322 //
1323 // highpass(ftop), MthOctaveBands(M,N-2,ftop), dcBand(ftop*2^(-M*(N-1)))
1324 //
1325 // A FILTER-BANK is defined here as a signal bandsplitter having the
1326 // property that summing its output signals gives an allpass-filtered
1327 // version of the filter-bank input signal. A more conventional term for
1328 // this is an "allpass-complementary filter bank". If the allpass filter
1329 // is a pure delay (and possible scaling), the filter bank is said to be
1330 // a "perfect-reconstruction filter bank" (see Vaidyanathan-1993 cited
1331 // below for details). A "graphic equalizer", in which band signals
1332 // are scaled by gains and summed, should be based on a filter bank.
1333 //
1334 // A SPECTRUM-ANALYZER is defined here as any band-split whose bands span
1335 // the relevant spectrum, but whose band-signals do not
1336 // necessarily sum to the original signal, either exactly or to within an
1337 // allpass filtering. Spectrum analyzer outputs are normally at least nearly
1338 // "power complementary", i.e., the power spectra of the individual bands
1339 // sum to the original power spectrum (to within some negligible tolerance).
1340 //
1341 // The filter-banks below are implemented as Butterworth or Elliptic
1342 // spectrum-analyzers followed by delay equalizers that make them
1343 // allpass-complementary.
1344 //
1345 // INCREASING CHANNEL ISOLATION
1346 // Go to higher filter orders - see Regalia et al. or Vaidyanathan (cited
1347 // below) regarding the construction of more aggressive recursive
1348 // filter-banks using elliptic or Chebyshev prototype filters.
1349 //
1350 // REFERENCES
1351 // - "Tree-structured complementary filter banks using all-pass sections",
1352 // Regalia et al., IEEE Trans. Circuits & Systems, CAS-34:1470-1484, Dec. 1987
1353 // - "Multirate Systems and Filter Banks", P. Vaidyanathan, Prentice-Hall, 1993
1354 // - Elementary filter theory: https://ccrma.stanford.edu/~jos/filters/
1355 //
1356 //------------------------- mth_octave_analyzer ----------------------------
1357 //
1358 // USAGE
1359 // _ : mth_octave_analyzer(O,M,ftop,N) : par(i,N,_); // Oth-order Butterworth
1360 // _ : mth_octave_analyzer6e(M,ftop,N) : par(i,N,_); // 6th-order elliptic
1361 //
1362 // where
1363 // O = order of filter used to split each frequency band into two
1364 // M = number of band-slices per octave
1365 // ftop = highest band-split crossover frequency (e.g., 20 kHz)
1366 // N = total number of bands (including dc and Nyquist)
1367 //
1368 // ACKNOWLEDGMENT
1369 // Recursive band-splitting formulation improved by Yann Orlarey.
1371 mth_octave_analyzer6e(M,ftop,N) = _ <: bsplit(N-1) with {
1372 fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
1373 lp(n) = lowpass6e(fc(n)); // 6th-order elliptic - see other choices above
1374 hp(n) = highpass6e(fc(n)); // (search for lowpass* and highpass*)
1375 bsplit(0) = _;
1376 bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
1377 };
1379 // Butterworth analyzers may be cascaded with allpass
1380 // delay-equalizers to make (allpass-complementary) filter banks:
1382 mth_octave_analyzer(O,M,ftop,N) = _ <: bsplit(N-1) with {
1383 fc(n) = ftop * 2^(float(n-N+1)/float(M));
1384 lp(n) = lowpass(O,fc(n)); // Order O Butterworth
1385 hp(n) = highpass(O,fc(n));
1386 bsplit(0) = _;
1387 bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
1388 };
1390 mth_octave_analyzer3(M,ftop,N) = mth_octave_analyzer(3,M,ftop,N);
1391 mth_octave_analyzer5(M,ftop,N) = mth_octave_analyzer(5,M,ftop,N);
1392 mth_octave_analyzer_default = mth_octave_analyzer6e; // default analyzer
1394 //------------------------ mth_octave_filterbank -------------------------
1395 // Allpass-complementary filter banks based on Butterworth band-splitting.
1396 // For Butterworth band-splits, the needed delay equalizer is easily found.
1398 mth_octave_filterbank(O,M,ftop,N) =
1399 mth_octave_analyzer(O,M,ftop,N) :
1400 delayeq(N) with {
1401 fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
1402 ap(n) = highpass_plus_lowpass(O,fc(n)); // delay-equalizing allpass
1403 delayeq(N) = par(i,N-2,apchain(i+1)), _, _;
1404 apchain(i) = seq(j,N-1-i,ap(j+1));
1405 };
1407 // dc-inverted version. This reduces the delay-equalizer order for odd O.
1408 // Negating the input signal makes the dc band noninverting
1409 // and all higher bands sign-inverted (if preferred).
1410 mth_octave_filterbank_alt(O,M,ftop,N) =
1411 mth_octave_analyzer(O,M,ftop,N) : delayeqi(O,N) with {
1412 fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
1413 ap(n) = highpass_minus_lowpass(O,fc(n)); // half the order of 'plus' case
1414 delayeqi(N) = par(i,N-2,apchain(i+1)), _, *(-1.0);
1415 apchain(i) = seq(j,N-1-i,ap(j+1));
1416 };
1418 // Note that even-order cases require complex coefficients.
1419 // See Vaidyanathan 1993 and papers cited there for more info.
1421 mth_octave_filterbank3(M,ftop,N) = mth_octave_filterbank_alt(3,M,ftop,N);
1422 mth_octave_filterbank5(M,ftop,N) = mth_octave_filterbank(5,M,ftop,N);
1424 mth_octave_filterbank_default = mth_octave_filterbank5;
1426 //======================= Mth-Octave Spectral Level =========================
1427 // Spectral Level: Display (in bar graphs) the average signal level in each
1428 // spectral band.
1429 //
1430 //------------------------ mth_octave_spectral_level -------------------------
1431 // USAGE: _ : mth_octave_spectral_level(M,ftop,NBands,tau,dB_offset);
1432 // where
1433 // M = bands per octave
1434 // ftop = lower edge frequency of top band
1435 // NBands = number of passbands (including highpass and dc bands),
1436 // tau = spectral display averaging-time (time constant) in seconds,
1437 // dB_offset = constant dB offset in all band level meters.
1438 //
1439 mth_octave_spectral_level6e(M,ftop,N,tau,dB_offset) = _<:
1440 _,mth_octave_analyzer6e(M,ftop,N) :
1441 _,(display:>_):attach with {
1442 display = par(i,N,dbmeter(i));
1443 dbmeter(i) = abs : smooth(tau2pole(tau)) : linear2db : +(dB_offset) :
1444 meter(N-i-1);
1445 meter(i) = speclevel_group(vbargraph("[%2i] [unit:dB]
1446 [tooltip: Spectral Band Level in dB]", -50, 10));
1447 // Can M be included in the label string somehow?
1448 speclevel_group(x) = hgroup(" CONSTANT-Q SPECTRUM ANALYZER (6E)
1449 [tooltip: See Faust's filter.lib for documentation and references]", x);
1450 };
1452 mth_octave_spectral_level_default = mth_octave_spectral_level6e;
1453 spectral_level = mth_octave_spectral_level(2,10000,20); // simplest case
1455 //---------------------- mth_octave_spectral_level_demo ----------------------
1456 // Demonstrate mth_octave_spectral_level in a standalone GUI.
1457 //
1458 // USAGE: _ : mth_octave_spectral_level_demo(BandsPerOctave);
1460 mth_octave_spectral_level_demo(M) =
1461 mth_octave_spectral_level_default(M,ftop,N,tau,dB_offset)
1462 with {
1463 // Span nearly 10 octaves so that lowest band-edge is at
1464 // ftop*2^(-Noct+2) = 40 Hz when ftop=10 kHz:
1465 N = int(10*M); // without 'int()', segmentation fault observed for M=1.67
1466 ftop = 10000;
1467 ctl_group(x) = hgroup(" SPECTRUM ANALYZER CONTROLS", x);
1468 tau = ctl_group(hslider(" Level Averaging Time [unit:sec]
1469 [tooltip: band-level averaging time in seconds]",
1470 0.1,0,1,0.01));
1471 dB_offset = ctl_group(hslider(" Level dB Offset [unit:dB]
1472 [tooltip: Level offset in decibels]",
1473 50,0,100,1));
1474 };
1476 spectral_level_demo = mth_octave_spectral_level_demo(1.5); // 2/3 octave
1478 //---------------- (third|half)_octave_(analyzer|filterbank) -----------------
1480 // Named special cases of mth_octave_* with defaults filled in:
1482 third_octave_analyzer(N) = mth_octave_analyzer_default(3,10000,N);
1483 third_octave_filterbank(N) = mth_octave_filterbank_default(3,10000,N);
1484 // Third-Octave Filter-Banks have been used in audio for over a century.
1485 // See, e.g.,
1486 // Acoustics [the book], by L. L. Beranek
1487 // Amer. Inst. Physics for the Acoustical Soc. America,
1488 // http://asa.aip.org/publications.html, 1986 (1st ed.1954)
1490 // Third-octave bands across the audio spectrum are too wide for current
1491 // typical computer screens, so half-octave bands are the default:
1492 half_octave_analyzer(N) = mth_octave_analyzer_default(2,10000,N);
1493 half_octave_filterbank(N) = mth_octave_filterbank_default(2,10000,N);
1495 octave_filterbank(N) = mth_octave_filterbank_default(1,10000,N);
1496 octave_analyzer(N) = mth_octave_analyzer_default(1,10000,N);
1498 //=========================== Filter-Bank Demos ==============================
1499 // Graphic Equalizer: Each filter-bank output signal routes through a fader.
1500 //
1501 // USAGE: _ : mth_octave_filterbank_demo(M) : _
1502 // where
1503 // M = number of bands per octave
1505 mth_octave_filterbank_demo(M) = bp1(bp,mthoctavefilterbankdemo) with {
1506 bp1 = component("effect.lib").bypass1;
1507 mofb_group(x) = vgroup("CONSTANT-Q FILTER BANK (Butterworth dyadic tree)
1508 [tooltip: See Faust's filter.lib for documentation and references]", x);
1509 bypass_group(x) = mofb_group(hgroup("", x));
1510 slider_group(x) = mofb_group(hgroup("", x));
1511 N = 10*M; // total number of bands (highpass band, octave-bands, dc band)
1512 ftop = 10000;
1513 mthoctavefilterbankdemo = chan;
1514 chan = mth_octave_filterbank_default(M,ftop,N) :
1517 [tooltip: Bandpass filter gain in dB]", -10, -70, 10, 0.1)) :
1518 smooth(0.999);
1519 bp = bypass_group(checkbox(" Bypass
1520 [tooltip: When this is checked, the filter-bank has no effect]"));
1521 };
1523 filterbank_demo = mth_octave_filterbank_demo(1); // octave-bands = default
1525 //=========== Arbritary-Crossover Filter-Banks and Spectrum Analyzers ========
1526 // These are similar to the Mth-octave filter-banks above, except that the
1527 // band-split frequencies are passed explicitly as arguments.
1528 //
1529 // USAGE:
1530 // _ : filterbank (O,freqs) : par(i,N,_); // Butterworth band-splits
1531 // _ : filterbanki(O,freqs) : par(i,N,_); // Inverted-dc version
1532 // _ : analyzer (O,freqs) : par(i,N,_); // No delay equalizer
1533 //
1534 // where
1535 // O = band-split filter order (ODD integer required for filterbank[i])
1536 // freqs = (fc1,fc2,...,fcNs) [in numerically ascending order], where
1537 // Ns=N-1 is the number of octave band-splits
1538 // (total number of bands N=Ns+1).
1539 //
1540 // If frequencies are listed explicitly as arguments, enclose them in parens:
1541 //
1542 // _ : filterbank(3,(fc1,fc2)) : _,_,_
1543 //
1544 // ACKNOWLEDGMENT
1545 // Technique for processing a variable number of signal arguments due
1546 // to Yann Orlarey (as is the entire Faust framework!)
1547 //
1548 //------------------------------ analyzer --------------------------------------
1549 analyzer(O,lfreqs) = _ <: bsplit(nb) with
1550 {
1551 nb = count(lfreqs);
1552 fc(n) = take(n, lfreqs);
1553 lp(n) = lowpass(O,fc(n));
1554 hp(n) = highpass(O,fc(n));
1555 bsplit(0) = _;
1556 bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
1557 };
1559 //----------------------------- filterbank -------------------------------------
1560 filterbank(O,lfreqs) = analyzer(O,lfreqs) : delayeq with
1561 {
1562 nb = count(lfreqs);
1563 fc(n) = take(n, lfreqs);
1564 ap(n) = highpass_plus_lowpass(O,fc(n));
1565 delayeq = par(i,nb-1,apchain(nb-1-i)),_,_;
1566 apchain(0) = _;
1567 apchain(i) = ap(i) : apchain(i-1);
1568 };
1570 //----------------------------- filterbanki ------------------------------------
1571 filterbanki(O,lfreqs) = _ <: bsplit(nb) with
1572 {
1573 fc(n) = take(n, lfreqs);
1574 lp(n) = lowpass(O,fc(n));
1575 hp(n) = highpass(O,fc(n));
1576 ap(n) = highpass_minus_lowpass(O,fc(n));
1577 bsplit(0) = *(-1.0);
1578 bsplit(i) = (hp(i) : delayeq(i-1)), (lp(i) <: bsplit(i-1));
1579 delayeq(0) = _; // moving the *(-1) here inverts all outputs BUT dc
1580 delayeq(i) = ap(i) : delayeq(i-1);
1581 };