New function: nyquist for lti-systems

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.

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):
+        """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)

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()