zpk2sos not working correctly when butterworth band-pass filter is synthesized

[todork]

My issue is about zpk2sos conversion I presume. The power of the numerators of second order sections are evaluated not correctly

Reproducing code example:

Sample code to reproduce the problem

Jupyter QtConsole 5.0.3
Python 3.9.5 (default, May  9 2021, 14:00:28) 
Type 'copyright', 'credits' or 'license' for more information
IPython 7.23.1 -- An enhanced Interactive Python. Type '?' for help.

import math
from scipy import signal
import sympy

s = sympy.Symbol('s')
filterOrder = 3
passBand = list(map(lambda x: x*2*math.pi, [1000, 7000]))
filterType = 'bandpass'
z,p,k = signal.butter( filterOrder, passBand, filterType, analog=True, output='zpk')
b,a = signal.zpk2tf(z,p,k)
sos = signal.zpk2sos(z,p,k)
print("Transfer Function from z,p,k")
n,d = signal.zpk2tf(z,p,k)
print("H(s) =", sympy.Poly(n,s)/sympy.Poly(d,s))
print()
print("Transfer Function from sos")
n1,d1 = signal.sos2tf(sos)
print("H(s) =", sympy.Poly(n1,s)/sympy.Poly(d1,s))
Transfer Function from z,p,k
H(s) = 53578846103558.1*s**3/(1.0*s**6 + 75398.223686155*s**5 + 3671492837.20524*s**4 + 95251281961881.0*s**3 + 1.01461309221017e+18*s**2 + 5.75806638891985e+21*s + 2.11044155773651e+25)

Transfer Function from sos
H(s) = 53578846103558.1*s**6/(1.0*s**6 + 75398.223686155*s**5 + 3671492837.20524*s**4 + 95251281961881.0*s**3 + 1.01461309221017e+18*s**2 + 5.75806638891985e+21*s + 2.11044155773651e+25)

It can be seen that the numerators are different in the power of s.

Scipy/Numpy/Python version information:

1.6.3 1.20.3 sys.version_info(major=3, minor=9, micro=5, releaselevel='final', serial=0)

[ilayn]

Hi @todork do you mind evaluating the results over a signal or their impulse response. Changing to transfer function of such filters is desperately ill-conditioned problems and also I wouldn't say we have the most robust conversion algorithm in place. So if you can avoid working in TF form anyways but missing 3 zeros should definitely make a difference if they are indeed inherently different.

[todork]

Hi @ilayn thank you for you prompt reply.

I work in SOS area

sos = signal.butter( filterOrder, passBand, filterType, analog=True, output='sos')

and I made the above example just to show where I see the problem. When we use Butterworth approximation of third order we should have three zeros and the obtained SOS biquads seems to have six.

sos = array([ [5.35788461e+13, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 3.25936385e+04, 1.76422758e+09], [1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 3.76991118e+04, 2.76348923e+08], [1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 5.10547335e+03, 4.32873446e+07] ])

I would expect

sos = array([ [0.00000000e+00, 5.35788461e+13, 0.00000000e+00, 1.00000000e+00, 3.25936385e+04, 1.76422758e+09], [0.00000000e+00, 1.00000000e+00, 0.00000000e+00, 1.00000000e+00, 3.76991118e+04, 2.76348923e+08], [0.00000000e+00, 1.00000000e+00, 0.00000000e+00, 1.00000000e+00, 5.10547335e+03, 4.32873446e+07] ])

Or may be I'm wrong?

Regards, Todor

[ilayn]

Ah OK. My bad. I wasn't clear enough. I was trying to convey the message that maybe conversion to transfer function was the faulty one but I'll try to check. Thanks again.

[roryyorke]

This is likely a duplicate of https://github.com/scipy/scipy/issues/5668

Out of curiosity, why do you want SOS for analog filters?

[todork]

Hi @roryyorke,

Yes the idea behind is the same. The issue was opened in 2016, is there a chance to solve it?

I would like to realize RC biquad sections with op amps.

Regards, T.

[roryyorke]

I won't be able to look at this soon, sorry, but maybe someone else can. It has come up quite a few times, so we must fix it.

Thanks for the info on using this for the op-amp filter design.

[roryyorke]

Closed via gh-15085.