idsp 0.22.0

DSP algorithms for embedded, mostly integer math
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

//! Frequency-response utilities for second-order sections.

use crate::Complex;
#[cfg(not(feature = "std"))]
#[allow(unused_imports)]
use num_traits::float::Float as _;
use num_traits::{AsPrimitive, Float, FloatConst};

use crate::iir::Biquad;

fn polyval<T: Float>(p: &[T], z: Complex<T>) -> Complex<T> {
    p.iter()
        .rev()
        .fold(Complex::new(T::zero(), T::zero()), |a, &p| {
            a * z + Complex::new(p, T::zero())
        })
}

/// Evaluate a `[b, a]` transfer function on the unit circle.
///
/// Evaluates
/// `H(z) = (sum_k b[k] z^-k) / (sum_k a[k] z^-k)` at `z = exp(-j 2 pi frequency)`.
///
/// `frequency` is relative to the sample rate.
///
/// `a` and `b` use the literature/cookbook sign convention, not the internal
/// [`Biquad`] storage convention. `a` must be non-empty.
pub fn freqz<T: Float + FloatConst>(b: &[T], a: &[T], frequency: T) -> Complex<T> {
    assert!(!a.is_empty(), "freqz requires a non-empty denominator");
    let (i, r) = (-T::TAU() * frequency).sin_cos();
    let z = Complex::new(r, i);
    let q = polyval(a, z);
    polyval(b, z) * q.conj() / (q.re() * q.re() + q.im() * q.im())
}

impl<C> Biquad<C> {
    /// Evaluate the transfer function on the unit circle.
    ///
    /// `frequency` is relative to the sample rate.
    pub fn freqz<T: 'static + Float + FloatConst>(&self, frequency: T) -> Complex<T>
    where
        C: Copy + AsPrimitive<T>,
    {
        let [b0, b1, b2, a1, a2] = self.ba.map(AsPrimitive::as_);
        freqz(&[b0, b1, b2], &[T::one(), -a1, -a2], frequency)
    }
}