Skip to content

Standardize squeeze processing in frequency response functions #507

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
murrayrm merged 5 commits into from
Jan 19, 2021
Merged

Standardize squeeze processing in frequency response functions #507

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
848112d
unify use of squeeze keyword in frequency response functions
murrayrm Jan 13, 2021
f367cf6
unify frequency response processing + unit tests
murrayrm Jan 17, 2021
451d6d2
change control.squeeze to control.squeeze_frequency_response
murrayrm Jan 17, 2021
90da4fb
update squeeze processing for freq response + unit tests, doc updates
murrayrm Jan 17, 2021
b338e32
address @sawyerbfuller review comments
murrayrm Jan 18, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion control/config.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,8 @@

# Package level default values
_control_defaults = {
'control.default_dt':0
'control.default_dt': 0,
'control.squeeze_frequency_response': None
}
defaults = dict(_control_defaults)

Expand Down
72 changes: 47 additions & 25 deletions control/frdata.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,8 @@
from numpy import angle, array, empty, ones, \
real, imag, absolute, eye, linalg, where, dot, sort
from scipy.interpolate import splprep, splev
from .lti import LTI
from .lti import LTI, _process_frequency_response
from . import config

__all__ = ['FrequencyResponseData', 'FRD', 'frd']

Expand Down Expand Up @@ -343,7 +344,7 @@ def __pow__(self, other):
# G(s) for a transfer function and G(omega) for an FRD object.
# 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):
def eval(self, omega, squeeze=None):
"""Evaluate a transfer function at angular frequency omega.

Note that a "normal" FRD only returns values for which there is an
Expand All @@ -352,19 +353,33 @@ def eval(self, omega, squeeze=True):

Parameters
----------
omega : float or array_like
omega : float or 1D 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.
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'].

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``.
fresp : complex ndarray
The frequency response of the system. If the system is SISO and
squeeze is not True, the shape of the array matches the shape of
omega. If the system is not SISO or squeeze is False, the first
two dimensions of the array are indices for the output and input
and the remaining dimensions match omega. If ``squeeze`` is True
then single-dimensional axes are removed.

"""
omega_array = np.array(omega, ndmin=1) # array-like version of omega

# Make sure that we are operating on a simple list
if len(omega_array.shape) > 1:
raise ValueError("input list must be 1D")

# Make sure that frequencies are all real-valued
if any(omega_array.imag > 0):
raise ValueError("FRD.eval can only accept real-valued omega")

Expand All @@ -384,16 +399,12 @@ def eval(self, omega, squeeze=True):
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

def __call__(self, s, squeeze=True):

return _process_frequency_response(self, omega, out, squeeze=squeeze)

def __call__(self, s, squeeze=None):
"""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.
Expand All @@ -403,18 +414,24 @@ def __call__(self, s, squeeze=True):

Parameters
----------
s : complex scalar or array_like
s : complex scalar or 1D 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.
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'].

Returns
-------
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``.

fresp : complex ndarray
The frequency response of the system. If the system is SISO and
squeeze is not True, the shape of the array matches the shape of
omega. If the system is not SISO or squeeze is False, the first
two dimensions of the array are indices for the output and input
and the remaining dimensions match omega. If ``squeeze`` is True
then single-dimensional axes are removed.

Raises
------
Expand All @@ -423,9 +440,14 @@ def __call__(self, s, squeeze=True):
:class:`FrequencyDomainData` systems are only defined at imaginary
frequency values.
"""
# Make sure that we are operating on a simple list
if len(np.array(s, ndmin=1).shape) > 1:

@sawyerbfuller sawyerbfuller Jan 18, 2021

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if len(np.array(s, ndmin=1).shape) > 1:
if len(np.atleast1d(s).shape) > 1:

raise ValueError("input list must be 1D")

if any(abs(np.array(s, ndmin=1).real) > 0):

@sawyerbfuller sawyerbfuller Jan 18, 2021 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if any(abs(np.array(s, ndmin=1).real) > 0):
if any(abs(np.atleast1d(s).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)
Expand Down
120 changes: 85 additions & 35 deletions control/lti.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
import numpy as np
from numpy import absolute, real, angle, abs
from warnings import warn
from . import config

__all__ = ['issiso', 'timebase', 'common_timebase', 'timebaseEqual',
'isdtime', 'isctime', 'pole', 'zero', 'damp', 'evalfr',
Expand Down Expand Up @@ -111,7 +112,7 @@ def damp(self):
Z = -real(splane_poles)/wn
return wn, Z, poles

def frequency_response(self, omega, squeeze=True):
def frequency_response(self, omega, squeeze=None):
"""Evaluate the linear time-invariant system at an array of angular
frequencies.

Expand All @@ -124,30 +125,36 @@ def frequency_response(self, omega, squeeze=True):

G(exp(j*omega*dt)) = mag*exp(j*phase).

In general the system may be multiple input, multiple output (MIMO), where
`m = self.inputs` number of inputs and `p = self.outputs` number of
outputs.
In general the system may be multiple input, multiple output (MIMO),
where `m = self.inputs` number of inputs and `p = self.outputs` number
of outputs.

Parameters
----------
omega : float or array_like
omega : float or 1D array_like
A list, tuple, array, or scalar value of frequencies in
radians/sec at which the system will be evaluated.
squeeze : bool, optional (default=True)
If True and the system is single input single output (SISO), i.e. `m=1`,
`p=1`, return a 1D array rather than a 3D array.
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'].

Returns
-------
mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray
mag : ndarray
The magnitude (absolute value, not dB or log10) of the system
frequency response. Array is ``(len(omega), )`` if
and only if system is SISO and ``squeeze=True``.
phase : (p, m, len(omega)) ndarray or (len(omega),) ndarray
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.

"""
omega = np.sort(np.array(omega, ndmin=1))
if isdtime(self, strict=True):
Expand Down Expand Up @@ -463,9 +470,8 @@ def damp(sys, doprint=True):
(p.real, p.imag, d, w))
return wn, damping, poles

def evalfr(sys, x, squeeze=True):
"""
Evaluate the transfer function of an LTI system for complex frequency x.
def evalfr(sys, x, squeeze=None):
"""Evaluate the transfer function of an LTI system for complex frequency x.

Returns the complex frequency response `sys(x)` where `x` is `s` for
continuous-time systems and `z` for discrete-time systems, with
Expand All @@ -481,17 +487,24 @@ def evalfr(sys, x, squeeze=True):
----------
sys: StateSpace or TransferFunction
Linear system
x : complex scalar or array_like
x : complex scalar or 1D array_like
Complex frequency(s)
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.
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'].

Returns
-------
fresp : (p, m, len(x)) complex ndarray 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``.
fresp : complex ndarray
The frequency response of the system. If the system is SISO and
squeeze is not True, the shape of the array matches the shape of
omega. If the system is not SISO or squeeze is False, the first two
dimensions of the array are indices for the output and input and the
remaining dimensions match omega. If ``squeeze`` is True then
single-dimensional axes are removed.

See Also
--------
Expand All @@ -511,13 +524,13 @@ def evalfr(sys, x, squeeze=True):
>>> # This is the transfer function matrix evaluated at s = i.

.. todo:: Add example with MIMO system

"""
return sys.__call__(x, squeeze=squeeze)

def freqresp(sys, omega, squeeze=True):
"""
Frequency response of an LTI system at multiple angular frequencies.

def freqresp(sys, omega, squeeze=None):
"""Frequency response of an LTI system at multiple angular frequencies.

In general the system may be multiple input, multiple output (MIMO), where
`m = sys.inputs` number of inputs and `p = sys.outputs` number of
outputs.
Expand All @@ -526,22 +539,27 @@ def freqresp(sys, omega, squeeze=True):
----------
sys: StateSpace or TransferFunction
Linear system
omega : float or array_like
omega : float or 1D 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.
squeeze : bool, optional (default=True)
If True and `sys` is single input, single output (SISO), returns
1D array rather than a 3D array.
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'].

Returns
-------
mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray
mag : ndarray
The magnitude (absolute value, not dB or log10) of the system
frequency response. Array is ``(len(omega), )`` if and only if system
is SISO and ``squeeze=True``.

phase : (p, m, len(omega)) ndarray or (len(omega),) ndarray
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

@sawyerbfuller sawyerbfuller Jan 18, 2021

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In ‘notes’ on line 575ish below, this function is actually a wrapper for :meth:’lti.frequency_response’

@murrayrm murrayrm Jan 18, 2021

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because of the way the documentation is built, the LTI class does not show up, so the documentation appears in the StateSpace and TransferFunction classes. Leaving unchanged.

The wrapped phase in radians of the system frequency response.
omega : ndarray
The list of sorted frequencies at which the response was
Expand Down Expand Up @@ -579,6 +597,7 @@ def freqresp(sys, omega, squeeze=True):
#>>> # input to the 1st output, and the phase (in radians) of the
#>>> # frequency response from the 1st input to the 2nd output, for
#>>> # s = 0.1i, i, 10i.

"""
return sys.frequency_response(omega, squeeze=squeeze)

Expand All @@ -593,3 +612,34 @@ def dcgain(sys):
at the origin
"""
return sys.dcgain()


# Process frequency responses in a uniform way
def _process_frequency_response(sys, omega, out, squeeze=None):
# Set value of squeeze argument if not set
if squeeze is None:
squeeze = config.defaults['control.squeeze_frequency_response']

if not hasattr(omega, '__len__'):
# received a scalar x, squeeze down the array along last dim
out = np.squeeze(out, axis=2)

#
# Get rid of unneeded dimensions
#
# There are three possible values for the squeeze keyword at this point:
#
# squeeze=None: squeeze input/output axes iff SISO
# squeeze=True: squeeze all single dimensional axes (ala numpy)
# squeeze-False: don't squeeze any axes
#
if squeeze is True:
# Squeeze everything that we can if that's what the user wants
return np.squeeze(out)
elif squeeze is None and sys.issiso():
# SISO system output squeezed unless explicitly specified otherwise
return out[0][0]
elif squeeze is False or squeeze is None:
return out
else:
raise ValueError("unknown squeeze value")
Loading