BUG: Signal.butter in analog mode has all the gain in the first sos

[Pi9GW]

Describe your issue.

Filter coefficients in analog mode should be normalized to section gains of one (1). Practically they are not realizable and plotting the section frequency responses without normalization is also not too nice. Plot figures attached, Y in dB. Top default, bottom normalized.

Reproducing Code Example

>>> signal.butter(2, 200, analog=True, output='sos')
array([[0.00000000e+00, 0.00000000e+00, 4.00000000e+04, 1.00000000e+00,
        2.82842712e+02, 4.00000000e+04]])
>>> signal.butter(4, 200, analog=True, output='sos')
array([[0.00000000e+00, 0.00000000e+00, 1.60000000e+09, 1.00000000e+00,
        3.69551813e+02, 4.00000000e+04],
       [0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 1.00000000e+00,
        1.53073373e+02, 4.00000000e+04]])
>>> signal.butter(6, 200, analog=True, output='sos')
array([[0.00000000e+00, 0.00000000e+00, 6.40000000e+13, 1.00000000e+00,
        3.86370331e+02, 4.00000000e+04],
       [0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 1.00000000e+00,
        2.82842712e+02, 4.00000000e+04],
       [0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 1.00000000e+00,
        1.03527618e+02, 4.00000000e+04]])

Error message

None

SciPy/NumPy/Python version information

1.8.0 1.21.5 sys.version_info(major=3, minor=10, micro=6, releaselevel='final', serial=0)

[Pi9GW]

There is more in this bug than already told. Gains in lowpass zpk format are off completely. Highpass works and the original bug is not there. Something else is broken as freqs picks about 1000 times too low frequency range. Plot shows default behavior. All this is present in SciPy 1.10.0 as well.

>>> butter(2, 1050, btype='low', analog=True, output='zpk')
(array([], dtype=float64), array([-742.46212025+742.46212025j, -742.46212025-742.46212025j]), 1102500.0)
>>> butter(2, 1050, btype='high', analog=True, output='zpk')
(array([0., 0.]), array([-742.46212025-742.46212025j, -742.46212025+742.46212025j]), 1.0)

All analog gains should be one including Chebyshev and Elliptic...

>>> butter(2, 200, btype='high', analog=True, output='sos')
array([[1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
        2.82842712e+02, 4.00000000e+04]])
>>> butter(4, 200, btype='high', analog=True, output='sos')
array([[1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
        3.69551813e+02, 4.00000000e+04],
       [1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
        1.53073373e+02, 4.00000000e+04]])
>>> butter(6, 200, btype='high', analog=True, output='sos')
array([[1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
        3.86370331e+02, 4.00000000e+04],
       [1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
        2.82842712e+02, 4.00000000e+04],
       [1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00,
        1.03527618e+02, 4.00000000e+04]])

[roryyorke]

There is more in this bug than already told. Gains in lowpass zpk format are off completely. Highpass works and the original bug is not there. Something else is broken as freqs picks about 1000 times too low frequency range.

I don't agree: that non-unity gain is how the zpk representation works, i.e.:

$$ g(s) = k\frac{\prod_i s-z_i}{\prod_j s-p_j} $$

so, to get unity gain at DC (where s=0), we need $k = \prod_j -p_j$, as it is in your example:

In [1]: import scipy

In [2]: from scipy import signal

In [3]: scipy.__version__
Out[3]: '1.10.0'

In [4]: signal.tf2zpk([10],[1,10])
Out[4]: (array([], dtype=float64), array([-10.]), 10.0)

In [6]: zpk = signal.butter(2, 1050, btype='low', analog=True, output='zpk')

In [8]: signal.zpk2tf(*zpk)
Out[8]: (array([1102500.]), array([1.00000000e+00, 1.48492424e+03, 1.10250000e+06]

One might wonder why we can't have a convention of unity gain at DC, but then we wouldn't be able to have poles or zeros at the origin.

[roryyorke]

"All the gain in the first second-order section" is true for all combinations of (analog, digital) and (butter, cheby, cheby2, ellip) (I'm not sure about bessel, but likely: think they all use zpk2sos).

Consider, e.g., the numerator of magnitude $10^{-14}$ in this elliptical filter:

In [18]: signal.iirdesign(0.09, 0.11, 0.1, 300, ftype='ellip', output='sos')
Out[18]: 
array([[ 1.64825222e-14,  6.03310059e-15,  1.64825222e-14,
         1.00000000e+00, -9.67090642e-01,  0.00000000e+00],
       [ 1.00000000e+00,  1.00000000e+00,  0.00000000e+00,
         1.00000000e+00, -1.93371066e+00,  9.36370425e-01],
       [ 1.00000000e+00, -9.10214404e-01,  1.00000000e+00,
         1.00000000e+00, -1.93236515e+00,  9.39544725e-01],
       [ 1.00000000e+00, -1.40194610e+00,  1.00000000e+00,
         1.00000000e+00, -1.93032634e+00,  9.44394208e-01],
       [ 1.00000000e+00, -1.61526902e+00,  1.00000000e+00,
         1.00000000e+00, -1.92784655e+00,  9.50377038e-01],
       [ 1.00000000e+00, -1.72276085e+00,  1.00000000e+00,
         1.00000000e+00, -1.92519540e+00,  9.56922819e-01],
       [ 1.00000000e+00, -1.78326738e+00,  1.00000000e+00,
         1.00000000e+00, -1.92261599e+00,  9.63532375e-01],
       [ 1.00000000e+00, -1.82008832e+00,  1.00000000e+00,
         1.00000000e+00, -1.92030038e+00,  9.69834981e-01],
       [ 1.00000000e+00, -1.84374256e+00,  1.00000000e+00,
         1.00000000e+00, -1.91838485e+00,  9.75602650e-01],
       [ 1.00000000e+00, -1.85949227e+00,  1.00000000e+00,
         1.00000000e+00, -1.91695883e+00,  9.80734540e-01],
       [ 1.00000000e+00, -1.87018494e+00,  1.00000000e+00,
         1.00000000e+00, -1.91608075e+00,  9.85227466e-01],
       [ 1.00000000e+00, -1.87745292e+00,  1.00000000e+00,
         1.00000000e+00, -1.91579498e+00,  9.89144839e-01],
       [ 1.00000000e+00, -1.88226630e+00,  1.00000000e+00,
         1.00000000e+00, -1.91614707e+00,  9.92590829e-01],
       [ 1.00000000e+00, -1.88520530e+00,  1.00000000e+00,
         1.00000000e+00, -1.91719603e+00,  9.95692024e-01],
       [ 1.00000000e+00, -1.88659980e+00,  1.00000000e+00,
         1.00000000e+00, -1.91902361e+00,  9.98586109e-01]])

It would be interesting to know if this is a problem in application of the filter; my first, unexamined guess is no, at least for digital filters, since it applies to all elements of the numerator, so the "down-scaling" is uniform for the first second-order section.

One case the scaling could be a real problem is in cascaded second-order electronic filters: in that case one would probably want to keep the pass-band gain close to 1 to avoid hitting power-supply limits on the one hand, and descending below the noise floor on the other hand.

I think the numerators could be exactly rescaled using numpy.frexp and numpy.ldexp to distribute the exponent between the sections ("exactly" here is from a floating-point point-of-view.)

Here's a proof-of-concept:

import numpy as np


def rescalesos(sos):
    # Strictly proof-of-concept
    #  - doesn't handle general analog (analog 1st order section has leading 0 in numerator)
    #    - maybe check for largest absolute exponent in whole numerator?
    #  - assumes SOS is of form
    #     [k, .... ],
    #     [1, .... ],
    #     [1, .... ],
    #     [..., .... ],
    #     [1, .... ],
    #     - fix this by summing over the numerators
    sos = sos.copy()
    
    mnum, enum = np.frexp(sos[:,:3])
    q, r = divmod(enum[0,0], sos.shape[0])

    enum[0] -= enum[0,0]
    enum[:r] += q+1
    enum[r:] += q

    sos[:,:3] = np.ldexp(mnum, enum)

    return sos


if __name__ == '__main__':
    from scipy import signal
    sos = signal.iirdesign(0.09, 0.11, 0.1, 300, ftype='ellip', output='sos')
    resos = rescalesos(sos)

    z, p, k = signal.sos2zpk(sos)
    z2, p2, k2 = signal.sos2zpk(resos)
    
    assert np.all(z == z2)
    assert np.all(p == p2)
    assert np.all(k == k2)    

    with np.printoptions(precision=3):
        print(resos)

which on my system gives

[[ 0.072  0.027  0.072  1.    -0.967  0.   ]
 [ 0.125  0.125  0.     1.    -1.934  0.936]
 [ 0.125 -0.114  0.125  1.    -1.932  0.94 ]
 [ 0.125 -0.175  0.125  1.    -1.93   0.944]
 [ 0.125 -0.202  0.125  1.    -1.928  0.95 ]
 [ 0.125 -0.215  0.125  1.    -1.925  0.957]
 [ 0.125 -0.223  0.125  1.    -1.923  0.964]
 [ 0.125 -0.228  0.125  1.    -1.92   0.97 ]
 [ 0.125 -0.23   0.125  1.    -1.918  0.976]
 [ 0.125 -0.232  0.125  1.    -1.917  0.981]
 [ 0.125 -0.234  0.125  1.    -1.916  0.985]
 [ 0.125 -0.235  0.125  1.    -1.916  0.989]
 [ 0.125 -0.235  0.125  1.    -1.916  0.993]
 [ 0.125 -0.236  0.125  1.    -1.917  0.996]
 [ 0.125 -0.236  0.125  1.    -1.919  0.999]]

[roryyorke]

FWIW, to whoever can triage this, I see it as more of an enhancement than a defect.

[Pi9GW]

To be fair I am not an expert of any kind in this, just trying to get some real world designs done for amateur radio. I have obviously stated some things without complete understanding of the subject.

Regarding the bug, why the behavior is different for highpass (that seems proper to me) and lowpass?

I am using Scipy for cascaded electronic filter design, so it would make sense to have default section (DC) gains of one.

To save some amplifiers I also combine second order sections to fourth order ones, with combined sections being optimized for dynamic range. E.g. for 8th order filter with two 4th order sections, one with sos 2, 3 and the other with 1, 4. Actual component values calculated by Differential Evolution optimizer for the transfer function (that is from the two combined sos). So far the results have been amazing and totally not possible without Scipy, for me at least. Thanks for everyone involved!

[roryyorke]

Regarding the bug, why the behavior is different for highpass (that seems proper to me) and lowpass?

The IIR filter designs all start by computing zeros, poles, and a gain k. This is converted to SOS if requested by the user, using sos2zpk, and sos2zpk always puts the gain in the first second-order section. For a high-pass filter $g(s)$, we expect $\lim_{\omega\rightarrow\infty} g(jw) \approx 1$, and given the definition of SOS I gave before, for a system where len(z)=len(p), which we expect for a HPF, $\lim_{\omega\rightarrow\infty} g(jw) = k$, so we expect k to be about 1. We can check that:

In [343]: for ftype in ['butter', 'cheby1', 'cheby2', 'ellip']: print(signal.iirdesign(10, 9, 0.1, 30, ftype=ftype, output='zpk', analog=True)[2])
1.0000000000000004
1.0
1.0
0.988553094656939

compare to the complementary low-pass filters, which have wildly varying gains between the types:

In [344]: for ftype in ['butter', 'cheby1', 'cheby2', 'ellip']: print(signal.iirdesign(9, 10, 0.1, 30, ftype=ftype, output='zpk', analog=True)[2])
3.039172424287326e+49
4066118519.7896285
4.106622065499582
0.03162277660168413

I am using Scipy for cascaded electronic filter design, so it would make sense to have default section (DC) gains of one.

This is a good idea, but it doesn't work in general; in particular, not all LPFs have a DC gain of 1 (ellip and cheby1 have gain <1 for even orders). And it definitely won't work for filters that aren't LPFs. However, you can do this normalization yourself after calling iirdesign, etc.

We could special-case butter and cheby2 LPFs to do this normalization, but I think the complexity cost for testing and maintenance would not justify it. I think there is a case to be made to distribute the gains between the second-order sections automatically, especially for analog filters. Support for analog SOS is relatively recent, so it's not too surprising that we're missing a few features. (I see, e.g., that there's no sos_freqs.)

FWIW, the method I posted 3 weeks ago is almost certainly overkill; one could just geometrically divided k between the sections using an nth-root, something like this (untested):

    k = sos[0,2]
    sos[0,:2] /= k
    sos[:,:2] *= k**(1/sos.shape[0])

[Pi9GW]

Regarding DC gains I meant 0 Hz for LP, ∞Hz for HP, center gain for BP, etc.

Implementation in analog electronics makes use of e.g. amplifiers with (voltage) gain of one, an example being Sallen–Key architecture [1]. With such, the LP 0Hz gain is one no matter the function, so any ripples will be on the above one gain side. Depending on the application the section gains may be scaled, but quite often a compromise can be found by rearranging the order, especially if the interfering signal properties are somewhat known. Quite often the DC gain is most important e.g. for measurements near 0 Hz (I consider DC a theoretical signal).

For digital implementation many issues of the analog world disappear or at least change markedly. I guess that is why the gain in first stage is maximized. I still do not get how it would not interfere with filtering as it will limit the input dynamic range greatly. Is it expected all input is below one? Overflows will result otherwise. In the end just the same problems as in analog side, just with stable coefficients ;^)

[1] https://en.wikipedia.org/wiki/Sallen%E2%80%93Key_topology

[Pi9GW]

Anyone interested to work on this or providing some guidance how to proceed for a fix/patch? For example GNU Octave provides the sos and gain separately from zpk2sos... Also one more question, is there interest for SciPy application level code within SciPy distribution? E.g. a filter design tool?

[Pi9GW]

Rtype ellip provides section gains in inverted order compared to butter.

[roryyorke]

For digital implementation many issues of the analog world disappear or at least change markedly.

I mostly use the signal processing tools for post-processing data, in which case here I use double precision. The range of values is about $10^{-300}$ to $10^{300}$, so analog issues like absolute noise floor or clipping are mostly irrelevant.

guess that is why the gain in first stage is maximized.

It's not maximized - it just is the "accumulated" gain (see my first comment in this thread). For digital filters, it's often very small (see my second comment); my guess is that for a very frequency low-pass filter ( << 1 rad/time unit), analog filters would also have a small gain.

Anyone interested to work on this or providing some guidance how to proceed for a fix/patch?

I'm not a core dev, but I think commits making scipy.signal better for analog filter design would be welcome. The change to zpk2sos to allow for analog filters at all was relatively recent, and I think it would be good to build on that.

I like your idea of normalizing the gain of each section of a SOS at a chosen frequency; it should be straightforward to write a post-processing function to do exactly that. I don't know if it would be accepted in scipy, but the SOS format isn't going to change anytime, so you could keep it as part of your own toolbox.

is there interest for SciPy application level code within SciPy distribution? E.g. a filter design tool?

I don't think Scipy packages any applications.

For discrete time filters there is https://github.com/chipmuenk/pyFDA which builds on scipy and looks very nice, though I've never used it. It looks like it also handles fixed-point filters, and here I think there is a need to distribute gains to avoid overflow and underflow, so maybe it has the tools needed for this sort of thing. It might also be a good place to start with for a more complicated analog filter design tool.

Rtype ellip provides section gains in inverted order compared to butter.

Could you give some detail, e.g., example code and results? All SOS designs use zpk2sos under the hood, so should all have the same behaviour.

[Pi9GW]

Thanks for patience in explaining and advise.

I have prepared an example case for mentioned gain distribution bug and the difference of elliptic case. One can clearly see that all high pass cases are roughly fine, but the only lp filter that is at 1 rad/s works as expected. It also shows elliptic has negative gain at first SOS. It still seems like a bug somewhere. Is there a way to discuss with the author or the maintainer for the background? I would like to have an idea if the analog side has seen much real design usage or has it mostly been a numerical exercise tool, where these issues do not matter at all.

There are two easily expected gain normalization points for analog filters, 0/inf frequency and peak response. Those two would nicely cover most of the cases. Every analog filter designer expects the normalization at 0/inf frequency.

For the filter designer a parametric pole-zero fitter would be a nice addition. The one that allows free number of poles and zeros, both complex and real to be fitted optimally for a specified response mask. Without the limitations of the classical filter functions.

I knew pyFDA, but unfortunately it seems there is no interest for the analog side.

from scipy.signal import butter, ellip, sos2tf, freqs
from scipy import __version__
import matplotlib.pyplot as plt
import numpy as np

if __version__ != '1.10.1':
	print('Warning: SciPy version '+__version__+' does not match the original bug report.')

n = 4 # Plots for orders of 3, 4, 5 and 6 are very informative

fig, axs = plt.subplots(2,4)
axs[0,0].title.set_text('Elliptic LP')
axs[0,1].title.set_text('Butterworth LP')
axs[0,2].title.set_text('Elliptic HP')
axs[0,3].title.set_text('Butterworth HP')

w0 = 100

sos=ellip(n, .1, 60, w0, 'lp', analog=True, output='sos')

for s in sos:
	w, h = freqs(sos2tf([s])[0], sos2tf([s])[1]);axs[0,0].semilogx(w, 20 * np.log10(abs(h)));

sos=butter(n, w0, 'lp', analog=True, output='sos')
for s in sos:
	w, h = freqs(sos2tf([s])[0], sos2tf([s])[1]);axs[0,1].semilogx(w, 20 * np.log10(abs(h)));

sos=ellip(n, .1, 60, w0, 'hp', analog=True, output='sos')
for s in sos:
	w, h = freqs(sos2tf([s])[0], sos2tf([s])[1]);axs[0,2].semilogx(w, 20 * np.log10(abs(h)));

sos=butter(n, w0, 'hp', analog=True, output='sos')
for s in sos:
	w, h = freqs(sos2tf([s])[0], sos2tf([s])[1]);axs[0,3].semilogx(w, 20 * np.log10(abs(h)));

w0 = 1

sos=ellip(n, .1, 60, w0, 'lp', analog=True, output='sos')
for s in sos:
	w, h = freqs(sos2tf([s])[0], sos2tf([s])[1]);axs[1,0].semilogx(w, 20 * np.log10(abs(h)));

sos=butter(n, w0, 'lp', analog=True, output='sos')
for s in sos:
	w, h = freqs(sos2tf([s])[0], sos2tf([s])[1]);axs[1,1].semilogx(w, 20 * np.log10(abs(h)));

sos=ellip(n, .1, 60, w0, 'hp', analog=True, output='sos')
for s in sos:
	w, h = freqs(sos2tf([s])[0], sos2tf([s])[1]);axs[1,2].semilogx(w, 20 * np.log10(abs(h)));

sos=butter(n, w0, 'hp', analog=True, output='sos')
for s in sos:
	w, h = freqs(sos2tf([s])[0], sos2tf([s])[1]);axs[1,3].semilogx(w, 20 * np.log10(abs(h)));

axs[0,0].grid();axs[0,1].grid();axs[0,2].grid();axs[0,3].grid()
axs[1,0].grid();axs[1,1].grid();axs[1,2].grid();axs[1,3].grid()

plt.show()

# ?? Conjugate pole-pair ordering...
# iirdesign output zpk has different conjugate pair ordering than tf2zpk that has neat ordered pairing