Skip to content

Loading…

New function: nyquist for lti-systems #418

Merged
merged 9 commits into from

7 participants

@Niklatz

First time participating to a open source project. I hope I didn't do too many mistakes.

In the scipy.signal module I missed the ability to generate a Nyquist plot, so I made a quick implementation and figured it might be useful to others.

Code based on the lti.bode() function.
n=10000 decided through trial and error. Might be somewhat high for simple systems.

Like I said, it's the first time I participate in a open source project and I don't really know if I made any mistakes. I hope the code is useful even if it doesn't make it into scipy. At least I tried.

As attachment the plot that is given out by the example code in the comment:
figure_1

@Niklatz Niklatz Added function to determine Nyquist plot
Code based on the lti.bode() function.
n=10000 decided through trial and error. Might be somewhat high for
simple systems.
6608fa8
@rgommers
SciPy member

Hi, welcome. Code looks pretty clean and docs look good. What is missing is unit tests for the function - have a look in scipy/signal/tests/test_ltisys.py how that's done for other functions in this module.

Whether this function should be added here and if the API is optimal I don't know. Probably good to ask that on the scipy-dev mailing list.

@ewmoore
SciPy member

Why return complex data as two parts instead of one?

Even though you won't see it if you make a plot, the function probably ought to return the frequency array too. A bode member function is part of the lti class, perhaps this should be added there too?

If this function is made to return the frequency array, bode could become a simple wrapper around it.. I don't know that it's really necessary though, as the amount of duplication is quite small.

@Niklatz Niklatz improved nyquist() function
added testing and changed return values as suggested by rgommers and
ewmoore
362790b
@pv
SciPy member
pv commented

From non-expert POV, looks good.

@pierre-haessig pierre-haessig commented on an outdated diff
scipy/signal/ltisys.py
@@ -930,3 +933,61 @@ def bode(system, w=None, n=100):
mag = 20.0 * numpy.log10(abs(y))
phase = numpy.arctan2(y.imag, y.real) * 180.0 / numpy.pi
return w, mag, phase
+
+
+def nyquist(system, w=None, n=10000):

As I just posted on the ML, I would rather call this function freqresp or some similar name conveying the idea of frequency response. I would change the docstring accordingly by renaming the nyq variable into H (or G since G(s) is mentionned in the Example)

reference in Matlab documentation : http://www.mathworks.fr/fr/help/control/ref/freqresp.html

Otherwise, looks good to me !

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@Niklatz Niklatz changed name to freqresp()
Changed the name as suggested on the Mailing List and adjusted the
docstrings
be4697b
@WarrenWeckesser

Are you aware of scipy.signal.freqs? It already provides almost the same functionality as freqresp, and includes code that attempts to figure out a reasonable frequency range if no frequencies are given.

This snippet produces the same plot as the example in the freqresp docstring:

import numpy as np
from scipy.signal import lti, freqs
import matplotlib.pyplot as plt

system = lti([], [1,1,1], [5])
w, h = freqs(system.num, system.den)
plt.plot(h.real, h.imag, 'b')
plt.plot(h.real, -h.imag, 'r')
plt.grid()
plt.show()

Perhaps it would be better to add a freqresp method to the lti class, and implement it using freqs.

@Niklatz

Well, I was not aware of that function but thanks for pointing to it.
I just had a look into the code of freqs and it seems like it uses the exact same code I do... Seems a bit redundant to keep the same code in the scipy package twice. The only difference I see is the fact that freqresp takes lti-systems and does not require you to pass numerator and denumerator.
If you consider it better practice that i implement freqresp by calling freqs that's no problem.

@WarrenWeckesser

Sorry, I just noticed that this PR already adds freqresp as a method of lti.

It feels like too much duplication to have both freqs and freqresp. Anyone see a clean, backwards-compatible change to the API that would allow freqs to also accept an lti instance? Then bode should follow the same convention. (Welcome to the joys of maintaining a clean, consistent API in a long-term project with many authors :)

@WarrenWeckesser

I thought a bit more about it. Perhaps maintaining freqs and freqresp isn't too bad. But, to eliminate duplicated code, I think freqresp should be implemented with freqs, and bode should be implemented with freqresp (or freqs).

Maybe something like:

def freqresp(system, w=None, n=None):
    worN = w
    if w is None:
        worN = n
    return freqs(system.num, system.den, worN=worN)

Implementing bode using freqresp or freqs would mean eliminating the duplicated functionality provided by filter_design.findfreqs (used by freqs) and ltisys._default_response_frequencies (used by bode). They don't behave exactly the same, so either freqs or bode would end up with a slightly different default frequency range.

@Niklatz

I don't really know how to decide what the frequency range for these functions should be, so I can't decide which one should be used. I can only say that in my opinion _default_response_frequenciesused by bode is easier to understand and better documented.
I just tested both of them however and they both delivered useful frequency ranges to work with (I generated a quick bode plot for some transfer functions to test this).

@bjornfor

Hi, I'm the guy that wrote _default_response_frequencies. I tried to make it behave exactly like Matlab bode. After lots of trial and error think I got it right; the LTI systems I've thrown at it get the same frequency range as Matlab. I didn't know about filter_design.findfreqs until now.

@WarrenWeckesser WarrenWeckesser commented on the diff
scipy/signal/ltisys.py
@@ -362,6 +362,9 @@ def bode(self, w=None, n=100):
"""
return bode(self, w=w, n=n)
+ def freqresp(self, w=None, n=10000):

This method needs a docstring. At a minimum, something analogous to the docstring of the method bode.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@WarrenWeckesser

@Niklatz, @bjornfor: I agree that _default_response_frequencies is much easier to read, and it works well for real poles. However, it ignores the complex parts of the poles, and it ignores the zeros entirely, so it can make some bad choices for the frequency range. I'll create a new pull request to address this soon.

Because of that, I think freqresp should use scipy.signal.freqs as I suggested above, but with the added code to check for system being a tuple rather than an lti instance (as it and bode already do).

Also, I noted inline that lti.freqresp needs a docstring.

@WarrenWeckesser

@Niklatz: P.S.: I should have said this in my first comment: thank you for the contribution! It is great to have more people working on scipy.

One thing you should do differently in future pull requests is to create a branch in your git repo in which to work, and submit that branch as the pull request. This will make it easier to rebase, or delete the branch, or work on other features while the pull request is being discussed.

@Niklatz Niklatz closed this
@Niklatz Niklatz reopened this
@Niklatz

whoops, clicked there by accident. Didn't mean to close the request, sorry about that.

@WarrenWeckesser Thanks for the hint about working with git. First time I do this so this is very helpful.

About the docstring: Do I really need to add another one in line 365 or is the one in line 939 in the function definition enough? I was confused why bode has one in both places. Is there a reason behind this?

About implementing freqresp with the help of scipy.signal.freqs, I will get to that soon and push the changes when I'm done.

-Niklas.

@WarrenWeckesser

@Niklatz: Yes, lti.freqresp() needs a docstring. A user who is exploring the lti class might not even know that the corresponding function exists. It is frustrating to be in ipython and do something like:

In [6]: s = lti([],[1,1,1],[5])

In [7]: s.freqresp?
Type:       instancemethod
String Form:<bound method lti.freqresp of <scipy.signal.ltisys.lti object at 0x4660250>>
File:       /home/warren/local_scipy/lib/python2.7/site-packages/scipy/signal/ltisys.py
Definition: s.freqresp(self, w=None, n=10000)
Docstring:  <no docstring>

What!? No docstring!? $%@$$@*&#$!

@Niklatz

@WarrenWeckesser: convincing explanation. I will add a docstring.

@pv
SciPy member
pv commented

Please also add it to the list in the docstring of scipy/signal/__init__.py --- so that it will be included in the online documentation.

@WarrenWeckesser

@Niklatz, @bjornfor, and anyone else interested: see #432

Niklatz added some commits
@Niklatz Niklatz freqresp() now uses freqs()
Also changed order of return values to match that of freqs() to keep it
consistent
342fb51
@Niklatz Niklatz fixed error in docstring
forgot to adjust the example to the changed order of return values

Also i just saw, that bode has already been added to __init__.py so I
took that back out... Sorry
8d969cd
@Niklatz Niklatz fixed testing for frequency range
Also adjusted docstring and examples to have consistent names with
freqs()
022ad07
@Niklatz

Alright, I hope I've got everything. Thank you all for your patience with me. If there are still things missing/wrong, just leave a comment.

@WarrenWeckesser WarrenWeckesser commented on an outdated diff
scipy/signal/ltisys.py
@@ -19,11 +19,12 @@
import scipy.integrate as integrate
import scipy.linalg as linalg
from scipy.lib.six.moves import xrange
+from scipy.signal.filter_design import freqs

Instead of this, just add freqs to the existing import from .filter_design in line 14 above.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@WarrenWeckesser WarrenWeckesser commented on an outdated diff
scipy/signal/ltisys.py
((5 lines not shown))
from numpy import r_, eye, real, atleast_1d, atleast_2d, poly, \
squeeze, diag, asarray
__all__ = ['tf2ss', 'ss2tf', 'abcd_normalize', 'zpk2ss', 'ss2zpk', 'lti',
- 'lsim', 'lsim2', 'impulse', 'impulse2', 'step', 'step2', 'bode']
+ 'lsim', 'lsim2', 'impulse', 'impulse2', 'step', 'step2', 'bode', 'freqresp']

Put 'freqresp' in a new line, to conform to PEP-8 (maximum of 79 characters in a line).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@WarrenWeckesser WarrenWeckesser commented on an outdated diff
scipy/signal/ltisys.py
@@ -362,6 +363,15 @@ def bode(self, w=None, n=100):
"""
return bode(self, w=w, n=n)
+ def freqresp(self, w=None, n=10000):
+ """Calculate the frequency response of a continuous-time system.
+
+ Returns a 2-tuple containing arrays of frequencies [rad/s] and complex magnitude.

PEP8: line is too long and has trailiing whitespace.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@WarrenWeckesser WarrenWeckesser commented on an outdated diff
scipy/signal/ltisys.py
((28 lines not shown))
+
+ Returns
+ -------
+ w : 1D ndarray
+ Frequency array [rad/s]
+ H : 1D ndarray
+ Array of complex magnitude values
+
+ Example
+ -------
+ # Generating the Nyquist plot of a transfer function
+
+ >>> from scipy import signal
+ >>> import matplotlib.pyplot as plt
+
+ >>> s1 = signal.lti([], [1, 1, 1], [5]) #transfer function: H(s) = 5 / (s-1)^3

PEP8: line too long. Also, put a space after the comment character.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@WarrenWeckesser WarrenWeckesser commented on an outdated diff
scipy/signal/tests/test_ltisys.py
@@ -5,7 +5,7 @@
import numpy as np
from numpy.testing import assert_almost_equal, assert_equal, run_module_suite
-from scipy.signal.ltisys import ss2tf, lsim2, impulse2, step2, lti, bode
+from scipy.signal.ltisys import ss2tf, lsim2, impulse2, step2, lti, bode, freqresp

PEP8: line too long.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@WarrenWeckesser WarrenWeckesser commented on an outdated diff
scipy/signal/tests/test_ltisys.py
((46 lines not shown))
+ w = [0.1, 1, 10, 100]
+ w, H = freqresp(system, w=w)
+ jw = w * 1j
+ y = np.polyval(system.num, jw) / np.polyval(system.den, jw)
+ expected_im = y.imag
+ assert_almost_equal(H.imag, expected_im)
+
+ def test_05(self):
+ """Test that freqresp() finds a reasonable frequency range."""
+ # 1st order low-pass filter: H(s) = 1 / (s + 1)
+ # Expected range is from 0.01 to 10.
+ system = lti([1], [1, 1])
+ vals = linalg.eigvals(system.A)
+ minpole = min(abs(np.real(vals)))
+ maxpole = max(abs(np.real(vals)))
+ n = 10;

Remove the semicolon.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@Niklatz Niklatz PEP-8 conformity
Thank you WarrenWeckesser for pointing those out.
33869b5
@WarrenWeckesser WarrenWeckesser commented on an outdated diff
scipy/signal/tests/test_ltisys.py
((6 lines not shown))
+
+ def test_01(self):
+ """Test freqresp() real part calculation (manual sanity check)."""
+ # 1st order low-pass filter: H(s) = 1 / (s + 1),
+ # re(H(s=0.1)) ~= 0.99
+ # re(H(s=1)) ~= 0.5
+ # re(H(s=10)) ~= 0.0099
+ system = lti([1], [1, 1])
+ w = [0.1, 1, 10]
+ w, H = freqresp(system, w=w)
+ expected_re = [0.99, 0.5, 0.0099]
+ assert_almost_equal(H.real, expected_re, decimal=1)
+
+ def test_02(self):
+ """Test freqresp() imaginary part calculation (manual sanity
+ check)."""

This breaks the single line summary into two lines, which is not ideal. A common variation that you could use here is

"""
Test freqresp() imaginary part calculation (manual sanity check).
"""
@rgommers SciPy member

Tests shouldn't have docstrings at all. Use descriptive names, and comments where needed.

@Niklatz
Niklatz added a note

So should I just change the docstrings to simple comment lines?

@rgommers SciPy member

yes please - that makes it quite a bit easier to see what test is failing.

@Niklatz
Niklatz added a note

Alright, will do that.

Heh, learn something new every day. I have been annoyed about not seeing the name of the function when an error occurred in a function with a docstring, but I didn't realize it was policy to not have docstrings.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@WarrenWeckesser

Thanks for the quick update! I just made one more style suggestion. Other than that, looks good to me.

@rgommers
SciPy member

@WarrenWeckesser please merge this today if the last commits look good as well - tomorrow I'll branch 0.12.x

@WarrenWeckesser WarrenWeckesser merged commit 56c3e4b into scipy:master
@WarrenWeckesser

Merged. Thanks Niklas!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Feb 5, 2013
  1. @Niklatz

    Added function to determine Nyquist plot

    Niklatz committed
    Code based on the lti.bode() function.
    n=10000 decided through trial and error. Might be somewhat high for
    simple systems.
Commits on Feb 6, 2013
  1. @Niklatz

    improved nyquist() function

    Niklatz committed
    added testing and changed return values as suggested by rgommers and
    ewmoore
Commits on Feb 8, 2013
  1. @Niklatz

    changed name to freqresp()

    Niklatz committed
    Changed the name as suggested on the Mailing List and adjusted the
    docstrings
Commits on Feb 12, 2013
  1. @Niklatz

    freqresp() now uses freqs()

    Niklatz committed
    Also changed order of return values to match that of freqs() to keep it
    consistent
  2. @Niklatz

    fixed error in docstring

    Niklatz committed
    forgot to adjust the example to the changed order of return values
    
    Also i just saw, that bode has already been added to __init__.py so I
    took that back out... Sorry
  3. @Niklatz

    fixed testing for frequency range

    Niklatz committed
    Also adjusted docstring and examples to have consistent names with
    freqs()
  4. @Niklatz

    PEP-8 conformity

    Niklatz committed
    Thank you WarrenWeckesser for pointing those out.
  5. @Niklatz

    PEP-8 Style change

    Niklatz committed
  6. @Niklatz

    removed docstring from tests

    Niklatz committed
    As asked by rgommers
This page is out of date. Refresh to see the latest.
Showing with 146 additions and 4 deletions.
  1. +1 −0 scipy/signal/__init__.py
  2. +74 −3 scipy/signal/ltisys.py
  3. +71 −1 scipy/signal/tests/test_ltisys.py
View
1 scipy/signal/__init__.py
@@ -116,6 +116,7 @@
.. autosummary::
:toctree: generated/
+ freqresp -- frequency response of a continuous-time LTI system.
lti -- linear time invariant system object.
lsim -- continuous-time simulation of output to linear system.
lsim2 -- like lsim, but `scipy.integrate.odeint` is used.
View
77 scipy/signal/ltisys.py
@@ -11,7 +11,7 @@
# Rewrote lsim2 and added impulse2.
#
-from .filter_design import tf2zpk, zpk2tf, normalize
+from .filter_design import tf2zpk, zpk2tf, normalize, freqs
import numpy
from numpy import product, zeros, array, dot, transpose, ones, \
nan_to_num, zeros_like, linspace
@@ -23,7 +23,8 @@
squeeze, diag, asarray
__all__ = ['tf2ss', 'ss2tf', 'abcd_normalize', 'zpk2ss', 'ss2zpk', 'lti',
- 'lsim', 'lsim2', 'impulse', 'impulse2', 'step', 'step2', 'bode']
+ 'lsim', 'lsim2', 'impulse', 'impulse2', 'step', 'step2', 'bode',
+ 'freqresp']
def tf2ss(num, den):
@@ -362,6 +363,16 @@ def bode(self, w=None, n=100):
"""
return bode(self, w=w, n=n)
+ def freqresp(self, w=None, n=10000):

This method needs a docstring. At a minimum, something analogous to the docstring of the method bode.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ """Calculate the frequency response of a continuous-time system.
+
+ Returns a 2-tuple containing arrays of frequencies [rad/s] and
+ complex magnitude.
+ See scipy.signal.freqresp for details.
+
+ """
+ return freqresp(self, w=w, n=n)
+
def lsim2(system, U=None, T=None, X0=None, **kwargs):
"""
@@ -611,7 +622,7 @@ def _default_response_frequencies(A, n):
A : ndarray
The system matrix, which is square.
n : int
- The number of time samples to generate.
+ The number of frequency samples to generate.
Returns
-------
@@ -930,3 +941,63 @@ def bode(system, w=None, n=100):
mag = 20.0 * numpy.log10(abs(y))
phase = numpy.arctan2(y.imag, y.real) * 180.0 / numpy.pi
return w, mag, phase
+
+
+def freqresp(system, w=None, n=10000):
+ """Calculate the frequency response of a continuous-time system.
+
+ Parameters
+ ----------
+ system : an instance of the LTI class or a tuple describing the system.
+ The following gives the number of elements in the tuple and
+ the interpretation:
+
+ * 2 (num, den)
+ * 3 (zeros, poles, gain)
+ * 4 (A, B, C, D)
+
+ w : array_like, optional
+ Array of frequencies (in rad/s). Magnitude and phase data is
+ calculated for every value in this array. If not given a reasonable
+ set will be calculated.
+ n : int, optional
+ Number of frequency points to compute if `w` is not given. The `n`
+ frequencies are logarithmically spaced in the range from two orders of
+ magnitude before the minimum (slowest) pole to two orders of magnitude
+ after the maximum (fastest) pole.
+
+ Returns
+ -------
+ w : 1D ndarray
+ Frequency array [rad/s]
+ H : 1D ndarray
+ Array of complex magnitude values
+
+ Example
+ -------
+ # Generating the Nyquist plot of a transfer function
+
+ >>> from scipy import signal
+ >>> import matplotlib.pyplot as plt
+
+ >>> s1 = signal.lti([], [1, 1, 1], [5])
+ # transfer function: H(s) = 5 / (s-1)^3
+
+ >>> w, H = signal.freqresp(s1)
+
+ >>> plt.figure()
+ >>> plt.plot(H.real, H.imag, "b")
+ >>> plt.plot(H.real, -H.imag, "r")
+ >>> plt.show()
+ """
+ if isinstance(system, lti):
+ sys = system
+ else:
+ sys = lti(*system)
+
+ if w is not None:
+ worN = w
+ else:
+ worN = n
+
+ return freqs(sys.num, sys.den, worN=worN)
View
72 scipy/signal/tests/test_ltisys.py
@@ -5,7 +5,8 @@
import numpy as np
from numpy.testing import assert_almost_equal, assert_equal, run_module_suite
-from scipy.signal.ltisys import ss2tf, lsim2, impulse2, step2, lti, bode
+from scipy.signal.ltisys import ss2tf, lsim2, impulse2, step2, lti, bode, \
+ freqresp
from scipy.signal.filter_design import BadCoefficients
import scipy.linalg as linalg
@@ -307,5 +308,74 @@ def test_06(self):
w, mag, phase = bode(system, n=2)
assert_equal(w[0], 0.01) # a fail would give not-a-number
+
+class Test_freqresp(object):
+
+ def test_real_part_manual(self):
+ # Test freqresp() real part calculation (manual sanity check).
+ # 1st order low-pass filter: H(s) = 1 / (s + 1),
+ # re(H(s=0.1)) ~= 0.99
+ # re(H(s=1)) ~= 0.5
+ # re(H(s=10)) ~= 0.0099
+ system = lti([1], [1, 1])
+ w = [0.1, 1, 10]
+ w, H = freqresp(system, w=w)
+ expected_re = [0.99, 0.5, 0.0099]
+ assert_almost_equal(H.real, expected_re, decimal=1)
+
+ def test_imag_part_manual(self):
+ # Test freqresp() imaginary part calculation (manual sanity check).
+ # 1st order low-pass filter: H(s) = 1 / (s + 1),
+ # im(H(s=0.1)) ~= -0.099
+ # im(H(s=1)) ~= -0.5
+ # im(H(s=10)) ~= -0.099
+ system = lti([1], [1, 1])
+ w = [0.1, 1, 10]
+ w, H = freqresp(system, w=w)
+ expected_im = [-0.099, -0.5, -0.099]
+ assert_almost_equal(H.imag, expected_im, decimal=1)
+
+ def test_real_part(self):
+ # Test freqresp() real part calculation.
+ # 1st order low-pass filter: H(s) = 1 / (s + 1)
+ system = lti([1], [1, 1])
+ w = [0.1, 1, 10, 100]
+ w, H = freqresp(system, w=w)
+ jw = w * 1j
+ y = np.polyval(system.num, jw) / np.polyval(system.den, jw)
+ expected_re = y.real
+ assert_almost_equal(H.real, expected_re)
+
+ def test_imag_part(self):
+ # Test freqresp() imaginary part calculation.
+ # 1st order low-pass filter: H(s) = 1 / (s + 1)
+ system = lti([1], [1, 1])
+ w = [0.1, 1, 10, 100]
+ w, H = freqresp(system, w=w)
+ jw = w * 1j
+ y = np.polyval(system.num, jw) / np.polyval(system.den, jw)
+ expected_im = y.imag
+ assert_almost_equal(H.imag, expected_im)
+
+ def test_freq_range(self):
+ # Test that freqresp() finds a reasonable frequency range.
+ # 1st order low-pass filter: H(s) = 1 / (s + 1)
+ # Expected range is from 0.01 to 10.
+ system = lti([1], [1, 1])
+ vals = linalg.eigvals(system.A)
+ minpole = min(abs(np.real(vals)))
+ maxpole = max(abs(np.real(vals)))
+ n = 10
+ expected_w = np.logspace(-2, 1, n)
+ w, H = freqresp(system, n=n)
+ assert_almost_equal(w, expected_w)
+
+ def test_pole_zero(self):
+ # Test that freqresp() doesn't fail on a system with a pole at 0.
+ # integrator, pole at zero: H(s) = 1 / s
+ system = lti([1], [1, 0])
+ w, H = freqresp(system, n=2)
+ assert_equal(w[0], 0.01) # a fail would give not-a-number
+
if __name__ == "__main__":
run_module_suite()
Something went wrong with that request. Please try again.