python-control · murrayrm · Jan 8, 2021 · Aug 14, 2020 · Aug 16, 2020 · Aug 16, 2020
diff --git a/control/bench/time_freqresp.py b/control/bench/time_freqresp.py
@@ -8,7 +8,7 @@
 sys_tf = tf(sys)
 w = logspace(-1,1,50)
 ntimes = 1000
-time_ss = timeit("sys.freqresp(w)", setup="from __main__ import sys, w", number=ntimes)
-time_tf = timeit("sys_tf.freqresp(w)", setup="from __main__ import sys_tf, w", number=ntimes)
+time_ss = timeit("sys.freqquency_response(w)", setup="from __main__ import sys, w", number=ntimes)
+time_tf = timeit("sys_tf.frequency_response(w)", setup="from __main__ import sys_tf, w", number=ntimes)
 print("State-space model on %d states: %f" % (nstates, time_ss))
 print("Transfer-function model on %d states: %f" % (nstates, time_tf))
diff --git a/control/frdata.py b/control/frdata.py
@@ -48,7 +48,7 @@
 from warnings import warn
 import numpy as np
 from numpy import angle, array, empty, ones, \
-    real, imag, absolute, eye, linalg, where, dot
+    real, imag, absolute, eye, linalg, where, dot, sort
 from scipy.interpolate import splprep, splev
 from .lti import LTI
 
@@ -100,24 +100,18 @@ def __init__(self, *args, **kwargs):
         object, other than an FRD, call FRD(sys, omega)
 
         """
+        # TODO: discrete-time FRD systems?
         smooth = kwargs.get('smooth', False)
 
         if len(args) == 2:
             if not isinstance(args[0], FRD) and isinstance(args[0], LTI):
                 # not an FRD, but still a system, second argument should be
                 # the frequency range
                 otherlti = args[0]
-                self.omega = array(args[1], dtype=float)
-                self.omega.sort()
+                self.omega = sort(np.asarray(args[1], dtype=float))
                 numfreq = len(self.omega)
-
                 # calculate frequency response at my points
-                self.fresp = empty(
-                    (otherlti.outputs, otherlti.inputs, numfreq),
-                    dtype=complex)
-                for k, w in enumerate(self.omega):
-                    self.fresp[:, :, k] = otherlti._evalfr(w)
-
+                self.fresp = otherlti(1j * self.omega, squeeze=False)
             else:
                 # The user provided a response and a freq vector
                 self.fresp = array(args[0], dtype=complex)
@@ -141,7 +135,7 @@ def __init__(self, *args, **kwargs):
             self.omega = args[0].omega
             self.fresp = args[0].fresp
         else:
-            raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
+            raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
         # create interpolation functions
         if smooth:
@@ -163,17 +157,17 @@ def __str__(self):
         mimo = self.inputs > 1 or self.outputs > 1
         outstr = ['Frequency response data']
 
-        mt, pt, wt = self.freqresp(self.omega)
         for i in range(self.inputs):
             for j in range(self.outputs):
                 if mimo:
                     outstr.append("Input %i to output %i:" % (i + 1, j + 1))
                 outstr.append('Freq [rad/s]  Response')
                 outstr.append('------------  ---------------------')
                 outstr.extend(
-                    ['%12.3f  %10.4g%+10.4gj' % (w, m, p)
-                     for m, p, w in zip(real(self.fresp[j, i, :]),
-                                        imag(self.fresp[j, i, :]), wt)])
+                    ['%12.3f  %10.4g%+10.4gj' % (w, re, im)
+                     for w, re, im in zip(self.omega,
+                                          real(self.fresp[j, i, :]),
+                                          imag(self.fresp[j, i, :]))])
 
         return '\n'.join(outstr)
 
@@ -342,110 +336,116 @@ def __pow__(self, other):
             return (FRD(ones(self.fresp.shape), self.omega) / self) * \
                 (self**(other+1))
 
-    def evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self._evalfr(omega) returns the value of the frequency response
-        at frequency omega.
-
-        Note that a "normal" FRD only returns values for which there is an
-        entry in the omega vector. An interpolating FRD can return
-        intermediate values.
-
-        """
-        warn("FRD.evalfr(omega) will be deprecated in a future release "
-             "of python-control; use sys.eval(omega) instead",
-             PendingDeprecationWarning)         # pragma: no coverage
-        return self._evalfr(omega)
-
     # Define the `eval` function to evaluate an FRD at a given (real)
     # frequency.  Note that we choose to use `eval` instead of `evalfr` to
     # avoid confusion with :func:`evalfr`, which takes a complex number as its
     # argument.  Similarly, we don't use `__call__` to avoid confusion between
     # G(s) for a transfer function and G(omega) for an FRD object.
-    def eval(self, omega):
-        """Evaluate a transfer function at a single angular frequency.
-
-        self.evalfr(omega) returns the value of the frequency response
-        at frequency omega.
+    # update Sawyer B. Fuller 2020.08.14: __call__ added to provide a uniform
+    # interface to systems in general and the lti.frequency_response method
+    def eval(self, omega, squeeze=True):
+        """Evaluate a transfer function at angular frequency omega.
 
         Note that a "normal" FRD only returns values for which there is an
         entry in the omega vector. An interpolating FRD can return
         intermediate values.
 
+        Parameters
+        ----------
+        omega : float or array_like
+            Frequencies in radians per second
+        squeeze : bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), returns a
+            1D array rather than a 3D array.
+
+        Returns
+        -------
+        fresp : (self.outputs, self.inputs, len(x)) or (len(x), ) complex ndarray
+            The frequency response of the system. Array is ``len(x)`` if and only
+            if system is SISO and ``squeeze=True``.
         """
-        return self._evalfr(omega)
-
-    # Internal function to evaluate the frequency responses
-    def _evalfr(self, omega):
-        """Evaluate a transfer function at a single angular frequency."""
-        # Preallocate the output.
-        if getattr(omega, '__iter__', False):
-            out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
-        else:
-            out = empty((self.outputs, self.inputs), dtype=complex)
+        omega_array = np.array(omega, ndmin=1) # array-like version of omega
+        if any(omega_array.imag > 0):
+            raise ValueError("FRD.eval can only accept real-valued omega")
 
         if self.ifunc is None:
-            try:
-                out = self.fresp[:, :, where(self.omega == omega)[0][0]]
-            except Exception:
+            elements = np.isin(self.omega, omega) # binary array
+            if sum(elements) < len(omega_array):
                 raise ValueError(
-                    "Frequency %f not in frequency list, try an interpolating"
-                    " FRD if you want additional points" % omega)
-        else:
-            if getattr(omega, '__iter__', False):
-                for i in range(self.outputs):
-                    for j in range(self.inputs):
-                        for k, w in enumerate(omega):
-                            frraw = splev(w, self.ifunc[i, j], der=0)
-                            out[i, j, k] = frraw[0] + 1.0j * frraw[1]
+                    "not all frequencies omega are in frequency list of FRD "
+                    "system. Try an interpolating FRD for additional points.")
             else:
-                for i in range(self.outputs):
-                    for j in range(self.inputs):
-                        frraw = splev(omega, self.ifunc[i, j], der=0)
-                        out[i, j] = frraw[0] + 1.0j * frraw[1]
-
+                out = self.fresp[:, :, elements]
+        else:
+            out = empty((self.outputs, self.inputs, len(omega_array)),
+                         dtype=complex)
+            for i in range(self.outputs):
+                for j in range(self.inputs):
+                    for k, w in enumerate(omega_array):
+                        frraw = splev(w, self.ifunc[i, j], der=0)
+                        out[i, j, k] = frraw[0] + 1.0j * frraw[1]
+        if not hasattr(omega, '__len__'):
+            # omega is a scalar, squeeze down array along last dim
+            out = np.squeeze(out, axis=2)
+        if squeeze and self.issiso():
+            out = out[0][0]
         return out
 
-    # Method for generating the frequency response of the system
-    def freqresp(self, omega):
-        """Evaluate the frequency response at a list of angular frequencies.
+    def __call__(self, s, squeeze=True):
+        """Evaluate system's transfer function at complex frequencies.
+
+        Returns the complex frequency response `sys(s)` of system `sys` with
+        `m = sys.inputs` number of inputs and `p = sys.outputs` number of
+        outputs.
 
-        Reports the value of the magnitude, phase, and angular frequency of
-        the requency response evaluated at omega, where omega is a list of
-        angular frequencies, and is a sorted version of the input omega.
+        To evaluate at a frequency omega in radians per second, enter
+        ``s = omega * 1j`` or use ``sys.eval(omega)``
 
         Parameters
         ----------
-        omega : array_like
-            A list of frequencies in radians/sec at which the system should be
-            evaluated. The list can be either a python list or a numpy array
-            and will be sorted before evaluation.
+        s : complex scalar or array_like
+            Complex frequencies
+        squeeze : bool, optional (default=True)
+            If True and `sys` is single input single output (SISO), i.e. `m=1`,
+            `p=1`, return a 1D array rather than a 3D array.
 
         Returns
         -------
-        mag : (self.outputs, self.inputs, len(omega)) ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.
-        phase : (self.outputs, self.inputs, len(omega)) ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray or list or tuple
-            The list of sorted frequencies at which the response was
-            evaluated.
-        """
-        # Preallocate outputs.
-        numfreq = len(omega)
-        mag = empty((self.outputs, self.inputs, numfreq))
-        phase = empty((self.outputs, self.inputs, numfreq))
+        fresp : (p, m, len(s)) complex ndarray or (len(s),) complex ndarray 
+            The frequency response of the system. Array is ``(len(s), )`` if
+            and only if system is SISO and ``squeeze=True``.
 
-        omega.sort()
 
-        for k, w in enumerate(omega):
-            fresp = self._evalfr(w)
-            mag[:, :, k] = abs(fresp)
-            phase[:, :, k] = angle(fresp)
+        Raises
+        ------
+        ValueError
+            If `s` is not purely imaginary, because
+            :class:`FrequencyDomainData` systems are only defined at imaginary
+            frequency values.
+        """
+        if any(abs(np.array(s, ndmin=1).real) > 0):
+            raise ValueError("__call__: FRD systems can only accept "
+                            "purely imaginary frequencies")
+        # need to preserve array or scalar status
+        if hasattr(s, '__len__'):
+            return self.eval(np.asarray(s).imag, squeeze=squeeze)
+        else:
+            return self.eval(complex(s).imag, squeeze=squeeze)
 
-        return mag, phase, omega
+    def freqresp(self, omega):
+        """(deprecated) Evaluate transfer function at complex frequencies.
+
+        .. deprecated::0.9.0
+            Method has been given the more pythonic name
+            :meth:`FrequencyResponseData.frequency_response`. Or use
+            :func:`freqresp` in the MATLAB compatibility module.
+        """
+        warn("FrequencyResponseData.freqresp(omega) will be removed in a "
+             "future release of python-control; use "
+             "FrequencyResponseData.frequency_response(omega), or "
+             "freqresp(sys, omega) in the MATLAB compatibility module "
+             "instead", DeprecationWarning)
+        return self.frequency_response(omega)
 
     def feedback(self, other=1, sign=-1):
         """Feedback interconnection between two FRD objects."""
@@ -515,11 +515,10 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
             "Frequency ranges of FRD do not match, conversion not implemented")
 
     elif isinstance(sys, LTI):
-        omega.sort()
-        fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
-        for k, w in enumerate(omega):
-            fresp[:, :, k] = sys._evalfr(w)
-
+        omega = np.sort(omega)
+        fresp = sys(1j * omega)
+        if len(fresp.shape) == 1:
+            fresp = fresp[np.newaxis, np.newaxis, :]
         return FRD(fresp, omega, smooth=True)
 
     elif isinstance(sys, (int, float, complex, np.number)):

diff --git a/control/freqplot.py b/control/freqplot.py
@@ -151,14 +151,14 @@ def bode_plot(syslist, omega=None,
 
     Notes
     -----
-    1. Alternatively, you may use the lower-level method
-       ``(mag, phase, freq) = sys.freqresp(freq)`` to generate the frequency
-       response for a system, but it returns a MIMO response.
+    1. Alternatively, you may use the lower-level methods
+       :meth:`LTI.frequency_response` or ``sys(s)`` or ``sys(z)`` or to
+       generate the frequency response for a single system.
 
     2. If a discrete time model is given, the frequency response is plotted
-       along the upper branch of the unit circle, using the mapping z = exp(j
-       \\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
-       timebase.  If not timebase is specified (dt = True), dt is set to 1.
+       along the upper branch of the unit circle, using the mapping ``z = exp(1j
+       * omega * dt)`` where `omega` ranges from 0 to `pi/dt` and `dt` is the discrete
+       timebase.  If timebase not specified (``dt=True``), `dt` is set to 1.
 
     Examples
     --------
@@ -228,7 +228,7 @@ def bode_plot(syslist, omega=None,
                 nyquistfrq = None
 
             # Get the magnitude and phase of the system
-            mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega_sys)
+            mag_tmp, phase_tmp, omega_sys = sys.frequency_response(omega_sys)
             mag = np.atleast_1d(np.squeeze(mag_tmp))
             phase = np.atleast_1d(np.squeeze(phase_tmp))
 
@@ -588,7 +588,7 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
                 "Nyquist is currently only implemented for SISO systems.")
         else:
             # Get the magnitude and phase of the system
-            mag_tmp, phase_tmp, omega = sys.freqresp(omega)
+            mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
             mag = np.squeeze(mag_tmp)
             phase = np.squeeze(phase_tmp)
 
@@ -718,7 +718,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     omega_plot = omega / (2. * math.pi) if Hz else omega
 
     # TODO: Need to add in the mag = 1 lines
-    mag_tmp, phase_tmp, omega = S.freqresp(omega)
+    mag_tmp, phase_tmp, omega = S.frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['s'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -728,7 +728,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['s'].tick_params(labelbottom=False)
     plot_axes['s'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = (P * S).freqresp(omega)
+    mag_tmp, phase_tmp, omega = (P * S).frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['ps'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -738,7 +738,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['ps'].set_ylabel("$|PS|$" + " (dB)" if dB else "")
     plot_axes['ps'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = (C * S).freqresp(omega)
+    mag_tmp, phase_tmp, omega = (C * S).frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['cs'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -749,7 +749,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
     plot_axes['cs'].set_ylabel("$|CS|$" + " (dB)" if dB else "")
     plot_axes['cs'].grid(grid, which='both')
 
-    mag_tmp, phase_tmp, omega = T.freqresp(omega)
+    mag_tmp, phase_tmp, omega = T.frequency_response(omega)
     mag = np.squeeze(mag_tmp)
     if dB:
         plot_axes['t'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)