'freqplot.feature_periphery_decades': 1,

'freqplot.number_of_samples': 1000,

'freqplot.dB': False, # Plot gain in dB

'freqplot.deg': True, # Plot phase in degrees

'freqplot.Hz': False, # Plot frequency in Hertz

'freqplot.grid': True, # Turn on grid for gain and phase

'freqplot.wrap_phase': False, # Wrap the phase plot at a given value

'freqplot.freq_label': "Frequency [{units}]",

'freqplot.share_magnitude': 'row',

'freqplot.share_phase': 'row',

'freqplot.share_frequency': 'col',

'freqplot.title_frame': 'axes',

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

if plot_type == None:

for response in self:

if plot_type is not None and response.plot_type != plot_type:

raise TypeError(

"inconsistent plot_types in data; set plot_type "

"to 'bode', 'nichols', or 'svplot'")

plot_type = response.plot_type

# Use FRD plot method, which can handle lists via plot functions

return FrequencyResponseData.plot(

data, omega=None, *fmt, ax=None, omega_limits=None, omega_num=None,

plot=None, plot_magnitude=True, plot_phase=None,

overlay_outputs=None, overlay_inputs=None, phase_label=None,

magnitude_label=None, label=None, display_margins=None,

margins_method='best', title=None, sharex=None, sharey=None, **kwargs):

"""Bode plot for a system.

#

# Process keywords and set defaults

#

# Make a copy of the kwargs dictionary since we will modify it

kwargs = dict(kwargs)

# Get values for params (and pop from list to allow keyword use in plot)

dB = config._get_param(

'freqplot', 'dB', kwargs, _freqplot_defaults, pop=True)

deg = config._get_param(

'freqplot', 'deg', kwargs, _freqplot_defaults, pop=True)

Hz = config._get_param(

'freqplot', 'Hz', kwargs, _freqplot_defaults, pop=True)

grid = config._get_param(

'freqplot', 'grid', kwargs, _freqplot_defaults, pop=True)

wrap_phase = config._get_param(

'freqplot', 'wrap_phase', kwargs, _freqplot_defaults, pop=True)

initial_phase = config._get_param(

'freqplot', 'initial_phase', kwargs, None, pop=True)

rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True)

title_frame = config._get_param(

'freqplot', 'title_frame', kwargs, _freqplot_defaults, pop=True)

# Set the default labels

freq_label = config._get_param(

'freqplot', 'freq_label', kwargs, _freqplot_defaults, pop=True)

if magnitude_label is None:

magnitude_label = "Magnitude [dB]" if dB else "Magnitude"

if phase_label is None:

phase_label = "Phase [deg]" if deg else "Phase [rad]"

# Use sharex and sharey as proxies for share_{magnitude, phase, frequency}

if sharey is not None:

if 'share_magnitude' in kwargs or 'share_phase' in kwargs:

ValueError(

"sharey cannot be present with share_magnitude/share_phase")

kwargs['share_magnitude'] = sharey

kwargs['share_phase'] = sharey

if sharex is not None:

if 'share_frequency' in kwargs:

ValueError(

"sharex cannot be present with share_frequency")

kwargs['share_frequency'] = sharex

# Legacy keywords for margins

display_margins = config._process_legacy_keyword(

kwargs, 'margins', 'display_margins', display_margins)

if kwargs.pop('margin_info', False):

warnings.warn(

"keyword 'margin_info' is deprecated; "

"use 'display_margins='overlay'")

if display_margins is False:

raise ValueError(

"conflicting_keywords: `display_margins` and `margin_info`")

margins_method = config._process_legacy_keyword(

kwargs, 'method', 'margins_method', margins_method)

if not isinstance(data, (list, tuple)):

data = [data]

#

# Pre-process the data to be plotted (unwrap phase, limit frequencies)

#

# To maintain compatibility with legacy uses of bode_plot(), we do some

# initial processing on the data, specifically phase unwrapping and

# setting the initial value of the phase. If bode_plot is called with

# plot == False, then these values are returned to the user (instead of

# the list of lines created, which is the new output for _plot functions.

#

# If we were passed a list of systems, convert to data

if any([isinstance(

sys, (StateSpace, TransferFunction)) for sys in data]):

data = frequency_response(

data, omega=omega, omega_limits=omega_limits,

omega_num=omega_num, Hz=Hz)

else:

# Generate warnings if frequency keywords were given

if omega_num is not None:

warnings.warn("`omega_num` ignored when passed response data")

elif omega is not None:

warnings.warn("`omega` ignored when passed response data")

# Check to make sure omega_limits is sensible

if omega_limits is not None and \

(len(omega_limits) != 2 or omega_limits[1] <= omega_limits[0]):

raise ValueError(f"invalid limits: {omega_limits=}")

# If plot_phase is not specified, check the data first, otherwise true

if plot_phase is None:

plot_phase = True if data[0].plot_phase is None else data[0].plot_phase

if not plot_magnitude and not plot_phase:

raise ValueError(

"plot_magnitude and plot_phase both False; no data to plot")

mag_data, phase_data, omega_data = [], [], []

for response in data:

noutputs, ninputs = response.noutputs, response.ninputs

if initial_phase is None:

# Start phase in the range 0 to -360 w/ initial phase = 0

# TODO: change this to 0 to 270 (?)

# If wrap_phase is true, use 0 instead (phase \in (-pi, pi])

initial_phase_value = -math.pi if wrap_phase is not True else 0

elif isinstance(initial_phase, (int, float)):

# Allow the user to override the default calculation

if deg:

initial_phase_value = initial_phase/180. * math.pi

else:

initial_phase_value = initial_phase

else:

raise ValueError("initial_phase must be a number.")

# Reshape the phase to allow standard indexing

phase = response.phase.copy().reshape((noutputs, ninputs, -1))

# Shift and wrap the phase

for i, j in itertools.product(range(noutputs), range(ninputs)):

# Shift the phase if needed

if abs(phase[i, j, 0] - initial_phase_value) > math.pi:

phase[i, j] -= 2*math.pi * round(

(phase[i, j, 0] - initial_phase_value) / (2*math.pi))

# Phase wrapping

if wrap_phase is False:

phase[i, j] = unwrap(phase[i, j]) # unwrap the phase

elif wrap_phase is True:

pass # default calc OK

elif isinstance(wrap_phase, (int, float)):

phase[i, j] = unwrap(phase[i, j]) # unwrap phase first

if deg:

wrap_phase *= math.pi/180.

# Shift the phase if it is below the wrap_phase

phase[i, j] += 2*math.pi * np.maximum(

0, np.ceil((wrap_phase - phase[i, j])/(2*math.pi)))

else:

raise ValueError("wrap_phase must be bool or float.")

# Put the phase back into the original shape

phase = phase.reshape(response.magnitude.shape)

# Save the data for later use (legacy return values)

mag_data.append(response.magnitude)

phase_data.append(phase)

omega_data.append(response.omega)

#

# Process `plot` keyword

#

# We use the `plot` keyword to track legacy usage of `bode_plot`.

# Prior to v0.10, the `bode_plot` command returned mag, phase, and

# omega. Post v0.10, we return an array with the same shape as the

# axes we use for plotting, with each array element containing a list

# of lines drawn on that axes.

#

# There are three possibilities at this stage in the code:

#

# * plot == True: set explicitly by the user. Return mag, phase, omega,

# with a warning.

#

# * plot == False: set explicitly by the user. Return mag, phase,

# omega, with a warning.

#

# * plot == None: this is the new default setting. Return an array of

# lines that were drawn.

#

# If `bode_plot` was called with no `plot` argument and the return

# values were used, the new code will cause problems (you get an array

# of lines instead of magnitude, phase, and frequency). To recover the

# old behavior, call `bode_plot` with `plot=True`.

#

# All of this should be removed in v0.11+ when we get rid of deprecated

# code.

#

if plot is not None:

warnings.warn(

"`bode_plot` return values of mag, phase, omega is deprecated; "

"use frequency_response()", DeprecationWarning)

if plot is False:

# Process the data to match what we were sent

for i in range(len(mag_data)):

mag_data[i] = _process_frequency_response(

data[i], omega_data[i], mag_data[i], squeeze=data[i].squeeze)

phase_data[i] = _process_frequency_response(

data[i], omega_data[i], phase_data[i], squeeze=data[i].squeeze)

if len(data) == 1:

return mag_data[0], phase_data[0], omega_data[0]

else:

return mag_data, phase_data, omega_data

#

# Find/create axes

#

# Data are plotted in a standard subplots array, whose size depends on

# which signals are being plotted and how they are combined. The

# baseline layout for data is to plot everything separately, with

# the magnitude and phase for each output making up the rows and the

# columns corresponding to the different inputs.

#

# Input 0 Input m

# +---------------+ +---------------+

# | mag H_y0,u0 | ... | mag H_y0,um |

# +---------------+ +---------------+

# +---------------+ +---------------+

# | phase H_y0,u0 | ... | phase H_y0,um |

# +---------------+ +---------------+

# : :

# +---------------+ +---------------+

# | mag H_yp,u0 | ... | mag H_yp,um |

# +---------------+ +---------------+

# +---------------+ +---------------+

# | phase H_yp,u0 | ... | phase H_yp,um |

# +---------------+ +---------------+

#

# Several operations are available that change this layout.

#

# * Omitting: either the magnitude or the phase plots can be omitted

# using the plot_magnitude and plot_phase keywords.

#

# * Overlay: inputs and/or outputs can be combined onto a single set of

# axes using the overlay_inputs and overlay_outputs keywords. This

# basically collapses data along either the rows or columns, and a

# legend is generated.

#

# Decide on the maximum number of inputs and outputs

ninputs, noutputs = 0, 0

for response in data: # TODO: make more pythonic/numpic

ninputs = max(ninputs, response.ninputs)

noutputs = max(noutputs, response.noutputs)

# Figure how how many rows and columns to use + offsets for inputs/outputs

if overlay_outputs and overlay_inputs:

nrows = plot_magnitude + plot_phase

ncols = 1

elif overlay_outputs:

nrows = plot_magnitude + plot_phase

ncols = ninputs

elif overlay_inputs:

nrows = (noutputs if plot_magnitude else 0) + \

(noutputs if plot_phase else 0)

ncols = 1

else:

nrows = (noutputs if plot_magnitude else 0) + \

(noutputs if plot_phase else 0)

ncols = ninputs

if ax is None:

# Set up default sharing of axis limits if not specified

for kw in ['share_magnitude', 'share_phase', 'share_frequency']:

if kw not in kwargs or kwargs[kw] is None:

kwargs[kw] = config.defaults['freqplot.' + kw]

fig, ax_array = _process_ax_keyword(

ax, (nrows, ncols), squeeze=False, rcParams=rcParams, clear_text=True)

legend_loc, legend_map, show_legend = _process_legend_keywords(

kwargs, (nrows,ncols), 'center right')

# Get the values for sharing axes limits

share_magnitude = kwargs.pop('share_magnitude', None)

share_phase = kwargs.pop('share_phase', None)

share_frequency = kwargs.pop('share_frequency', None)

# Set up axes variables for easier access below

if plot_magnitude and not plot_phase:

mag_map = np.empty((noutputs, ninputs), dtype=tuple)

for i in range(noutputs):

for j in range(ninputs):

if overlay_outputs and overlay_inputs:

mag_map[i, j] = (0, 0)

elif overlay_outputs:

mag_map[i, j] = (0, j)

elif overlay_inputs:

mag_map[i, j] = (i, 0)

else:

mag_map[i, j] = (i, j)

phase_map = np.full((noutputs, ninputs), None)

share_phase = False

elif plot_phase and not plot_magnitude:

phase_map = np.empty((noutputs, ninputs), dtype=tuple)

for i in range(noutputs):

for j in range(ninputs):

if overlay_outputs and overlay_inputs:

phase_map[i, j] = (0, 0)

elif overlay_outputs:

phase_map[i, j] = (0, j)

elif overlay_inputs:

phase_map[i, j] = (i, 0)

else:

phase_map[i, j] = (i, j)

mag_map = np.full((noutputs, ninputs), None)

share_magnitude = False

else:

mag_map = np.empty((noutputs, ninputs), dtype=tuple)

phase_map = np.empty((noutputs, ninputs), dtype=tuple)

for i in range(noutputs):

for j in range(ninputs):

if overlay_outputs and overlay_inputs:

mag_map[i, j] = (0, 0)

phase_map[i, j] = (1, 0)

elif overlay_outputs:

mag_map[i, j] = (0, j)

phase_map[i, j] = (1, j)

elif overlay_inputs:

mag_map[i, j] = (i*2, 0)

phase_map[i, j] = (i*2 + 1, 0)

else:

mag_map[i, j] = (i*2, j)

phase_map[i, j] = (i*2 + 1, j)

# Identity map needed for setting up shared axes

ax_map = np.empty((nrows, ncols), dtype=tuple)

for i, j in itertools.product(range(nrows), range(ncols)):

ax_map[i, j] = (i, j)

#

# Set up axes limit sharing

#

# This code uses the share_magnitude, share_phase, and share_frequency

# keywords to decide which axes have shared limits and what ticklabels

# to include. The sharing code needs to come before the plots are

# generated, but additional code for removing tick labels needs to come

# *during* and *after* the plots are generated (see below).

#

# Note: if the various share_* keywords are None then a previous set of

# axes are available and no updates should be made.

#

# Utility function to turn off sharing

def _share_axes(ref, share_map, axis):

ref_ax = ax_array[ref]

for index in np.nditer(share_map, flags=["refs_ok"]):

if index.item() == ref:

continue

if axis == 'x':

ax_array[index.item()].sharex(ref_ax)

elif axis == 'y':

ax_array[index.item()].sharey(ref_ax)

else:

raise ValueError("axis must be 'x' or 'y'")

# Process magnitude, phase, and frequency axes

for name, value, map, axis in zip(

['share_magnitude', 'share_phase', 'share_frequency'],

[ share_magnitude, share_phase, share_frequency],

[ mag_map, phase_map, ax_map],

[ 'y', 'y', 'x']):

if value in [True, 'all']:

_share_axes(map[0 if axis == 'y' else -1, 0], map, axis)

elif axis == 'y' and value in ['row']:

for i in range(noutputs if not overlay_outputs else 1):

_share_axes(map[i, 0], map[i], 'y')

elif axis == 'x' and value in ['col']:

for j in range(ncols):

_share_axes(map[-1, j], map[:, j], 'x')

elif value in [False, 'none']:

# TODO: turn off any sharing that is on

pass

elif value is not None:

raise ValueError(

f"unknown value for `{name}`: '{value}'")

#

# Plot the data

#

# The mag_map and phase_map arrays have the indices axes needed for

# making the plots. Labels are used on each axes for later creation of

# legends. The generic labels if of the form:

#

# To output label, From input label, system name

#

# The input and output labels are omitted if overlay_inputs or

# overlay_outputs is False, respectively. The system name is always

# included, since multiple calls to plot() will require a legend that

# distinguishes which system signals are plotted. The system name is

# stripped off later (in the legend-handling code) if it is not needed.

#

# Note: if we are building on top of an existing plot, tick labels

# should be preserved from the existing axes. For log scale axes the

# tick labels seem to appear no matter what => we have to detect if

# they are present at the start and, it not, remove them after calling

# loglog or semilogx.

#

# Create a list of lines for the output

out = np.empty((nrows, ncols), dtype=object)

for i in range(nrows):

for j in range(ncols):

out[i, j] = [] # unique list in each element

# Process label keyword

line_labels = _process_line_labels(label, len(data), ninputs, noutputs)

# Utility function for creating line label

def _make_line_label(response, output_index, input_index):

label = "" # start with an empty label

# Add the output name if it won't appear as an axes label

if noutputs > 1 and overlay_outputs:

label += response.output_labels[output_index]

# Add the input name if it won't appear as a column label

if ninputs > 1 and overlay_inputs:

label += ", " if label != "" else ""

label += response.input_labels[input_index]

# Add the system name (will strip off later if redundant)

label += ", " if label != "" else ""

label += f"{response.sysname}"

return label

for index, response in enumerate(data):

# Get the (pre-processed) data in fully indexed form

mag = mag_data[index].reshape((noutputs, ninputs, -1))

phase = phase_data[index].reshape((noutputs, ninputs, -1))

omega_sys, sysname = omega_data[index], response.sysname

for i, j in itertools.product(range(noutputs), range(ninputs)):

# Get the axes to use for magnitude and phase

ax_mag = ax_array[mag_map[i, j]]

ax_phase = ax_array[phase_map[i, j]]

# Get the frequencies and convert to Hz, if needed

omega_plot = omega_sys / (2 * math.pi) if Hz else omega_sys

if response.isdtime(strict=True):

nyq_freq = (0.5/response.dt) if Hz else (math.pi/response.dt)

# Save the magnitude and phase to plot

mag_plot = 20 * np.log10(mag[i, j]) if dB else mag[i, j]

phase_plot = phase[i, j] * 180. / math.pi if deg else phase[i, j]

# Generate a label

if line_labels is None:

label = _make_line_label(response, i, j)

else:

label = line_labels[index, i, j]

# Magnitude

if plot_magnitude:

pltfcn = ax_mag.semilogx if dB else ax_mag.loglog

# Plot the main data

lines = pltfcn(

omega_plot, mag_plot, *fmt, label=label, **kwargs)

out[mag_map[i, j]] += lines

# Save the information needed for the Nyquist line

if response.isdtime(strict=True):

ax_mag.axvline(

nyq_freq, color=lines[0].get_color(), linestyle='--',

label='_nyq_mag_' + sysname)

# Add a grid to the plot

ax_mag.grid(grid and not display_margins, which='both')

# Phase

if plot_phase:

lines = ax_phase.semilogx(

omega_plot, phase_plot, *fmt, label=label, **kwargs)

out[phase_map[i, j]] += lines

# Save the information needed for the Nyquist line

if response.isdtime(strict=True):

ax_phase.axvline(

nyq_freq, color=lines[0].get_color(), linestyle='--',

label='_nyq_phase_' + sysname)

# Add a grid to the plot

ax_phase.grid(grid and not display_margins, which='both')

#

# Display gain and phase margins (SISO only)

#

if display_margins:

if ninputs > 1 or noutputs > 1:

raise NotImplementedError(

"margins are not available for MIMO systems")

# Compute stability margins for the system

margins = stability_margins(response, method=margins_method)

gm, pm, Wcg, Wcp = (margins[i] for i in [0, 1, 3, 4])

# Figure out sign of the phase at the first gain crossing

# (needed if phase_wrap is True)

phase_at_cp = phase[

0, 0, (np.abs(omega_data[0] - Wcp)).argmin()]

if phase_at_cp >= 0.:

phase_limit = 180.

else:

phase_limit = -180.

if Hz:

Wcg, Wcp = Wcg/(2*math.pi), Wcp/(2*math.pi)

# Draw lines at gain and phase limits

if plot_magnitude:

ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':',

zorder=-20)

mag_ylim = ax_mag.get_ylim()

if plot_phase:

ax_phase.axhline(y=phase_limit if deg else

math.radians(phase_limit),

color='k', linestyle=':', zorder=-20)

phase_ylim = ax_phase.get_ylim()

# Annotate the phase margin (if it exists)

if plot_phase and pm != float('inf') and Wcp != float('nan'):

# Draw dotted lines marking the gain crossover frequencies

if plot_magnitude:

ax_mag.axvline(Wcp, color='k', linestyle=':', zorder=-30)

ax_phase.axvline(Wcp, color='k', linestyle=':', zorder=-30)

# Draw solid segments indicating the margins

if deg:

ax_phase.semilogx(

[Wcp, Wcp], [phase_limit + pm, phase_limit],

color='k', zorder=-20)

else:

ax_phase.semilogx(

[Wcp, Wcp], [math.radians(phase_limit) +

math.radians(pm),

math.radians(phase_limit)],

color='k', zorder=-20)

# Annotate the gain margin (if it exists)

if plot_magnitude and gm != float('inf') and \

Wcg != float('nan'):

# Draw dotted lines marking the phase crossover frequencies

ax_mag.axvline(Wcg, color='k', linestyle=':', zorder=-30)

if plot_phase:

ax_phase.axvline(Wcg, color='k', linestyle=':', zorder=-30)

# Draw solid segments indicating the margins

if dB:

ax_mag.semilogx(

[Wcg, Wcg], [0, -20*np.log10(gm)],

color='k', zorder=-20)

else:

ax_mag.loglog(

[Wcg, Wcg], [1., 1./gm], color='k', zorder=-20)

if display_margins == 'overlay':

# TODO: figure out how to handle case of multiple lines

# Put the margin information in the lower left corner

if plot_magnitude:

ax_mag.text(

0.04, 0.06,

'G.M.: %.2f %s\nFreq: %.2f %s' %

(20*np.log10(gm) if dB else gm,

'dB ' if dB else '',

Wcg, 'Hz' if Hz else 'rad/s'),

horizontalalignment='left',

verticalalignment='bottom',

transform=ax_mag.transAxes,

fontsize=8 if int(mpl.__version__[0]) == 1 else 6)

if plot_phase:

ax_phase.text(

0.04, 0.06,

'P.M.: %.2f %s\nFreq: %.2f %s' %

(pm if deg else math.radians(pm),

'deg' if deg else 'rad',

Wcp, 'Hz' if Hz else 'rad/s'),

horizontalalignment='left',

verticalalignment='bottom',

transform=ax_phase.transAxes,

fontsize=8 if int(mpl.__version__[0]) == 1 else 6)

else:

# Put the title underneath the suptitle (one line per system)

ax = ax_mag if ax_mag else ax_phase

axes_title = ax.get_title()

if axes_title is not None and axes_title != "":

axes_title += "\n"

with plt.rc_context(rcParams):

ax.set_title(

axes_title + f"{sysname}: "

"Gm = %.2f %s(at %.2f %s), "

"Pm = %.2f %s (at %.2f %s)" %

(20*np.log10(gm) if dB else gm,

'dB ' if dB else '',

Wcg, 'Hz' if Hz else 'rad/s',

pm if deg else math.radians(pm),

'deg' if deg else 'rad',

Wcp, 'Hz' if Hz else 'rad/s'))

#

# Finishing handling axes limit sharing

#

# This code handles labels on Bode plots and also removes tick labels

# on shared axes. It needs to come *after* the plots are generated,

# in order to handle two things:

#

# * manually generated labels and grids need to reflect the limits for

# shared axes, which we don't know until we have plotted everything;

#

# * the loglog and semilog functions regenerate the labels (not quite

# sure why, since using sharex and sharey in subplots does not have

# this behavior).

#

# Note: as before, if the various share_* keywords are None then a

# previous set of axes are available and no updates are made. (TODO: true?)

#

for i in range(noutputs):

for j in range(ninputs):

# Utility function to generate phase labels

def gen_zero_centered_series(val_min, val_max, period):

v1 = np.ceil(val_min / period - 0.2)

v2 = np.floor(val_max / period + 0.2)

return np.arange(v1, v2 + 1) * period

# Label the phase axes using multiples of 45 degrees

if plot_phase:

ax_phase = ax_array[phase_map[i, j]]

# Set the labels

if deg:

ylim = ax_phase.get_ylim()

num = np.floor((ylim[1] - ylim[0]) / 45)

factor = max(1, np.round(num / (32 / nrows)) * 2)

ax_phase.set_yticks(gen_zero_centered_series(

ylim[0], ylim[1], 45 * factor))

ax_phase.set_yticks(gen_zero_centered_series(

ylim[0], ylim[1], 15 * factor), minor=True)

else:

ylim = ax_phase.get_ylim()

num = np.ceil((ylim[1] - ylim[0]) / (math.pi/4))

factor = max(1, np.round(num / (36 / nrows)) * 2)

ax_phase.set_yticks(gen_zero_centered_series(

ylim[0], ylim[1], math.pi / 4. * factor))

ax_phase.set_yticks(gen_zero_centered_series(

ylim[0], ylim[1], math.pi / 12. * factor), minor=True)

# Turn off y tick labels for shared axes

for i in range(0, noutputs):

for j in range(1, ncols):

if share_magnitude in [True, 'all', 'row']:

ax_array[mag_map[i, j]].tick_params(labelleft=False)

if share_phase in [True, 'all', 'row']:

ax_array[phase_map[i, j]].tick_params(labelleft=False)

# Turn off x tick labels for shared axes

for i in range(0, nrows-1):

for j in range(0, ncols):

if share_frequency in [True, 'all', 'col']:

ax_array[i, j].tick_params(labelbottom=False)

# If specific omega_limits were given, use them

if omega_limits is not None:

for i, j in itertools.product(range(nrows), range(ncols)):

ax_array[i, j].set_xlim(omega_limits)

#

# Label the axes (including header labels)

#

# Once the data are plotted, we label the axes. The horizontal axes is

# always frequency and this is labeled only on the bottom most row. The

# vertical axes can consist either of a single signal or a combination

# of signals (when overlay_inputs or overlay_outputs is True)

#

# Input/output signals are give at the top of columns and left of rows

# when these are individually plotted.

#

# Label the columns (do this first to get row labels in the right spot)

for j in range(ncols):

# If we have more than one column, label the individual responses

if (noutputs > 1 and not overlay_outputs or ninputs > 1) \

and not overlay_inputs:

with plt.rc_context(rcParams):

ax_array[0, j].set_title(f"From {data[0].input_labels[j]}")

# Label the frequency axis

ax_array[-1, j].set_xlabel(

freq_label.format(units="Hz" if Hz else "rad/s"))

# Label the rows

for i in range(noutputs if not overlay_outputs else 1):

if plot_magnitude:

ax_mag = ax_array[mag_map[i, 0]]

ax_mag.set_ylabel(magnitude_label)

if plot_phase:

ax_phase = ax_array[phase_map[i, 0]]

ax_phase.set_ylabel(phase_label)

if (noutputs > 1 or ninputs > 1) and not overlay_outputs:

if plot_magnitude and plot_phase:

# Get existing ylabel for left column and add a blank line

ax_mag.set_ylabel("\n" + ax_mag.get_ylabel())

ax_phase.set_ylabel("\n" + ax_phase.get_ylabel())

# Find the midpoint between the row axes (+ tight_layout)

_, ypos = _find_axes_center(fig, [ax_mag, ax_phase])

# Get the bounding box including the labels

inv_transform = fig.transFigure.inverted()

mag_bbox = inv_transform.transform(

ax_mag.get_tightbbox(fig.canvas.get_renderer()))

# Figure out location for the text (center left in figure frame)

xpos = mag_bbox[0, 0] # left edge

# Put a centered label as text outside the box

fig.text(

0.8 * xpos, ypos, f"To {data[0].output_labels[i]}\n",

rotation=90, ha='left', va='center',

fontsize=rcParams['axes.titlesize'])

else:

# Only a single axes => add label to the left

ax_array[i, 0].set_ylabel(

f"To {data[0].output_labels[i]}\n" +

ax_array[i, 0].get_ylabel())

#

# Update the plot title (= figure suptitle)

#

# If plots are built up by multiple calls to plot() and the title is

# not given, then the title is updated to provide a list of unique text

# items in each successive title. For data generated by the frequency

# response function this will generate a common prefix followed by a

# list of systems (e.g., "Step response for sys[1], sys[2]").

#

# Set the initial title for the data (unique system names, preserving order)

seen = set()

sysnames = [response.sysname for response in data \

if not (response.sysname in seen or seen.add(response.sysname))]

if ax is None and title is None:

if data[0].title is None:

title = "Bode plot for " + ", ".join(sysnames)

else:

# Allow data to set the title (used by gangof4)

title = data[0].title

_update_plot_title(title, fig, rcParams=rcParams, frame=title_frame)

elif ax is None:

_update_plot_title(

title, fig=fig, rcParams=rcParams, frame=title_frame,

use_existing=False)

#

# Create legends

#

# Legends can be placed manually by passing a legend_map array that

# matches the shape of the suplots, with each item being a string

# indicating the location of the legend for that axes (or None for no

# legend).

#

# If no legend spec is passed, a minimal number of legends are used so

# that each line in each axis can be uniquely identified. The details

# depends on the various plotting parameters, but the general rule is

# to place legends in the top row and right column.

#

# Because plots can be built up by multiple calls to plot(), the legend

# strings are created from the line labels manually. Thus an initial

# call to plot() may not generate any legends (eg, if no signals are

# overlaid), but subsequent calls to plot() will need a legend for each