Linear filter

[barakugav]

Hey, the library looks amazing, thank you for providing these capabilities!

I want to apply a high pass filter, and the design module looks perfect for it, but how do i apply the filter after designing it? Im used to Python scipy.signal.lfilter, does the library expose such function? sosfiltfilt_dyn looks like it does the job for second order filters only.

Thanks in advance

[trueb2]

Thanks!

We don't have lfilter_dyn implemented right now because cascaded biquads (second order sections) are more numerically stable to apply.

From SciPy (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html)

The function sosfilt (and filter design using output='sos') should be preferred over lfilter for most filtering tasks, as second-order sections have fewer numerical problems

You could use sosfilt_dyn if you don't care about the phase delay. If you do care for zero phase shift, then sosfiltfilt_dyn is the one to use.

If the filter was already designed into ba, then you would need lfilter_dyn implemented

[barakugav]

Thanks a lot, will use the SOS filters. I also required an option to pass the argument zi (argument of both lfilter and sosfilter) in scipy. Does the library support passing it?

[trueb2]

The second order sections (Sos instances) are the coefficients and internal state right now: https://docs.rs/sci-rs/latest/src/sci_rs/signal/filter/design/sos.rs.html#22-33

That means that you can initialize Sos instances (with zi set to 0). The instances hold the state from the last invocation. I haven't hit a use case where you want to start filtering from a known zi that isn't already based on previous passes, so there isn't currently an interface to construct with b, a, and zi.

https://docs.rs/sci-rs/latest/src/sci_rs/signal/filter/design/sos.rs.html#52-58 only takes the coefficients at the moment.

[barakugav]

I didnt realize the state is stored within the SOS instances. That fits my need exactly, thanks! As you said, i only need to initialize the state with zeros.

I used butter_dyn successfully with FilterBandType::Bandstop, but i get an error in another place in my code when designing a high pass filter using iirfilter_dyn:

let sample_rate = 600;
let highpass_filter_order = 5;
let highpass_freq = 2;
let nyq = sample_rate as f64 / 2.0;
let filter = iirfilter_dyn(
   highpass_filter_order as usize,
    vec![highpass_freq as f64 / nyq],
    None,
    None,
    Some(FilterBandType::Highpass),
    None,
    None,
    Some(FilterOutputType::Sos),
    None,
);

The panic i get is:

thread ... panicked at .../.cargo/registry/src/index.crates.io-6f17d22bba15001f/sci-rs-0.3.12/src/signal/filter/design/zpk2sos.rs:246:5:
assertion failed: p.is_empty()

Am i passing the wrong parameters?

Thanks again for helping, the library really solve my need.

[trueb2]

I don’t think the sample rates make sense there. Normally I would provide separate arguments for the cutoff (2Hz) and the sample rate (600Hz). No reason to do the division before the function.

I don’t normally use the iirfilter_dyn interface directly since it is called by butter_dyn.

The iirfilter function should mirror the behavior of the scipy iirfilter function

[barakugav]

Thanks for the quip response. I tried both passing the sample frequency to the iirfilter_dyn:

let filter = iirfilter_dyn(
    5,
    vec![2.0],
    None,
    None,
    Some(FilterBandType::Highpass),
    None,
    None,
    Some(FilterOutputType::Sos),
    Some(600.0),
);

And also using the butter_dyn function:

let filter = butter_dyn(
    5,
    vec![2.0],
    Some(FilterBandType::Highpass),
    Some(false),
    Some(FilterOutputType::Sos),
    Some(600.0),
);

Both panic exactly the same way. The filter is successfully created with the Zpk format but fails the conversion to sos.

I already use these exact arguments using the same functions call in Python and it works, so maybe there is a bug or i am passing the arguments in Rust wrong. I can also create the SOS coefficients in Python and hardcode it in Rust, but i prefer not to.

I appreciate the help!

[barakugav]

I found the bug, in zpk2sos_dyn, the line

if p1.im.is_zero() && p.iter().map(|pi| pi.re).sum::<F>().is_zero() {

Should be

if p1.im.is_zero() && p.iter().filter(|pi| pi.im.is_zero()).count() == 0 {

Matching the scipy code:

if np.isreal(p1) and np.isreal(p).sum() == 0:

In addition, i had to multiply the k of Zpk format by -1, didnt investigate why, but with both of these changes i get the same exact filtering as the scipy one.

I think the following python line was not translated correctly in lp2hp_zpk_dyn:

k_hp = k * real(prod(-z) / prod(-p))

[barakugav]

https://github.com/qsib-cbie/sci-rs/pull/41

[trueb2]

Here is the relevant code on the k value

    // # Cancel out gain change caused by inversion
    // k_hp = k * real(prod(-z) / prod(-p))

    let num = zpk
        .z
        .iter()
        .map(|zi| Complex::new(one_neg, F::zero()) + *zi)
        .fold(Complex::new(F::one(), F::zero()), |acc, zi| acc * zi);
    let denom = zpk
        .p
        .iter()
        .map(|pi| Complex::new(F::one(), F::zero()) * *pi)
        .fold(Complex::new(F::one(), F::zero()), |acc, pi| acc * pi);
    let k_hp = zpk.k * (num / denom).real();

Looks like the .map(|pi| Complex::new(F::one(), F::zero()) * *pi) is incorrect and should one_neg like above

[trueb2]

Indeed the p negation is an issue. I will add test cases and drive this and your PR to completion.

[trueb2]

Fixes are merged and tested against SciPy output. I released a new crate version (0.3.13)

[barakugav]

Amazing! @trueb2 Thanks for all the help!