"""FrequencyResponseData(frdata, omega[, smooth])

#

# Class attributes

#

# These attributes are defined as class attributes so that they are

# documented properly. They are "overwritten" in __init__.

#

#: Number of system inputs.

#:

#: :meta hide-value:

ninputs = 1

#: Number of system outputs.

#:

#: :meta hide-value:

noutputs = 1

#: Squeeze processing parameter.

#:

#: By default, if a system is single-input, single-output (SISO) then

#: the outputs (and inputs) are returned as a 1D array (indexed by

#: frequency) and if a system is multi-input or multi-output, then the

#: outputs are returned as a 2D array (indexed by output and frequency)

#: or a 3D array (indexed by output, trace, and frequency). If

#: `squeeze` = True, access to the output response will remove

#: single-dimensional entries from the shape of the inputs and outputs

#: even if the system is not SISO. If `squeeze` = False, the output is

#: returned as a 3D array (indexed by the output, input, and frequency)

#: even if the system is SISO. The default value can be set using

#: config.defaults['control.squeeze_frequency_response'].

#:

#: :meta hide-value:

squeeze = None

_epsw = 1e-8 #: Bound for exact frequency match

def __init__(self, *args, **kwargs):

"""FrequencyResponseData(response, omega[, dt])

smooth = kwargs.pop('smooth', False)

#

# Process positional arguments

#

if len(args) == 3:

# Discrete time transfer function

dt = args[-1]

if 'dt' in kwargs:

warn("received multiple dt arguments, "

"using positional arg dt = %s" % dt)

kwargs['dt'] = dt

args = args[:-1]

if len(args) == 2:

if not isinstance(args[0], FRD) and isinstance(args[0], LTI):

# not an FRD, but still an LTI system, second argument

# should be the frequency range

otherlti = args[0]

self.omega = sort(np.asarray(args[1], dtype=float))

# calculate frequency response at specified points

if otherlti.isctime():

s = 1j * self.omega

self.frdata = otherlti(s, squeeze=False)

else:

z = np.exp(1j * self.omega * otherlti.dt)

self.frdata = otherlti(z, squeeze=False)

arg_dt = otherlti.dt

# Copy over signal and system names, if not specified

kwargs['inputs'] = kwargs.get('inputs', otherlti.input_labels)

kwargs['outputs'] = kwargs.get(

'outputs', otherlti.output_labels)

if not otherlti._generic_name_check():

kwargs['name'] = kwargs.get('name', _extended_system_name(

otherlti.name, prefix_suffix_name='sampled'))

else:

# The user provided a response and a freq vector

self.frdata = array(args[0], dtype=complex, ndmin=1)

if self.frdata.ndim == 1:

self.frdata = self.frdata.reshape(1, 1, -1)

self.omega = array(args[1], dtype=float, ndmin=1)

if self.frdata.ndim != 3 or self.omega.ndim != 1 or \

self.frdata.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"

" size")

arg_dt = None

elif len(args) == 1:

# Use the copy constructor.

if not isinstance(args[0], FRD):

raise TypeError(

"The one-argument constructor can only take in"

" an FRD object. Received %s." % type(args[0]))

self.omega = args[0].omega

self.frdata = args[0].frdata

arg_dt = args[0].dt

# Copy over signal and system names, if not specified

kwargs['inputs'] = kwargs.get('inputs', args[0].input_labels)

kwargs['outputs'] = kwargs.get('outputs', args[0].output_labels)

else:

raise ValueError(

"Needs 1 or 2 arguments; received %i." % len(args))

#

# Process keyword arguments

#

# If data was generated by a system, keep track of that (used when

# plotting data). Otherwise, use the system name, if given.

self.sysname = kwargs.pop('sysname', kwargs.get('name', None))

# Keep track of default properties for plotting

self.plot_phase = kwargs.pop('plot_phase', None)

self.title = kwargs.pop('title', None)

self.plot_type = kwargs.pop('plot_type', 'bode')

# 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._return_singvals=kwargs.pop('_return_singvals', False)

# Determine whether to squeeze the output

self.squeeze=kwargs.pop('squeeze', None)

if self.squeeze not in (None, True, False):

raise ValueError("unknown squeeze value")

defaults = {

'inputs': self.frdata.shape[1] if not getattr(

self, 'input_index', None) else self.input_labels,

'outputs': self.frdata.shape[0] if not getattr(

self, 'output_index', None) else self.output_labels,

'name': getattr(self, 'name', None)}

if arg_dt is not None:

if isinstance(args[0], LTI):

arg_dt = common_timebase(args[0].dt, arg_dt)

kwargs['dt'] = arg_dt

# Process signal names

name, inputs, outputs, states, dt = _process_iosys_keywords(

kwargs, defaults)

InputOutputSystem.__init__(

self, name=name, inputs=inputs, outputs=outputs, dt=dt, **kwargs)

# create interpolation functions

if smooth:

# Set the order of the fit

if self.omega.size < 2:

raise ValueError("can't smooth with only 1 frequency")

degree = 3 if self.omega.size > 3 else self.omega.size - 1

self._ifunc = empty((self.frdata.shape[0], self.frdata.shape[1]),

dtype=tuple)

for i in range(self.frdata.shape[0]):

for j in range(self.frdata.shape[1]):

self._ifunc[i, j], u = splprep(

u=self.omega, x=[real(self.frdata[i, j, :]),

imag(self.frdata[i, j, :])],

w=1.0/(absolute(self.frdata[i, j, :]) + 0.001),

s=0.0, k=degree)

else:

self._ifunc = None

#

# Frequency response properties

#

# Different properties of the frequency response that can be used for

# analysis and characterization.

#

@property

def magnitude(self):

"""Magnitude of the frequency response.

frdata = _process_frequency_response(

self, self.omega, self.frdata, squeeze=self.squeeze)

return NamedSignal(

np.abs(frdata), self.output_labels, self.input_labels)

@property

def phase(self):

"""Phase of the frequency response.

frdata = _process_frequency_response(

self, self.omega, self.frdata, squeeze=self.squeeze)

return NamedSignal(

np.angle(frdata), self.output_labels, self.input_labels)

@property

def frequency(self):

"""Frequencies at which the response is evaluated.

return self.omega

@property

def complex(self):

"""Complex value of the frequency response.

frdata = _process_frequency_response(

self, self.omega, self.frdata, squeeze=self.squeeze)

return NamedSignal(

frdata, self.output_labels, self.input_labels)

@property

def response(self):

warn("response property is deprecated; use complex", FutureWarning)

return self.complex

@property

def fresp(self):

warn("fresp attribute is deprecated; use frdata", FutureWarning)

return self.frdata

def __str__(self):

"""String representation of the transfer function."""

mimo = self.ninputs > 1 or self.noutputs > 1

outstr = [f"{InputOutputSystem.__str__(self)}"]

nl = "\n " if mimo else "\n"

sp = " " if mimo else ""

for i in range(self.ninputs):

for j in range(self.noutputs):

if mimo:

outstr.append(

"\nInput %i to output %i:" % (i + 1, j + 1))

outstr.append(nl + 'Freq [rad/s] Response')

outstr.append(sp + '------------ ---------------------')

outstr.extend(

[sp + '%12.3f %10.4g%+10.4gj' % (w, re, im)

for w, re, im in zip(self.omega,

real(self.frdata[j, i, :]),

imag(self.frdata[j, i, :]))])

return '\n'.join(outstr)

def _repr_eval_(self):

# Loadable format

out = "FrequencyResponseData(\n{d},\n{w}{smooth}".format(

d=repr(self.frdata), w=repr(self.omega),

smooth=(self._ifunc and ", smooth=True") or "")

out += self._dt_repr()

if len(labels := self._label_repr()) > 0:

out += ",\n" + labels

out += ")"

return out

def __neg__(self):

"""Negate a transfer function."""

return FRD(-self.frdata, self.omega)

def __add__(self, other):

"""Add two LTI objects (parallel connection)."""

if isinstance(other, FRD):

# verify that the frequencies match

if len(other.omega) != len(self.omega) or \

(other.omega != self.omega).any():

warn("Frequency points do not match; expect "

"truncation and interpolation.")

# Convert the second argument to a frequency response function.

# or re-base the frd to the current omega (if needed)

if isinstance(other, (int, float, complex, np.number)):

other = _convert_to_frd(

other, omega=self.omega,

inputs=self.ninputs, outputs=self.noutputs)

else:

other = _convert_to_frd(other, omega=self.omega)

# Promote SISO object to compatible dimension

if self.issiso() and not other.issiso():

self = np.ones((other.noutputs, other.ninputs)) * self

elif not self.issiso() and other.issiso():

other = np.ones((self.noutputs, self.ninputs)) * other

# Check that the input-output sizes are consistent.

if self.ninputs != other.ninputs:

raise ValueError(

"The first summand has %i input(s), but the " \

"second has %i." % (self.ninputs, other.ninputs))

if self.noutputs != other.noutputs:

raise ValueError(

"The first summand has %i output(s), but the " \

"second has %i." % (self.noutputs, other.noutputs))

return FRD(self.frdata + other.frdata, other.omega)

def __radd__(self, other):

"""Right add two LTI objects (parallel connection)."""

return self + other

def __sub__(self, other):

"""Subtract two LTI objects."""

return self + (-other)

def __rsub__(self, other):

"""Right subtract two LTI objects."""

return other + (-self)

def __mul__(self, other):

"""Multiply two LTI objects (serial connection)."""

# Convert the second argument to a transfer function.

if isinstance(other, (int, float, complex, np.number)):

return FRD(self.frdata * other, self.omega,

smooth=(self._ifunc is not None))

else:

other = _convert_to_frd(other, omega=self.omega)

# Promote SISO object to compatible dimension

if self.issiso() and not other.issiso():

self = bdalg.append(*([self] * other.noutputs))

elif not self.issiso() and other.issiso():

other = bdalg.append(*([other] * self.ninputs))

# Check that the input-output sizes are consistent.

if self.ninputs != other.noutputs:

raise ValueError(

"H = G1*G2: input-output size mismatch: "

"G1 has %i input(s), G2 has %i output(s)." %

(self.ninputs, other.noutputs))

inputs = other.ninputs

outputs = self.noutputs

frdata = empty((outputs, inputs, len(self.omega)),

dtype=self.frdata.dtype)

for i in range(len(self.omega)):

frdata[:, :, i] = self.frdata[:, :, i] @ other.frdata[:, :, i]

return FRD(frdata, self.omega,

smooth=(self._ifunc is not None) and

(other._ifunc is not None))

def __rmul__(self, other):

"""Right Multiply two LTI objects (serial connection)."""

# Convert the second argument to an frd function.

if isinstance(other, (int, float, complex, np.number)):

return FRD(self.frdata * other, self.omega,

smooth=(self._ifunc is not None))

else:

other = _convert_to_frd(other, omega=self.omega)

# Promote SISO object to compatible dimension

if self.issiso() and not other.issiso():

self = bdalg.append(*([self] * other.ninputs))

elif not self.issiso() and other.issiso():

other = bdalg.append(*([other] * self.noutputs))

# Check that the input-output sizes are consistent.

if self.noutputs != other.ninputs:

raise ValueError(

"H = G1*G2: input-output size mismatch: "

"G1 has %i input(s), G2 has %i output(s)." %

(other.ninputs, self.noutputs))

inputs = self.ninputs

outputs = other.noutputs

frdata = empty((outputs, inputs, len(self.omega)),

dtype=self.frdata.dtype)

for i in range(len(self.omega)):

frdata[:, :, i] = other.frdata[:, :, i] @ self.frdata[:, :, i]

return FRD(frdata, self.omega,

smooth=(self._ifunc is not None) and

(other._ifunc is not None))

# TODO: Division of MIMO transfer function objects is not written yet.

def __truediv__(self, other):

"""Divide two LTI objects."""

if isinstance(other, (int, float, complex, np.number)):

return FRD(self.frdata * (1/other), self.omega,

smooth=(self._ifunc is not None))

else:

other = _convert_to_frd(other, omega=self.omega)

if (other.ninputs > 1 or other.noutputs > 1):

# FRD.__truediv__ is currently only implemented for SISO systems

return NotImplemented

return FRD(self.frdata/other.frdata, self.omega,

smooth=(self._ifunc is not None) and

(other._ifunc is not None))

# TODO: Division of MIMO transfer function objects is not written yet.

def __rtruediv__(self, other):

"""Right divide two LTI objects."""

if isinstance(other, (int, float, complex, np.number)):

return FRD(other / self.frdata, self.omega,

smooth=(self._ifunc is not None))

else:

other = _convert_to_frd(other, omega=self.omega)

if (self.ninputs > 1 or self.noutputs > 1):

# FRD.__rtruediv__ is currently only implemented for SISO systems

return NotImplemented

return other / self

def __pow__(self, other):

if not type(other) == int:

raise ValueError("Exponent must be an integer")

if other == 0:

return FRD(ones(self.frdata.shape), self.omega,

smooth=(self._ifunc is not None)) # unity

if other > 0:

return self * (self**(other-1))

if other < 0:

return (FRD(ones(self.frdata.shape), self.omega) / self) * \

(self**(other+1))

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

# 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=None):

"""Evaluate a transfer function at a frequency point.

omega_array = np.array(omega, ndmin=1) # array of frequencies

# 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("eval can only accept real-valued frequencies")

if self._ifunc is None:

elements = np.isin(self.omega, omega) # binary array

if sum(elements) < len(omega_array):

raise ValueError(

"not all frequencies are in frequency list of FRD "

"system. Try an interpolating FRD for additional points.")

else:

out = self.frdata[:, :, elements]

else:

out = empty((self.noutputs, self.ninputs, len(omega_array)),

dtype=complex)

for i in range(self.noutputs):

for j in range(self.ninputs):

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]

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

def __call__(self, x=None, squeeze=None, return_magphase=None):

"""Evaluate system transfer function at point in complex plane.

if x 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

if return_magphase is not None:

raise ValueError("return_magphase not allowed when x != None")

# Make sure that we are operating on a simple list

if len(np.atleast_1d(x).shape) > 1:

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

if any(abs(np.atleast_1d(x).real) > 0):

raise ValueError("__call__: FRD systems can only accept "

"purely imaginary frequencies")

# need to preserve array or scalar status

if hasattr(x, '__len__'):

return self.eval(np.asarray(x).imag, squeeze=squeeze)

else:

return self.eval(complex(x).imag, squeeze=squeeze)

# Implement iter to allow assigning to a tuple

def __iter__(self):

frdata = _process_frequency_response(

self, self.omega, self.frdata, squeeze=self.squeeze)

if self._return_singvals:

# Legacy processing for singular values

return iter((self.frdata[:, 0, :], self.omega))

elif not self.return_magphase:

return iter((self.omega, frdata))

return iter((np.abs(frdata), np.angle(frdata), self.omega))

def __getitem__(self, key):

if not isinstance(key, Iterable) or len(key) != 2:

# Implement (thin) getitem to allow access via legacy indexing

return list(self.__iter__())[key]

# Convert signal names to integer offsets (via NamedSignal object)

iomap = NamedSignal(

self.frdata[:, :, 0], self.output_labels, self.input_labels)

indices = iomap._parse_key(key, level=1) # ignore index checks

outdx, outputs = _process_subsys_index(indices[0], self.output_labels)

inpdx, inputs = _process_subsys_index(indices[1], self.input_labels)

# Create the system name

sysname = config.defaults['iosys.indexed_system_name_prefix'] + \

self.name + config.defaults['iosys.indexed_system_name_suffix']

return FrequencyResponseData(

self.frdata[outdx, :][:, inpdx], self.omega, self.dt,

inputs=inputs, outputs=outputs, name=sysname)

# 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.

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", FutureWarning)

return self.frequency_response(omega)

def feedback(self, other=1, sign=-1):

"""Feedback interconnection between two FRD objects.

other = _convert_to_frd(other, omega=self.omega)

if (self.noutputs != other.ninputs or self.ninputs != other.noutputs):

raise ValueError(

"FRD.feedback, inputs/outputs mismatch")

# TODO: handle omega re-mapping

# reorder array axes in order to leverage numpy broadcasting

myfrdata = np.moveaxis(self.frdata, 2, 0)

otherfrdata = np.moveaxis(other.frdata, 2, 0)

I_AB = eye(self.ninputs)[np.newaxis, :, :] + otherfrdata @ myfrdata

resfrdata = (myfrdata @ linalg.inv(I_AB))

frdata = np.moveaxis(resfrdata, 0, 2)

return FRD(frdata, other.omega, smooth=(self._ifunc is not None))

def append(self, other):

"""Append a second model to the present model.

other = _convert_to_frd(other, omega=self.omega, inputs=other.ninputs,

outputs=other.noutputs)

# TODO: handle omega re-mapping

new_frdata = np.zeros(

(self.noutputs + other.noutputs, self.ninputs + other.ninputs,

self.omega.shape[-1]), dtype=complex)

new_frdata[:self.noutputs, :self.ninputs, :] = np.reshape(

self.frdata, (self.noutputs, self.ninputs, -1))

new_frdata[self.noutputs:, self.ninputs:, :] = np.reshape(

other.frdata, (other.noutputs, other.ninputs, -1))

return FRD(new_frdata, self.omega, smooth=(self._ifunc is not None))

# Plotting interface

def plot(self, plot_type=None, *args, **kwargs):

"""Plot the frequency response using Bode or singular values plot.

from .freqplot import bode_plot, singular_values_plot

from .nichols import nichols_plot

if plot_type is None:

plot_type = self.plot_type

if plot_type == 'bode':

return bode_plot(self, *args, **kwargs)

elif plot_type == 'nichols':

return nichols_plot(self, *args, **kwargs)

elif plot_type == 'svplot':

return singular_values_plot(self, *args, **kwargs)

else:

raise ValueError(f"unknown plot type '{plot_type}'")

# Convert to pandas

def to_pandas(self):

"""Convert response data to pandas data frame.

if not pandas_check():

ImportError('pandas not installed')

import pandas

# Create a dict for setting up the data frame

data = {'omega': self.omega}

data.update(

{'H_{%s, %s}' % (out, inp): self.frdata[i, j] \

for i, out in enumerate(self.output_labels) \

for j, inp in enumerate(self.input_labels)})

"""Convert a system to frequency response data form (if needed).

if isinstance(sys, FRD):

omega.sort()

if len(omega) == len(sys.omega) and \

(abs(omega - sys.omega) < FRD._epsw).all():

# frequencies match, and system was already frd; simply use

return sys

raise NotImplementedError(

"Frequency ranges of FRD do not match, conversion not implemented")

elif isinstance(sys, LTI):

omega = np.sort(omega)