Improving scaling factors in IIR SOS sections in scipy.signal
Dear signal processing gurus, I'm using python for playing around, testing and quickly prototyping IIR filters for a piece of code which is otherwise written in C[++] & FPGA, but I miss one thing in zpk2sos, namely the complete scaling factor ends up in the first sos section, rather than being "evenly distributed" across all the sections. For floating point operations this doesn't make much of a difference, but not applying the scaling inside a sos section with integer arithmetics is potentially more problematic to work with (and I would like to have support for this to be able to compare calculations & results). Example:
signal.butter(8, 3E-4, output='sos') array([[ 2.42594124e-27, 4.85188247e-27, 2.42594124e-27, 1.00000000e+00, -1.99815208e+00, 9.98152971e-01], [ 1.00000000e+00, 2.00000000e+00, 1.00000000e+00, 1.00000000e+00, -1.99843306e+00, 9.98433944e-01], [ 1.00000000e+00, 2.00000000e+00, 1.00000000e+00, 1.00000000e+00, -1.99895244e+00, 9.98953323e-01], [ 1.00000000e+00, 2.00000000e+00, 1.00000000e+00, 1.00000000e+00, -1.99963144e+00, 9.99632331e-01]])
I would like to see the 2.4259e-27 part to be distributed among the four sos sections (rather than ending up with [1 2 1] everywhere else), which is also what Matlab returns. Since I heavily prefer using Python compared to Matlab (I don't even have the licence myself), I wanted to ask for your opinion before I start writing the code. Would such a change be acceptable upstream? And what would be the best way to handle "backward compatibility"? (I've never contributed to Python before, but I assume it's doable :) A bunch of unrelated questions: (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". That is: it seems to be using some points from the "future" to determine the result, similar to how a FIR filter would take some points from the past + some points from the future to calculate the current point. I can send some screenshot, but I'm obviously missing something important here. Is there another approach to exactly simulate out[i] = b[0] * in[i] + b[1] * in[i-1] + b[2] * in[i-2] - a[1] * out[i-1] - a[2] * out[i-2] for the output? (2) I would like to simulate the filter with integer (fixed point) arithmetics. For a low-pass filter with either very low or very high cut-off frequency it's way too easy to either get into integer overflows, or be forced into cutting a lot of significant bits off. What's the least painful way to perform multiplications with a slightly higher bit count? Sadly something like dtype=np.int128 doesn't seem to exist :( I'm currently having some troubles with LP filters with cut-off frequency very close to the Nyquist frequency, and I'm thinking of perhaps improving the signal package first before investigating the problem any further, as that would help me simulate and compare where things go wrong. (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'm sorry, I don't remember what exactly it says, a friend showed it to me a while ago, but we no longer have access to that computer with Matlab installed). Does anyone know anything about this? I suspect it tries to round complex numbers close to the unit circle in a way that ensures the points remain inside the circle, but that's pure speculation. (It would be nice to be able to offer something similar on the python side.) (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. Thank you very much, Mojca PS: I'm relatively new to the subject of filtering digital signals, so please bear with me :)
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
Hi, On Tue, 31 Mar 2020 at 16:32, Eric Larson wrote:
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.
In the meantime I managed to get access to a machine with Matlab installed and also tested what Octave does. I get the "correct" / full coefficients with fdatool (a graphical tool; see below for results), while > pkg load signal > [z,p,k] = butter(4, 3e-4); > sos = zp2sos(z,p,k); > sos sos = 4.9253e-14 9.8505e-14 4.9253e-14 1.0000e+00 -1.9983e+00 9.9826e-01 1.0000e+00 2.0000e+00 1.0000e+00 1.0000e+00 -1.9993e+00 9.9928e-01 sadly gives me the results equivalent to python :( With fdatool I can pick between "double precision", "single precision" and "fixed point" (with some additional settings), but let's say that ending up with double-precision is good enough for now. This is output from fdatool after some very basic point-and-click-ing: % Discrete-Time IIR Filter (real) % ------------------------------- % Filter Structure : Direct-Form II, Second-Order Sections % Number of Sections : 2 % Stable : Yes % Linear Phase : No SOS Matrix: 1 2 1 1 -1.9927241098599966 0.99281261642961327 1 2 1 1 -1.9826478027009247 0.98273586173273308 Scale Values: 0.000022126642404234759 0.000022014757952111739 You can see that 4.9e-14 is split into two sections (2.2e-7), with each section serving as its own standalone sos section with correct gain.
(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.
Sure, but generally FIR filters are (very) long and the result would be something like y[t] = a[- n] * x[t - n] + a[-n + 1] * x[t - n + 1] + ... + a[n] * x[t + n] with y[t] depending on values from the future. While one can have an IIR filter of even the first or the second order with y[t] = b[0] * x[t] + b[1] * x[t - 1] + b[2] * x[t - 2] - a[1] * y[t - 1] - a[2] * y[t - 2] which should depend exclusively on the values from the past. What I didn't particularly like was the following behaviour: sos = signal.butter(4, 3e-4, output='sos') plt.plot(signal.sosfiltfilt(sos, np.concatenate([np.zeros(2000), np.ones(100000)]))) plt.plot(signal.sosfiltfilt(sos, np.concatenate([np.zeros(10000), np.ones(100000)]))) plt.show() I wanted both to produce the same shape. But I believe I now found the solution at https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.sosfiltfil... with zi = x[:4].mean() * sosfilt_zi(sos) y2, zo = sosfilt(sos, x, zi=zi) When comparing the results with my own filter implementation and the results were wildly different, I forgot to search for the proper solution to correctly initialize the filter. (The result without initialisation looked more like a FIR filter.)
(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 :)
After staring at the Matlab interface for a while (today I had the first chance to access a machine) and seeing the gazillion of different option I tend to agree with this. There's probably a very large gap between my own ad-hoc implementation to get the job done (which is already working using just a few lines of code), and the feature-complete variant.
(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.
I now found for example Chapter 7 (Quantized Filter Analysis) in B. A. Shenoi: Introduction to Digital Signal Processing and Filter Design, and I guess there are a lot more books on this topic. I'll look through it, but it's certainly tons more work compared to just "fixing" the gain factors, so maybe too much for now. But just for completeness and to answer my own question. Here are the rounding options offered by Matlab: - Ceiling - Nearest - Nearest (convergent) - Round - Zero - Floor (there's a gazillion of other options to fine-tune the results). That above books says: "The convergent operation is the same as rounding except that in the case when the number is exactly halfway, it is rounded down if the penultimate bit is zero and rounded up if it is one." I'll come up with something for gain correction first and ask for feedback once I have some code to show. Mojca
participants (2)
-
Eric Larson -
Mojca Miklavec