0
$\begingroup$

I'm having a lot of fun writing signal processing code in Python / numpy, and I'm resisting the urge to pre-optimize the code. But my biquad implementation is slower than I want. Here's the inner loop -- it uses afor loop: is there a way to vectorize it? Alternatively, could I use a scipy filter, and if so, how do I translate a0, a1, a2, b0, b1, b2?

        # src_frames is an ndarray containing source samples
        # dst_frames is an equal size ndarray to receive the filtered samples
        # self._x[] and self._y[] are the feedforward and feedback delay elements
        for i, x in enumerate(src_frames):
            # compute one output sample
            y = (b0 * x + b1 * self._x[1] + b2 * self._x[2] -
                          a1 * self._y[1] - a2 * self._y[2])
            # shift delay elements
            self._x[2] = self._x[1]
            self._x[1] = x
            self._y[2] = self._y[1]
            self._y[1] = y
            dst_frames[i] = y

(Yes, I know I'm not using _x[0] and _y[0], but it seemed clearer to keep the indices numbered with the coefficients.)

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ If applying to 1D data, you can use scipy.signal.lfilter which allows you to specify the numerator and denominator coefficient vectors b and a, respectively. Keep in mind DSP.SE is more geared toward questions about DSP rather than getting coding tips. Stack overflow would be better suited for these questions in general. $\endgroup$
    – Ash
    Dec 11, 2022 at 3:37
  • $\begingroup$ Good comment -- thank you for both suggestions. $\endgroup$
    – fearless_fool
    Dec 11, 2022 at 12:04

1 Answer 1

Reset to default
2
$\begingroup$

IIR filters don't vectorize easily because they are recursive in nature.

Given that IIR filters are so important almost all high level language distributions come with a library that has a heavily optimized filter function. Most of this optimization happens close to the hardware level: utilizing the specific instruction set of the processor, pipeline optimization, memory access optimization, register usage, caching of state variables, loop unrolling, etc.

Options for high level optimizations are limited. Some examples are

  • Split into cascaded second order sections
  • On most processors Transposed Form II is a bit faster than Direct Form I since it has less state variables. (Do NOT use the other two)
  • Create multi channel filters: For example for a stereo file can process both channels at the same time.
  • Choose "processor friendly" frame sizes. Typically multiples of 4 or 8 work a little better

In your case it's probably best just go with https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.sosfilt.html

Improve this answer
$\endgroup$
2
  • $\begingroup$ I've been delving into scipy.signal.sosfilt, and I agree that's the most sensible way to go. I now understand how the ax and bx parameters get passed to sosfilt -- it's a direct mapping. $\endgroup$
    – fearless_fool
    Dec 11, 2022 at 12:03
  • $\begingroup$ Note that splitting the filter into cascaded second-order sections is something you want to do anyway, for numerical reasons. It would be interesting to benchmark a $n^{th}$-order filter using 2nd-order sections against an $n^{th}$ order state-space system that realizes the same transfer function. In something like Scipy it may actually be faster (it was on the TMS320F2812, up to somewhere between 4th and 8th-order). $\endgroup$
    – TimWescott
    Dec 12, 2022 at 1:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.