frequency_response() returns FRD; FRD allows return_magphase

murrayrm · murrayrm · commit 6c3c630dcd6f · 2022-04-01T21:19:49.000-07:00
diff --git a/control/frdata.py b/control/frdata.py
@@ -44,11 +44,14 @@
 """
 
 # External function declarations
+from copy import copy
 from warnings import warn
+
 import numpy as np
 from numpy import angle, array, empty, ones, \
     real, imag, absolute, eye, linalg, where, sort
 from scipy.interpolate import splprep, splev
+
 from .lti import LTI, _process_frequency_response
 from .exception import pandas_check
 from .namedio import _NamedIOSystem
@@ -141,7 +144,7 @@ def __init__(self, *args, **kwargs):
         The default constructor is FRD(d, w), where w is an iterable of
         frequency points, and d is the matching frequency data.
 
-        If d is a single list, 1d array, or tuple, a SISO system description
+        If d is a single list, 1D array, or tuple, a SISO system description
         is assumed. d can also be
 
         To call the copy constructor, call FRD(sys), where sys is a
@@ -170,13 +173,12 @@ def __init__(self, *args, **kwargs):
 
             else:
                 # The user provided a response and a freq vector
-                self.fresp = array(args[0], dtype=complex)
-                if len(self.fresp.shape) == 1:
-                    self.fresp = self.fresp.reshape(1, 1, len(args[0]))
-                self.omega = array(args[1], dtype=float)
-                if len(self.fresp.shape) != 3 or \
-                        self.fresp.shape[-1] != self.omega.shape[-1] or \
-                        len(self.omega.shape) != 1:
+                self.fresp = array(args[0], dtype=complex, ndmin=1)
+                if self.fresp.ndim == 1:
+                    self.fresp = self.fresp.reshape(1, 1, -1)
+                self.omega = array(args[1], dtype=float, ndmin=1)
+                if self.fresp.ndim != 3 or self.omega.ndim != 1 or \
+                        self.fresp.shape[-1] != self.omega.shape[-1]:
                     raise TypeError(
                         "The frequency data constructor needs a 1-d or 3-d"
                         " response data array and a matching frequency vector"
@@ -206,6 +208,12 @@ def __init__(self, *args, **kwargs):
 
         # Keep track of return type
         self.return_magphase=kwargs.pop('return_magphase', False)
+        if self.return_magphase not in (True, False):
+            raise ValueError("unknown return_magphase value")
+
+        self.squeeze=kwargs.pop('squeeze', None)
+        if self.squeeze not in (None, True, False):
+            raise ValueError("unknown squeeze value")
 
         # Make sure there were no extraneous keywords
         if kwargs:
@@ -477,7 +485,7 @@ def eval(self, omega, squeeze=None):
 
         return _process_frequency_response(self, omega, out, squeeze=squeeze)
 
-    def __call__(self, s, squeeze=None):
+    def __call__(self, s=None, squeeze=None, return_magphase=None):
         """Evaluate system's transfer function at complex frequencies.
 
         Returns the complex frequency response `sys(s)` of system `sys` with
@@ -490,17 +498,31 @@ def __call__(self, s, squeeze=None):
         For a frequency response data object, the argument must be an
         imaginary number (since only the frequency response is defined).
 
+        If ``s`` is not given, this function creates a copy of a frequency
+        response data object with a different set of output settings.
+
         Parameters
         ----------
         s : complex scalar or 1D array_like
-            Complex frequencies
-        squeeze : bool, optional (default=True)
+            Complex frequencies.  If not specified, return a copy of the
+            frequency response data object with updated settings for output
+            processing (``squeeze``, ``return_magphase``).
+
+        squeeze : bool, optional
             If squeeze=True, remove single-dimensional entries from the shape
             of the output even if the system is not SISO. If squeeze=False,
             keep all indices (output, input and, if omega is array_like,
             frequency) even if the system is SISO. The default value can be
             set using config.defaults['control.squeeze_frequency_response'].
 
+        return_magphase : bool, optional
+            If True, then a frequency response data object will enumerate as a
+            tuple of the form (mag, phase, omega) where where ``mag`` is the
+            magnitude (absolute value, not dB or log10) of the system
+            frequency response, ``phase`` is the wrapped phase in radians of
+            the system frequency response, and ``omega`` is the (sorted)
+            frequencies at which the response was evaluated.
+
         Returns
         -------
         fresp : complex ndarray
@@ -519,6 +541,17 @@ def __call__(self, s, squeeze=None):
             frequency values.
 
         """
+        if s is None:
+            # Create a copy of the response with new keywords
+            response = copy(self)
+
+            # Update any keywords that we were passed
+            response.squeeze = self.squeeze if squeeze is None else squeeze
+            response.return_magphase = self.return_magphase \
+                if return_magphase is None else return_magphase
+
+            return response
+
         # Make sure that we are operating on a simple list
         if len(np.atleast_1d(s).shape) > 1:
             raise ValueError("input list must be 1D")
@@ -533,6 +566,22 @@ def __call__(self, s, squeeze=None):
         else:
             return self.eval(complex(s).imag, squeeze=squeeze)
 
+    # Implement iter to allow assigning to a tuple
+    def __iter__(self):
+        fresp = _process_frequency_response(
+            self, self.omega, self.fresp, squeeze=self.squeeze)
+        if not self.return_magphase:
+            return iter((self.omega, fresp))
+        return iter((np.abs(fresp), np.angle(fresp), self.omega))
+
+    # Implement (thin) getitem to allow access via legacy indexing
+    def __getitem__(self, index):
+        return list(self.__iter__())[index]
+
+    # Implement (thin) len to emulate legacy testing interface
+    def __len__(self):
+        return 3 if self.return_magphase else 2
+
     def freqresp(self, omega):
         """(deprecated) Evaluate transfer function at complex frequencies.
 
diff --git a/control/lti.py b/control/lti.py
@@ -77,9 +77,9 @@ def _set_inputs(self, value):
         """
         Deprecated attribute; use :attr:`ninputs` instead.
 
-        The ``input`` attribute was used to store the number of system inputs.
-        It is no longer used.  If you need access to the number of inputs for
-        an LTI system, use :attr:`ninputs`.
+        The ``inputs`` attribute was used to store the number of system
+        inputs.  It is no longer used.  If you need access to the number
+        of inputs for an LTI system, use :attr:`ninputs`.
         """)
 
     def _get_outputs(self):
@@ -100,7 +100,7 @@ def _set_outputs(self, value):
         """
         Deprecated attribute; use :attr:`noutputs` instead.
 
-        The ``output`` attribute was used to store the number of system
+        The ``outputs`` attribute was used to store the number of system
         outputs.  It is no longer used.  If you need access to the number of
         outputs for an LTI system, use :attr:`noutputs`.
         """)
@@ -197,17 +197,21 @@ def frequency_response(self, omega, squeeze=None):
 
         Returns
         -------
-        mag : ndarray
-            The magnitude (absolute value, not dB or log10) of the system
-            frequency response.  If the system is SISO and squeeze is not
-            True, the array is 1D, indexed by frequency.  If the system is not
-            SISO or squeeze is False, the array is 3D, indexed by the output,
-            input, and frequency.  If ``squeeze`` is True then
-            single-dimensional axes are removed.
-        phase : ndarray
-            The wrapped phase in radians of the system frequency response.
-        omega : ndarray
-            The (sorted) frequencies at which the response was evaluated.
+        response : :class:`FrequencyReponseData`
+            Frequency response data object representing the frequency
+            response.  This object can be assigned to a tuple using
+
+                mag, phase, omega = response
+
+            where ``mag`` is the magnitude (absolute value, not dB or log10)
+            of the system frequency response, ``phase`` is the wrapped phase
+            in radians of the system frequency response, and ``omega`` is the
+            (sorted) frequencies at which the response was evaluated.  If the
+            system is SISO and squeeze is not True, ``mag`` and ``phase`` are
+            1D, indexed by frequency.  If the system is not SISO or squeeze is
+            False, the array is 3D, indexed by the output, input, and
+            frequency.  If ``squeeze`` is True then single-dimensional axes
+            are removed.
 
         """
         omega = np.sort(np.array(omega, ndmin=1))
@@ -218,8 +222,12 @@ def frequency_response(self, omega, squeeze=None):
             s = np.exp(1j * omega * self.dt)
         else:
             s = 1j * omega
-        response = self.__call__(s, squeeze=squeeze)
-        return abs(response), angle(response), omega
+
+        # Return the data as a frequency response data object
+        from .frdata import FrequencyResponseData
+        response = self.__call__(s)
+        return FrequencyResponseData(
+            response, omega, return_magphase=True, squeeze=squeeze)
 
     def dcgain(self):
         """Return the zero-frequency gain"""
diff --git a/control/tests/lti_test.py b/control/tests/lti_test.py
@@ -267,7 +267,7 @@ def test_squeeze_exceptions(self, fcn):
             sys = fcn(ct.rss(2, 1, 1))
 
         with pytest.raises(ValueError, match="unknown squeeze value"):
-            sys.frequency_response([1], squeeze=1)
+            resp = sys.frequency_response([1], squeeze='siso')
         with pytest.raises(ValueError, match="unknown squeeze value"):
             sys([1j], squeeze='siso')
         with pytest.raises(ValueError, match="unknown squeeze value"):
diff --git a/control/timeresp.py b/control/timeresp.py
@@ -463,7 +463,7 @@ def __call__(self, **kwargs):
         # Update any keywords that we were passed
         response.transpose = kwargs.pop('transpose', self.transpose)
         response.squeeze = kwargs.pop('squeeze', self.squeeze)
-        response.return_x = kwargs.pop('return_x', self.squeeze)
+        response.return_x = kwargs.pop('return_x', self.return_x)
 
         # Check for new labels
         input_labels = kwargs.pop('input_labels', None)
diff --git a/control/xferfcn.py b/control/xferfcn.py
@@ -61,6 +61,7 @@
 from re import sub
 from .lti import LTI, common_timebase, isdtime, _process_frequency_response
 from .exception import ControlMIMONotImplemented
+from .frdata import FrequencyResponseData
 from .namedio import _NamedIOSystem, _process_signal_list
 from . import config
 
@@ -1382,6 +1383,9 @@ def _convert_to_transfer_function(sys, inputs=1, outputs=1):
 
         return TransferFunction(num, den)
 
+    elif isinstance(sys, FrequencyResponseData):
+        raise TypeError("Can't convert given FRD to TransferFunction system.")
+
     # If this is array-like, try to create a constant feedthrough
     try:
         D = array(sys, ndmin=2)
