[SciPy-Dev] Improving scaling factors in IIR SOS sections in scipy.signal

[Eric Larson]

>
> Would such a change be acceptable upstream? And what would be the best way
> to handle "backward compatibility"?
>

I was originally thinking that a new kwarg for the SOS-generation functions
would be okay. But this would involve modifying potentially many functions,
and at the same time seems like a separable step. In other words, once you
get the `sos` array from the generation function (or even construct one
yourself), you can directly figure out the gain factors of the sections
(whatever they happen to be), compute the total gain, and then redistribute
it among the sections however you want.

Maybe a new `sos_even = distribute_gain(sos)` would be good. Then if there
are other gain-distribution methods that are relevant eventually we could
add a `method='even' | 'first'` kwarg. Perhaps we'd even want this `method`
kwarg from the start so you could test that round-trip calls with `'even'`
and then `'first'` give back the original `sos` you got from the
construction function.

I expect that this would only be a few lines (loop plus NumPy arithmetic),
but it sounds generally useful enough to include, especially if MATLAB (and
Octave) do this redistribution by default.

(1) I was playing with signal.sosfiltfilt(...) to apply an IIR filter
> to the input signal. I didn't check the source code, but the results
> look as if the filter was working partially "backwards".
>

An IIR filter can use the same parts that a FIR filter can (i.e., as many
values of the input as it wants) plus its own previous output values. If
you want the IIR filter to do only one or the other, you restrict the
numerator or denominator of the transfer function. Since you say you're a
newcomer to digital filtering, I'd recommend reading up a bit on
z-transform transfer function basics to get a handle on this if you're
interested.

(2) I would like to simulate the filter with integer (fixed point)
> arithmetics.
>
...

(4) I have an ad-hoc simulation for filters with integer arithmetics.
> But it might be nice to add this support to scipy.signal as well, as
> this functionality might be generally useful. I can work on this, but
> I wouldn't mind a little bit of guidance on the subject if I start
> working in that direction. Anyway, I would start by fixing the
> coefficients first, and leave this one for the end.


After a little bit of searching the most relevant thing I could find is
"demodel" for modeling fixed-point in Python, but I'm not sure how up to
date it is. This sort of functionality in general seems like it would be
useful, but perhaps better suited to its own package, at least to start, as
I suspect it will require a lot of (fixed-point) domain-specific knowledge
and testing to get right. For it to be included in SciPy we'd want to have
someone with this domain expertise (or willingness to acquire enough of it)
to commit to helping maintain the code -- I don't think any current
maintainers have it (please correct me if I'm wrong!), and we don't want to
be stuck in a position of maintaining code we don't understand :)

(3) Matlab seems to have some support for converting coefficients from
> sos sections into integers, and also tells you how many bits are
> needed to perform fixed-point arithmetics. It offers one specific way
> of rounding the coefficients to keep the calculations stable
>

I don't know anything about this, but in general a workable path has been:

1. Look at textbooks, MATLAB, or Octave for references (cannot look at
their code due to license restrictions) for a given algorithm
2. Look for a BSD-compatible implementation that exists in Python (maybe
not polished but could be adapted to SciPy)
3. If none exists, implement the algorithm in Python/SciPy from scratch
4. Test against values produced in MATLAB/Octave.

Eric
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20200331/b926adc4/attachment.html>