idsp/iir/
coefficients.rs

1use num_traits::{AsPrimitive, Float, FloatConst};
2use serde::{Deserialize, Serialize};
3
4/// Transition/corner shape
5#[derive(Copy, Clone, Debug, PartialEq, PartialOrd, Serialize, Deserialize)]
6pub enum Shape<T> {
7    /// Q, 1/sqrt(2) for critical
8    Q(T),
9    /// Relative bandwidth in octaves
10    Bandwidth(T),
11    /// Slope steepnes, 1 for critical
12    Slope(T),
13}
14
15impl<T: Float + FloatConst> Default for Shape<T> {
16    fn default() -> Self {
17        Self::Q(T::SQRT_2().recip())
18    }
19}
20
21/// Standard audio biquad filter builder
22///
23/// <https://www.w3.org/TR/audio-eq-cookbook/>
24#[derive(Copy, Clone, Debug, PartialEq, PartialOrd, Serialize, Deserialize)]
25pub struct Filter<T> {
26    /// Angular critical frequency (in units of sampling angular frequency)
27    /// Corner frequency, or 3dB cutoff frequency,
28    frequency: T,
29    /// Passband gain
30    gain: T,
31    /// Shelf gain (only for peaking, lowshelf, highshelf)
32    /// Relative to passband gain
33    shelf: T,
34    /// Transition/corner shape
35    shape: Shape<T>,
36}
37
38impl<T: Float + FloatConst> Default for Filter<T> {
39    fn default() -> Self {
40        Self {
41            frequency: T::zero(),
42            gain: T::one(),
43            shape: Shape::default(),
44            shelf: T::one(),
45        }
46    }
47}
48
49impl<T> Filter<T>
50where
51    T: 'static + Float + FloatConst,
52    f32: AsPrimitive<T>,
53{
54    /// Set crititcal frequency from absolute units.
55    ///
56    /// # Arguments
57    /// * `critical_frequency`: "Relevant" or "corner" or "center" frequency
58    ///   in the same units as `sample_frequency`
59    /// * `sample_frequency`: The sample frequency in the same units as `critical_frequency`.
60    ///   E.g. both in SI Hertz or `rad/s`.
61    pub fn frequency(&mut self, critical_frequency: T, sample_frequency: T) -> &mut Self {
62        self.critical_frequency(critical_frequency / sample_frequency)
63    }
64
65    /// Set relative critical frequency
66    ///
67    /// # Arguments
68    /// * `f0`: Relative critical frequency in units of the sample frequency.
69    ///   Must be `0 <= f0 <= 0.5`.
70    pub fn critical_frequency(&mut self, f0: T) -> &mut Self {
71        self.angular_critical_frequency(T::TAU() * f0)
72    }
73
74    /// Set relative critical angular frequency
75    ///
76    /// # Arguments
77    /// * `w0`: Relative critical angular frequency.
78    ///   Must be `0 <= w0 <= π`. Defaults to `0.0`.
79    pub fn angular_critical_frequency(&mut self, w0: T) -> &mut Self {
80        self.frequency = w0;
81        self
82    }
83
84    /// Set reference gain
85    ///
86    /// # Arguments
87    /// * `k`: Linear reference gain. Defaults to `1.0`.
88    pub fn gain(&mut self, k: T) -> &mut Self {
89        self.gain = k;
90        self
91    }
92
93    /// Set reference gain in dB
94    ///
95    /// # Arguments
96    /// * `k_db`: Reference gain in dB. Defaults to `0.0`.
97    pub fn gain_db(&mut self, k_db: T) -> &mut Self {
98        self.gain(10.0.as_().powf(k_db / 20.0.as_()))
99    }
100
101    /// Set linear shelf gain
102    ///
103    /// Used only for `peaking`, `highshelf`, `lowshelf` filters.
104    ///
105    /// # Arguments
106    /// * `a`: Linear shelf gain. Defaults to `1.0`.
107    pub fn shelf(&mut self, a: T) -> &mut Self {
108        self.shelf = a;
109        self
110    }
111
112    /// Set shelf gain in dB
113    ///
114    /// Used only for `peaking`, `highshelf`, `lowshelf` filters.
115    ///
116    /// # Arguments
117    /// * `a_db`: Linear shelf gain. Defaults to `0.0`.
118    pub fn shelf_db(&mut self, a_db: T) -> &mut Self {
119        self.shelf(10.0.as_().powf(a_db / 20.0.as_()))
120    }
121
122    /// Set inverse Q parameter of the filter
123    ///
124    /// The inverse "steepness"/"narrowness" of the filter transition.
125    /// Defaults `sqrt(2)` which is as steep as possible without overshoot.
126    ///
127    /// # Arguments
128    /// * `qi`: Inverse Q parameter.
129    pub fn inverse_q(&mut self, qi: T) -> &mut Self {
130        self.q(qi.recip())
131    }
132
133    /// Set Q parameter of the filter
134    ///
135    /// The "steepness"/"narrowness" of the filter transition.
136    /// Defaults `1/sqrt(2)` which is as steep as possible without overshoot.
137    ///
138    /// This affects the same parameter as `bandwidth()` and `shelf_slope()`.
139    /// Use only one of them.
140    ///
141    /// # Arguments
142    /// * `q`: Q parameter.
143    pub fn q(&mut self, q: T) -> &mut Self {
144        self.shape = Shape::Q(q);
145        self
146    }
147
148    /// Set the relative bandwidth
149    ///
150    /// This affects the same parameter as `inverse_q()` and `shelf_slope()`.
151    /// Use only one of them.
152    ///
153    /// # Arguments
154    /// * `bw`: Bandwidth in octaves
155    pub fn bandwidth(&mut self, bw: T) -> &mut Self {
156        self.shape = Shape::Bandwidth(bw);
157        self
158    }
159
160    /// Set the shelf slope.
161    ///
162    /// This affects the same parameter as `inverse_q()` and `bandwidth()`.
163    /// Use only one of them.
164    ///
165    /// # Arguments
166    /// * `s`: Shelf slope. A slope of `1.0` is maximally steep without overshoot.
167    pub fn shelf_slope(&mut self, s: T) -> &mut Self {
168        self.shape = Shape::Slope(s);
169        self
170    }
171
172    /// Set the shape
173    pub fn set_shape(&mut self, s: Shape<T>) -> &mut Self {
174        self.shape = s;
175        self
176    }
177
178    /// Get inverse Q
179    fn qi(&self) -> T {
180        match self.shape {
181            Shape::Q(qi) => qi.recip(),
182            Shape::Bandwidth(bw) => {
183                2.0.as_()
184                    * (T::LN_2() / 2.0.as_() * bw * self.frequency / self.frequency.sin()).sinh()
185            }
186            Shape::Slope(s) => {
187                ((self.gain + T::one() / self.gain) * (T::one() / s - T::one()) + 2.0.as_()).sqrt()
188            }
189        }
190    }
191
192    /// Get (cos(w0), alpha=sin(w0)/(2*q))
193    fn fcos_alpha(&self) -> (T, T) {
194        let (fsin, fcos) = self.frequency.sin_cos();
195        (fcos, 0.5.as_() * fsin * self.qi())
196    }
197
198    /// Low pass filter
199    ///
200    /// Builds second order biquad low pass filter coefficients.
201    ///
202    /// ```
203    /// use idsp::iir::*;
204    /// let ba = Filter::default()
205    ///     .critical_frequency(0.1)
206    ///     .gain(1000.0)
207    ///     .lowpass();
208    /// let iir = Biquad::<i32>::from(&ba);
209    /// let mut xy = [0; 4];
210    /// let x = vec![3, -4, 5, 7, -3, 2];
211    /// let y: Vec<_> = x.iter().map(|x0| iir.update(&mut xy, *x0)).collect();
212    /// assert_eq!(y, [5, 3, 9, 25, 42, 49]);
213    /// ```
214    pub fn lowpass(&self) -> [[T; 3]; 2] {
215        let (fcos, alpha) = self.fcos_alpha();
216        let b = self.gain * 0.5.as_() * (T::one() - fcos);
217        [
218            [b, (2.0).as_() * b, b],
219            [T::one() + alpha, (-2.0).as_() * fcos, T::one() - alpha],
220        ]
221    }
222
223    /// High pass filter
224    ///
225    /// Builds second order biquad high pass filter coefficients.
226    ///
227    /// ```
228    /// use idsp::iir::*;
229    /// let ba = Filter::default()
230    ///     .critical_frequency(0.1)
231    ///     .gain(1000.0)
232    ///     .highpass();
233    /// let iir = Biquad::<i32>::from(&ba);
234    /// let mut xy = [0; 4];
235    /// let x = vec![3, -4, 5, 7, -3, 2];
236    /// let y: Vec<_> = x.iter().map(|x0| iir.update(&mut xy, *x0)).collect();
237    /// assert_eq!(y, [5, -9, 11, 12, -1, 17]);
238    /// ```
239    pub fn highpass(&self) -> [[T; 3]; 2] {
240        let (fcos, alpha) = self.fcos_alpha();
241        let b = self.gain * 0.5.as_() * (T::one() + fcos);
242        [
243            [b, (-2.0).as_() * b, b],
244            [T::one() + alpha, (-2.0).as_() * fcos, T::one() - alpha],
245        ]
246    }
247
248    /// Band pass
249    ///
250    /// ```
251    /// use idsp::iir::*;
252    /// let ba = Filter::default()
253    ///     .frequency(1000.0, 48e3)
254    ///     .q(5.0)
255    ///     .gain_db(3.0)
256    ///     .bandpass();
257    /// println!("{ba:?}");
258    /// ```
259    pub fn bandpass(&self) -> [[T; 3]; 2] {
260        let (fcos, alpha) = self.fcos_alpha();
261        let b = self.gain * alpha;
262        [
263            [b, T::zero(), -b],
264            [T::one() + alpha, (-2.0).as_() * fcos, T::one() - alpha],
265        ]
266    }
267
268    /// A notch filter
269    ///
270    /// Has zero gain at the critical frequency.
271    pub fn notch(&self) -> [[T; 3]; 2] {
272        let (fcos, alpha) = self.fcos_alpha();
273        let f2 = (-2.0).as_() * fcos;
274        [
275            [self.gain, f2 * self.gain, self.gain],
276            [T::one() + alpha, f2, T::one() - alpha],
277        ]
278    }
279
280    /// An allpass filter
281    ///
282    /// Has constant `gain` at all frequency but a variable phase shift.
283    pub fn allpass(&self) -> [[T; 3]; 2] {
284        let (fcos, alpha) = self.fcos_alpha();
285        let f2 = (-2.0).as_() * fcos;
286        [
287            [
288                (T::one() - alpha) * self.gain,
289                f2 * self.gain,
290                (T::one() + alpha) * self.gain,
291            ],
292            [T::one() + alpha, f2, T::one() - alpha],
293        ]
294    }
295
296    /// A peaking/dip filter
297    ///
298    /// Has `gain*shelf_gain` at critical frequency and `gain` elsewhere.
299    pub fn peaking(&self) -> [[T; 3]; 2] {
300        let (fcos, alpha) = self.fcos_alpha();
301        let s = self.shelf.sqrt();
302        let f2 = (-2.0).as_() * fcos;
303        [
304            [
305                (T::one() + alpha * s) * self.gain,
306                f2 * self.gain,
307                (T::one() - alpha * s) * self.gain,
308            ],
309            [T::one() + alpha / s, f2, T::one() - alpha / s],
310        ]
311    }
312
313    /// Low shelf
314    ///
315    /// Approaches `gain*shelf_gain` below critical frequency and `gain` above.
316    ///
317    /// ```
318    /// use idsp::iir::*;
319    /// let ba = Filter::default()
320    ///     .frequency(1000.0, 48e3)
321    ///     .shelf_slope(2.0)
322    ///     .shelf_db(20.0)
323    ///     .lowshelf();
324    /// println!("{ba:?}");
325    /// ```
326    pub fn lowshelf(&self) -> [[T; 3]; 2] {
327        let (fcos, alpha) = self.fcos_alpha();
328        let s = self.shelf.sqrt();
329        let tsa = 2.0.as_() * s.sqrt() * alpha;
330        let sp1 = s + T::one();
331        let sm1 = s - T::one();
332        [
333            [
334                s * self.gain * (sp1 - sm1 * fcos + tsa),
335                2.0.as_() * s * self.gain * (sm1 - sp1 * fcos),
336                s * self.gain * (sp1 - sm1 * fcos - tsa),
337            ],
338            [
339                sp1 + sm1 * fcos + tsa,
340                (-2.0).as_() * (sm1 + sp1 * fcos),
341                sp1 + sm1 * fcos - tsa,
342            ],
343        ]
344    }
345
346    /// Low shelf
347    ///
348    /// Approaches `gain*shelf_gain` above critical frequency and `gain` below.
349    pub fn highshelf(&self) -> [[T; 3]; 2] {
350        let (fcos, alpha) = self.fcos_alpha();
351        let s = self.shelf.sqrt();
352        let tsa = 2.0.as_() * s.sqrt() * alpha;
353        let sp1 = s + T::one();
354        let sm1 = s - T::one();
355        [
356            [
357                s * self.gain * (sp1 + sm1 * fcos + tsa),
358                (-2.0).as_() * s * self.gain * (sm1 + sp1 * fcos),
359                s * self.gain * (sp1 + sm1 * fcos - tsa),
360            ],
361            [
362                sp1 - sm1 * fcos + tsa,
363                2.0.as_() * (sm1 - sp1 * fcos),
364                sp1 - sm1 * fcos - tsa,
365            ],
366        ]
367    }
368
369    /// I/HO
370    ///
371    /// Notch, integrating below, flat `shelf_gain` above
372    pub fn iho(&self) -> [[T; 3]; 2] {
373        let (fcos, alpha) = self.fcos_alpha();
374        let fsin = 0.5.as_() * self.frequency.sin();
375        let a = (T::one() + fcos) / (2.0.as_() * self.shelf);
376        [
377            [
378                self.gain * (T::one() + alpha),
379                (-2.0).as_() * self.gain * fcos,
380                self.gain * (T::one() - alpha),
381            ],
382            [a + fsin, (-2.0).as_() * a, a - fsin],
383        ]
384    }
385}
386
387// TODO
388// SOS cascades:
389// butterworth
390// elliptic
391// chebychev1/2
392// bessel
393
394#[cfg(test)]
395mod test {
396    use super::*;
397
398    use core::f64;
399    use num_complex::Complex64;
400
401    use crate::iir::*;
402
403    #[test]
404    #[ignore]
405    fn lowpass_noise_shaping() {
406        let ba = Biquad::<i32>::from(
407            &Filter::default()
408                .critical_frequency(1e-5f64)
409                .gain(1e3)
410                .lowpass(),
411        );
412        println!("{:?}", ba);
413        let mut xy = [0; 5];
414        for _ in 0..(1 << 24) {
415            ba.update(&mut xy, 1);
416        }
417        for _ in 0..10 {
418            ba.update(&mut xy, 1);
419            println!("{xy:?}");
420        }
421    }
422
423    fn polyval(p: &[f64], x: Complex64) -> Complex64 {
424        p.iter()
425            .fold(
426                (Complex64::default(), Complex64::new(1.0, 0.0)),
427                |(a, xi), pi| (a + xi * *pi, xi * x),
428            )
429            .0
430    }
431
432    fn freqz(b: &[f64], a: &[f64], f: f64) -> Complex64 {
433        let z = Complex64::new(0.0, -f64::consts::TAU * f).exp();
434        polyval(b, z) / polyval(a, z)
435    }
436
437    #[derive(Copy, Clone, Debug, PartialEq, PartialOrd)]
438    enum Tol {
439        GainDb(f64, f64),
440        GainBelowDb(f64),
441    }
442    impl Tol {
443        fn check(&self, h: Complex64) -> bool {
444            let g = 10.0 * h.norm_sqr().log10();
445            match self {
446                Self::GainDb(want, tol) => (g - want).abs() <= *tol,
447                Self::GainBelowDb(want) => g <= *want,
448            }
449        }
450    }
451
452    fn check_freqz(f: f64, g: Tol, ba: &[[f64; 3]; 2]) {
453        let h = freqz(&ba[0], &ba[1], f);
454        let hp = h.to_polar();
455        assert!(
456            g.check(h),
457            "freq {f}: response {h}={hp:?} does not meet {g:?}"
458        );
459    }
460
461    fn check_transfer(ba: &[[f64; 3]; 2], fg: &[(f64, Tol)]) {
462        println!("{ba:?}");
463
464        for (f, g) in fg {
465            check_freqz(*f, *g, ba);
466        }
467
468        // Quantize and back
469        let bai = (&Biquad::<i32>::from(ba)).into();
470        println!("{bai:?}");
471
472        for (f, g) in fg {
473            check_freqz(*f, *g, &bai);
474        }
475    }
476
477    #[test]
478    fn lowpass() {
479        check_transfer(
480            &Filter::default()
481                .critical_frequency(0.01)
482                .gain_db(20.0)
483                .lowpass(),
484            &[
485                (1e-3, Tol::GainDb(20.0, 0.01)),
486                (0.01, Tol::GainDb(17.0, 0.02)),
487                (4e-1, Tol::GainBelowDb(-40.0)),
488            ],
489        );
490    }
491
492    #[test]
493    fn highpass() {
494        check_transfer(
495            &Filter::default()
496                .critical_frequency(0.1)
497                .gain_db(-2.0)
498                .highpass(),
499            &[
500                (1e-3, Tol::GainBelowDb(-40.0)),
501                (0.1, Tol::GainDb(-5.0, 0.02)),
502                (4e-1, Tol::GainDb(-2.0, 0.01)),
503            ],
504        );
505    }
506
507    #[test]
508    fn bandpass() {
509        check_transfer(
510            &Filter::default()
511                .critical_frequency(0.02)
512                .bandwidth(2.0)
513                .gain_db(3.0)
514                .bandpass(),
515            &[
516                (1e-4, Tol::GainBelowDb(-35.0)),
517                (0.01, Tol::GainDb(0.0, 0.02)),
518                (0.02, Tol::GainDb(3.0, 0.01)),
519                (0.04, Tol::GainDb(0.0, 0.04)),
520                (4e-1, Tol::GainBelowDb(-25.0)),
521            ],
522        );
523    }
524
525    #[test]
526    fn allpass() {
527        check_transfer(
528            &Filter::default()
529                .critical_frequency(0.02)
530                .gain_db(-10.0)
531                .allpass(),
532            &[
533                (1e-4, Tol::GainDb(-10.0, 0.01)),
534                (0.01, Tol::GainDb(-10.0, 0.01)),
535                (0.02, Tol::GainDb(-10.0, 0.01)),
536                (0.04, Tol::GainDb(-10.0, 0.01)),
537                (4e-1, Tol::GainDb(-10.0, 0.01)),
538            ],
539        );
540    }
541
542    #[test]
543    fn notch() {
544        check_transfer(
545            &Filter::default()
546                .critical_frequency(0.02)
547                .bandwidth(2.0)
548                .notch(),
549            &[
550                (1e-4, Tol::GainDb(0.0, 0.01)),
551                (0.01, Tol::GainDb(-3.0, 0.02)),
552                (0.02, Tol::GainBelowDb(-140.0)),
553                (0.04, Tol::GainDb(-3.0, 0.02)),
554                (4e-1, Tol::GainDb(0.0, 0.01)),
555            ],
556        );
557    }
558
559    #[test]
560    fn peaking() {
561        check_transfer(
562            &Filter::default()
563                .critical_frequency(0.02)
564                .bandwidth(2.0)
565                .gain_db(-10.0)
566                .shelf_db(20.0)
567                .peaking(),
568            &[
569                (1e-4, Tol::GainDb(-10.0, 0.01)),
570                (0.01, Tol::GainDb(0.0, 0.04)),
571                (0.02, Tol::GainDb(10.0, 0.01)),
572                (0.04, Tol::GainDb(0.0, 0.04)),
573                (4e-1, Tol::GainDb(-10.0, 0.05)),
574            ],
575        );
576    }
577
578    #[test]
579    fn highshelf() {
580        check_transfer(
581            &Filter::default()
582                .critical_frequency(0.02)
583                .gain_db(-10.0)
584                .shelf_db(-20.0)
585                .highshelf(),
586            &[
587                (1e-6, Tol::GainDb(-10.0, 0.01)),
588                (1e-4, Tol::GainDb(-10.0, 0.01)),
589                (0.02, Tol::GainDb(-20.0, 0.01)),
590                (4e-1, Tol::GainDb(-30.0, 0.01)),
591            ],
592        );
593    }
594
595    #[test]
596    fn lowshelf() {
597        check_transfer(
598            &Filter::default()
599                .critical_frequency(0.02)
600                .gain_db(-10.0)
601                .shelf_db(-20.0)
602                .lowshelf(),
603            &[
604                (1e-6, Tol::GainDb(-30.0, 0.01)),
605                (1e-4, Tol::GainDb(-30.0, 0.01)),
606                (0.02, Tol::GainDb(-20.0, 0.01)),
607                (4e-1, Tol::GainDb(-10.0, 0.01)),
608            ],
609        );
610    }
611
612    #[test]
613    fn iho() {
614        check_transfer(
615            &Filter::default()
616                .critical_frequency(0.01)
617                .gain_db(-20.0)
618                .shelf_db(10.0)
619                .q(10.)
620                .iho(),
621            &[
622                (1e-5, Tol::GainDb(40.0, 0.01)),
623                (0.01, Tol::GainBelowDb(-40.0)),
624                (4.99e-1, Tol::GainDb(-10.0, 0.01)),
625            ],
626        );
627    }
628}
idsp/iir/coefficients.rs

idsp/iir/
coefficients.rs
