search

scipy / scipy

SciPy library main repository
https://scipy.org
BSD 3-Clause "New" or "Revised" License
13.08k stars 5.19k forks source link

scipy.signal.bode doesn't work for state space (A, B, C, D) argument (Trac #1863) #2382

Closed scipy-gitbot closed 11 years ago

scipy-gitbot commented 11 years ago

Original ticket http://projects.scipy.org/scipy/ticket/1863 on 2013-03-11 by @hazelnusse, assigned to @cournape.

See attached for an example of what fails. This is the error I get:

/home/hazelnusse/usr/lib64/python3.2/site-packages/scipy/signal/filter_design.py:306: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless "results may be meaningless", BadCoefficients) Traceback (most recent call last): File "test_bode.py", line 19, in bode((A, B, C, D)) File "/home/hazelnusse/usr/lib64/python3.2/site-packages/scipy/signal/ltisys.py", line 944, in bode y = numpy.polyval(sys.num, jw) / numpy.polyval(sys.den, jw) File "/home/hazelnusse/usr/lib64/python3.2/site-packages/numpy/lib/polynomial.py", line 671, in polyval y = x * y + p[i] ValueError: operands could not be broadcast together with shapes (100) (3)

scipy-gitbot commented 11 years ago

Attachment added by @hazelnusse on 2013-03-11: test_bode.py

scipy-gitbot commented 11 years ago

@hazelnusse wrote on 2013-03-13

PR gh-959 (https://github.com/scipy/scipy/pull/432) was recently mereged but this bug still remains.

hazelnusse commented 11 years ago

@WarrenWeckesser @bjornfor: I just wanted to give you guys a heads up that this issue is still arround. Am I using signal.bode() incorrectly or is it just that the state space interface hasn't been tested as much as the transfer function interface?

bjornfor commented 11 years ago

Hm, maybe I never tested state space? Not sure what happens here.

/home/bfo/.local/lib/python2.7/site-packages/scipy/signal/ltisys.pyc in bode(system, w, n)
    927 
    928     jw = w * 1j
--> 929     y = numpy.polyval(sys.num, jw) / numpy.polyval(sys.den, jw)
    930     mag = 20.0 * numpy.log10(abs(y))
    931     phase = numpy.arctan2(y.imag, y.real) * 180.0 / numpy.pi

/usr/lib/python2.7/dist-packages/numpy/lib/polynomial.pyc in polyval(p, x)
    640         y = NX.zeros_like(x)
    641     for i in range(len(p)):
--> 642         y = x * y + p[i]
    643     return y
    644 

ValueError: operands could not be broadcast together with shapes (100) (3) 
> /usr/lib/python2.7/dist-packages/numpy/lib/polynomial.py(642)polyval()
    641     for i in range(len(p)):
--> 642         y = x * y + p[i]
    643     return y
ipdb> p x * y
array([ 0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j])
ipdb> p p[i]
array([0, 0, 1])

An array of shape 100x1 doesn't add with an array of shape 3x1...

WarrenWeckesser commented 11 years ago

@hazelnusse: Thanks for reporting the problem. As you pointed out in https://github.com/scipy/scipy/pull/432, the problem is that when an lti instance is created from state space matrices, sys.num is a 2D array containing a single row. I would like to fix that problem, too, but that has been the behavior for a long time, and I don't know if there would be a reasonable way to deprecate it.

I have submitted a pull request that fixes bode (and freqresp,which has the same problem) here: https://github.com/scipy/scipy/issues/2382

rgommers commented 11 years ago

gh-2382 merged.

carloshernangarrido commented 2 years ago

I solved the problem by changing:

sys = sp.signal.lti(A, B, C, D) w, mag, _ = sys.bode(w=1000)

by

eps = 1e-6 num, den = sp.signal.ss2tf(A, B, C, D) num = [coef for coef in num[0] if (abs(coef) > eps)] sys = sp.signal.lti(num, den) w, mag, _ = sys.bode(w=1000)

i.e. you have to delete near zero coefficients.