idsp/iir/
coefficients.rs

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

idsp/iir/
coefficients.rs
