"""TransferFunction(num, den[, dt])

def __init__(self, *args):

"""TransferFunction(num, den[, dt])

args = deepcopy(args)

if len(args) == 2:

# The user provided a numerator and a denominator.

(num, den) = args

dt = config.defaults['xferfcn.default_dt']

elif len(args) == 3:

# Discrete time transfer function

(num, den, dt) = args

elif len(args) == 1:

# Use the copy constructor.

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

raise TypeError("The one-argument constructor can only take \

% type(args[0]))

num = args[0].num

den = args[0].den

# TODO: not sure this can ever happen since dt is always present

try:

dt = args[0].dt

except NameError: # pragma: no coverage

dt = config.defaults['xferfcn.default_dt']

else:

raise ValueError("Needs 1, 2 or 3 arguments; received %i."

% len(args))

num = _clean_part(num)

den = _clean_part(den)

inputs = len(num[0])

outputs = len(num)

# Make sure numerator and denominator matrices have consistent sizes

if inputs != len(den[0]):

raise ValueError(

"The numerator has %i input(s), but the denominator has "

"%i input(s)." % (inputs, len(den[0])))

if outputs != len(den):

raise ValueError(

"The numerator has %i output(s), but the denominator has "

"%i output(s)." % (outputs, len(den)))

# Additional checks/updates on structure of the transfer function

for i in range(outputs):

# Make sure that each row has the same number of columns

if len(num[i]) != inputs:

raise ValueError(

"Row 0 of the numerator matrix has %i elements, but row "

"%i has %i." % (inputs, i, len(num[i])))

if len(den[i]) != inputs:

raise ValueError(

"Row 0 of the denominator matrix has %i elements, but row "

"%i has %i." % (inputs, i, len(den[i])))

# Check for zeros in numerator or denominator

# TODO: Right now these checks are only done during construction.

# It might be worthwhile to think of a way to perform checks if the

# user modifies the transfer function after construction.

for j in range(inputs):

# Check that we don't have any zero denominators.

zeroden = True

for k in den[i][j]:

if k:

zeroden = False

break

if zeroden:

raise ValueError(

"Input %i, output %i has a zero denominator."

% (j + 1, i + 1))

# If we have zero numerators, set the denominator to 1.

zeronum = True

for k in num[i][j]:

if k:

zeronum = False

break

if zeronum:

den[i][j] = ones(1)

LTI.__init__(self, inputs, outputs, dt)

self.num = num

self.den = den

self._truncatecoeff()

def __call__(self, s):

"""Evaluate the system's transfer function for a complex variable

if self.issiso():

# return a scalar

return self.horner(s)[0][0]

else:

# return a matrix

return self.horner(s)

def _truncatecoeff(self):

"""Remove extraneous zero coefficients from num and den.

# Beware: this is a shallow copy. This should be okay.

data = [self.num, self.den]

for p in range(len(data)):

for i in range(self.outputs):

for j in range(self.inputs):

# Find the first nontrivial coefficient.

nonzero = None

for k in range(data[p][i][j].size):

if data[p][i][j][k]:

nonzero = k

break

if nonzero is None:

# The array is all zeros.

data[p][i][j] = zeros(1)

else:

# Truncate the trivial coefficients.

data[p][i][j] = data[p][i][j][nonzero:]

[self.num, self.den] = data

def __str__(self, var=None):

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

mimo = self.inputs > 1 or self.outputs > 1

if var is None:

# TODO: replace with standard calls to lti functions

var = 's' if self.dt is None or self.dt == 0 else 'z'

outstr = ""

for i in range(self.inputs):

for j in range(self.outputs):

if mimo:

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

# Convert the numerator and denominator polynomials to strings.

numstr = _tf_polynomial_to_string(self.num[j][i], var=var)

denstr = _tf_polynomial_to_string(self.den[j][i], var=var)

# Figure out the length of the separating line

dashcount = max(len(numstr), len(denstr))

dashes = '-' * dashcount

# Center the numerator or denominator

if len(numstr) < dashcount:

numstr = ' ' * ((dashcount - len(numstr)) // 2) + numstr

if len(denstr) < dashcount:

denstr = ' ' * ((dashcount - len(denstr)) // 2) + denstr

outstr += "\n" + numstr + "\n" + dashes + "\n" + denstr + "\n"

# See if this is a discrete time system with specific sampling time

if not (self.dt is None) and type(self.dt) != bool and self.dt > 0:

# TODO: replace with standard calls to lti functions

outstr += "\ndt = " + self.dt.__str__() + "\n"

return outstr

# represent to implement a re-loadable version

def __repr__(self):

"""Print transfer function in loadable form"""

if self.issiso():

return "TransferFunction({num}, {den}{dt})".format(

num=self.num[0][0].__repr__(), den=self.den[0][0].__repr__(),

dt=(isdtime(self, strict=True) and ', {}'.format(self.dt)) or '')

else:

return "TransferFunction({num}, {den}{dt})".format(

num=self.num.__repr__(), den=self.den.__repr__(),

dt=(isdtime(self, strict=True) and ', {}'.format(self.dt)) or '')

def _repr_latex_(self, var=None):

"""LaTeX representation of transfer function, for Jupyter notebook"""

mimo = self.inputs > 1 or self.outputs > 1

if var is None:

# ! TODO: replace with standard calls to lti functions

var = 's' if self.dt is None or self.dt == 0 else 'z'

out = ['$$']

if mimo:

out.append(r"\begin{bmatrix}")

for i in range(self.outputs):

for j in range(self.inputs):

# Convert the numerator and denominator polynomials to strings.

numstr = _tf_polynomial_to_string(self.num[i][j], var=var)

denstr = _tf_polynomial_to_string(self.den[i][j], var=var)

numstr = _tf_string_to_latex(numstr, var=var)

denstr = _tf_string_to_latex(denstr, var=var)

out += [r"\frac{", numstr, "}{", denstr, "}"]

if mimo and j < self.outputs - 1:

out.append("&")

if mimo:

out.append(r"\\")

if mimo:

out.append(r" \end{bmatrix}")

# See if this is a discrete time system with specific sampling time

if not (self.dt is None) and type(self.dt) != bool and self.dt > 0:

out += [r"\quad dt = ", str(self.dt)]

out.append("$$")

return ''.join(out)

def __neg__(self):

"""Negate a transfer function."""

num = deepcopy(self.num)

for i in range(self.outputs):

for j in range(self.inputs):

num[i][j] *= -1

return TransferFunction(num, self.den, self.dt)

def __add__(self, other):

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

from .statesp import StateSpace

# Convert the second argument to a transfer function.

if isinstance(other, StateSpace):

other = _convert_to_transfer_function(other)

elif not isinstance(other, TransferFunction):

other = _convert_to_transfer_function(other, inputs=self.inputs,

outputs=self.outputs)

# Check that the input-output sizes are consistent.

if self.inputs != other.inputs:

raise ValueError(

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

% (self.inputs, other.inputs))

if self.outputs != other.outputs:

raise ValueError(

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

% (self.outputs, other.outputs))

# Figure out the sampling time to use

if self.dt is None and other.dt is not None:

dt = other.dt # use dt from second argument

elif (other.dt is None and self.dt is not None) or \

(timebaseEqual(self, other)):

dt = self.dt # use dt from first argument

else:

raise ValueError("Systems have different sampling times")

# Preallocate the numerator and denominator of the sum.

num = [[[] for j in range(self.inputs)] for i in range(self.outputs)]

den = [[[] for j in range(self.inputs)] for i in range(self.outputs)]

for i in range(self.outputs):

for j in range(self.inputs):

num[i][j], den[i][j] = _add_siso(

self.num[i][j], self.den[i][j],

other.num[i][j], other.den[i][j])

return TransferFunction(num, den, dt)

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)):

other = _convert_to_transfer_function(other, inputs=self.inputs,

outputs=self.inputs)

else:

other = _convert_to_transfer_function(other)

# Check that the input-output sizes are consistent.

if self.inputs != other.outputs:

raise ValueError(

"C = A * B: A has %i column(s) (input(s)), but B has %i "

"row(s)\n(output(s))." % (self.inputs, other.outputs))

inputs = other.inputs

outputs = self.outputs

# Figure out the sampling time to use

if self.dt is None and other.dt is not None:

dt = other.dt # use dt from second argument

elif (other.dt is None and self.dt is not None) or \

(self.dt == other.dt):

dt = self.dt # use dt from first argument

else:

raise ValueError("Systems have different sampling times")

# Preallocate the numerator and denominator of the sum.

num = [[[0] for j in range(inputs)] for i in range(outputs)]

den = [[[1] for j in range(inputs)] for i in range(outputs)]

# Temporary storage for the summands needed to find the (i, j)th

# element of the product.

num_summand = [[] for k in range(self.inputs)]

den_summand = [[] for k in range(self.inputs)]

# Multiply & add.

for row in range(outputs):

for col in range(inputs):

for k in range(self.inputs):

num_summand[k] = polymul(

self.num[row][k], other.num[k][col])

den_summand[k] = polymul(

self.den[row][k], other.den[k][col])

num[row][col], den[row][col] = _add_siso(

num[row][col], den[row][col],

num_summand[k], den_summand[k])

return TransferFunction(num, den, dt)

def __rmul__(self, other):

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

# Convert the second argument to a transfer function.

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

other = _convert_to_transfer_function(other, inputs=self.inputs,

outputs=self.inputs)

else:

other = _convert_to_transfer_function(other)

# Check that the input-output sizes are consistent.

if other.inputs != self.outputs:

raise ValueError(

"C = A * B: A has %i column(s) (input(s)), but B has %i "

"row(s)\n(output(s))." % (other.inputs, self.outputs))

inputs = self.inputs

outputs = other.outputs

# Figure out the sampling time to use

if self.dt is None and other.dt is not None:

dt = other.dt # use dt from second argument

elif (other.dt is None and self.dt is not None) \

or (self.dt == other.dt):

dt = self.dt # use dt from first argument

else:

raise ValueError("Systems have different sampling times")

# Preallocate the numerator and denominator of the sum.

num = [[[0] for j in range(inputs)] for i in range(outputs)]

den = [[[1] for j in range(inputs)] for i in range(outputs)]

# Temporary storage for the summands needed to find the

# (i, j)th element

# of the product.

num_summand = [[] for k in range(other.inputs)]

den_summand = [[] for k in range(other.inputs)]

for i in range(outputs): # Iterate through rows of product.

for j in range(inputs): # Iterate through columns of product.

for k in range(other.inputs): # Multiply & add.

num_summand[k] = polymul(other.num[i][k], self.num[k][j])

den_summand[k] = polymul(other.den[i][k], self.den[k][j])

num[i][j], den[i][j] = _add_siso(

num[i][j], den[i][j],

num_summand[k], den_summand[k])

return TransferFunction(num, den, dt)

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

other = _convert_to_transfer_function(

other, inputs=self.inputs,

outputs=self.inputs)

else:

other = _convert_to_transfer_function(other)

if (self.inputs > 1 or self.outputs > 1 or

other.inputs > 1 or other.outputs > 1):

raise NotImplementedError(

"TransferFunction.__truediv__ is currently \

implemented only for SISO systems.")

# Figure out the sampling time to use

if self.dt is None and other.dt is not None:

dt = other.dt # use dt from second argument

elif (other.dt is None and self.dt is not None) or \

(self.dt == other.dt):

dt = self.dt # use dt from first argument

else:

raise ValueError("Systems have different sampling times")

num = polymul(self.num[0][0], other.den[0][0])

den = polymul(self.den[0][0], other.num[0][0])

return TransferFunction(num, den, dt)

# TODO: Remove when transition to python3 complete

def __div__(self, other):

return TransferFunction.__truediv__(self, other)

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

other = _convert_to_transfer_function(

other, inputs=self.inputs,

outputs=self.inputs)

else:

other = _convert_to_transfer_function(other)

if (self.inputs > 1 or self.outputs > 1 or

other.inputs > 1 or other.outputs > 1):

raise NotImplementedError(

"TransferFunction.__rtruediv__ is currently implemented only "

"for SISO systems.")

return other / self

# TODO: Remove when transition to python3 complete

def __rdiv__(self, other):

return TransferFunction.__rtruediv__(self, other)

def __pow__(self, other):

if not type(other) == int:

raise ValueError("Exponent must be an integer")

if other == 0:

return TransferFunction([1], [1]) # unity

if other > 0:

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

if other < 0:

return (TransferFunction([1], [1]) / self) * (self**(other + 1))

def __getitem__(self, key):

key1, key2 = key

# pre-process

if isinstance(key1, int):

key1 = slice(key1, key1 + 1, 1)

if isinstance(key2, int):

key2 = slice(key2, key2 + 1, 1)

# dim1

start1, stop1, step1 = key1.start, key1.stop, key1.step

if step1 is None:

step1 = 1

if start1 is None:

start1 = 0

if stop1 is None:

stop1 = len(self.num)

# dim1

start2, stop2, step2 = key2.start, key2.stop, key2.step

if step2 is None:

step2 = 1

if start2 is None:

start2 = 0

if stop2 is None:

stop2 = len(self.num[0])

num = []

den = []

for i in range(start1, stop1, step1):

num_i = []

den_i = []

for j in range(start2, stop2, step2):

num_i.append(self.num[i][j])

den_i.append(self.den[i][j])

num.append(num_i)

den.append(den_i)

if self.isctime():

return TransferFunction(num, den)

else:

return TransferFunction(num, den, self.dt)

def evalfr(self, omega):

"""Evaluate a transfer function at a single angular frequency.

warn("TransferFunction.evalfr(omega) will be deprecated in a "

"future release of python-control; use evalfr(sys, omega*1j) "

"instead", PendingDeprecationWarning)

return self._evalfr(omega)

def _evalfr(self, omega):

"""Evaluate a transfer function at a single angular frequency."""

# TODO: implement for discrete time systems

if isdtime(self, strict=True):

# Convert the frequency to discrete time

dt = timebase(self)

s = exp(1.j * omega * dt)

if np.any(omega * dt > pi):

warn("_evalfr: frequency evaluation above Nyquist frequency")

else:

s = 1.j * omega

return self.horner(s)

def horner(self, s):

"""Evaluate the systems's transfer function for a complex variable

# Preallocate the output.

if getattr(s, '__iter__', False):

out = empty((self.outputs, self.inputs, len(s)), dtype=complex)

else:

out = empty((self.outputs, self.inputs), dtype=complex)

for i in range(self.outputs):

for j in range(self.inputs):

out[i][j] = (polyval(self.num[i][j], s) /

polyval(self.den[i][j], s))

return out

def freqresp(self, omega):

"""Evaluate the transfer function at a list of angular frequencies.

# Preallocate outputs.

numfreq = len(omega)

mag = empty((self.outputs, self.inputs, numfreq))

phase = empty((self.outputs, self.inputs, numfreq))

# Figure out the frequencies

omega.sort()

if isdtime(self, strict=True):

dt = timebase(self)

slist = np.array([exp(1.j * w * dt) for w in omega])

if max(omega) * dt > pi:

warn("freqresp: frequency evaluation above Nyquist frequency")

else:

slist = np.array([1j * w for w in omega])

# Compute frequency response for each input/output pair

for i in range(self.outputs):

for j in range(self.inputs):

fresp = (polyval(self.num[i][j], slist) /

polyval(self.den[i][j], slist))

mag[i, j, :] = abs(fresp)

phase[i, j, :] = angle(fresp)

return mag, phase, omega

def pole(self):

"""Compute the poles of a transfer function."""

_, den, denorder = self._common_den(allow_nonproper=True)

rts = []

for d, o in zip(den, denorder):

rts.extend(roots(d[:o + 1]))

return np.array(rts)

def zero(self):

"""Compute the zeros of a transfer function."""

if self.inputs > 1 or self.outputs > 1:

raise NotImplementedError(

"TransferFunction.zero is currently only implemented "

"for SISO systems.")

else:

# for now, just give zeros of a SISO tf

return roots(self.num[0][0])

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

"""Feedback interconnection between two LTI objects."""

other = _convert_to_transfer_function(other)

if (self.inputs > 1 or self.outputs > 1 or

other.inputs > 1 or other.outputs > 1):

# TODO: MIMO feedback

raise NotImplementedError(

"TransferFunction.feedback is currently only implemented "

"for SISO functions.")

# Figure out the sampling time to use

if self.dt is None and other.dt is not None:

dt = other.dt # use dt from second argument

elif (other.dt is None and self.dt is not None) or \

(self.dt == other.dt):

dt = self.dt # use dt from first argument

else:

raise ValueError("Systems have different sampling times")

num1 = self.num[0][0]

den1 = self.den[0][0]

num2 = other.num[0][0]

den2 = other.den[0][0]

num = polymul(num1, den2)

den = polyadd(polymul(den2, den1), -sign * polymul(num2, num1))

return TransferFunction(num, den, dt)

# For MIMO or SISO systems, the analytic expression is

# self / (1 - sign * other * self)

# But this does not work correctly because the state size will be too

# large.

def minreal(self, tol=None):

"""Remove cancelling pole/zero pairs from a transfer function"""

# based on octave minreal

# default accuracy

from sys import float_info

sqrt_eps = sqrt(float_info.epsilon)

# pre-allocate arrays

num = [[[] for j in range(self.inputs)] for i in range(self.outputs)]

den = [[[] for j in range(self.inputs)] for i in range(self.outputs)]

for i in range(self.outputs):

for j in range(self.inputs):

# split up in zeros, poles and gain

newzeros = []

zeros = roots(self.num[i][j])

poles = roots(self.den[i][j])

gain = self.num[i][j][0] / self.den[i][j][0]

# check all zeros

for z in zeros:

t = tol or \

1000 * max(float_info.epsilon, abs(z) * sqrt_eps)

idx = where(abs(z - poles) < t)[0]

if len(idx):

# cancel this zero against one of the poles

poles = delete(poles, idx[0])

else:

# keep this zero

newzeros.append(z)

# poly([]) returns a scalar, but we always want a 1d array

num[i][j] = np.atleast_1d(gain * real(poly(newzeros)))

den[i][j] = np.atleast_1d(real(poly(poles)))

# end result

return TransferFunction(num, den, self.dt)

def returnScipySignalLTI(self):

"""Return a list of a list of :class:`scipy.signal.lti` objects.

# TODO: implement for discrete time systems

if self.dt != 0 and self.dt is not None:

raise NotImplementedError("Function not \

implemented in discrete time")

# Preallocate the output.

out = [[[] for j in range(self.inputs)] for i in range(self.outputs)]

for i in range(self.outputs):

for j in range(self.inputs):

out[i][j] = lti(self.num[i][j], self.den[i][j])

return out

def _common_den(self, imag_tol=None, allow_nonproper=False):

"""

# Machine precision for floats.

eps = finfo(float).eps

real_tol = sqrt(eps * self.inputs * self.outputs)

# The tolerance to use in deciding if a pole is complex

if (imag_tol is None):

imag_tol = 2 * real_tol

# A list to keep track of cumulative poles found as we scan

# self.den[..][..]

poles = [[] for j in range(self.inputs)]

# RvP, new implementation 180526, issue #194

# BG, modification, issue #343, PR #354

# pre-calculate the poles for all num, den

# has zeros, poles, gain, list for pole indices not in den,

# number of poles known at the time analyzed

# do not calculate minreal. Rory's hint .minreal()

poleset = []

for i in range(self.outputs):

poleset.append([])

for j in range(self.inputs):

if abs(self.num[i][j]).max() <= eps:

poleset[-1].append([array([], dtype=float),

roots(self.den[i][j]), 0.0, [], 0])

else:

z, p, k = tf2zpk(self.num[i][j], self.den[i][j])

poleset[-1].append([z, p, k, [], 0])

# collect all individual poles

for j in range(self.inputs):

for i in range(self.outputs):

currentpoles = poleset[i][j][1]

nothave = ones(currentpoles.shape, dtype=bool)

for ip, p in enumerate(poles[j]):

collect = (np.isclose(currentpoles.real, p.real,

atol=real_tol) &

np.isclose(currentpoles.imag, p.imag,

atol=imag_tol) &

nothave)

if np.any(collect):

# mark first found pole as already collected

nothave[nonzero(collect)[0][0]] = False

else:

# remember id of pole not in tf

poleset[i][j][3].append(ip)

for h, c in zip(nothave, currentpoles):

if h:

if abs(c.imag) < imag_tol:

c = c.real

poles[j].append(c)

# remember how many poles now known

poleset[i][j][4] = len(poles[j])

# figure out maximum number of poles, for sizing the den

maxindex = max([len(p) for p in poles])

den = zeros((self.inputs, maxindex + 1), dtype=float)

num = zeros((max(1, self.outputs, self.inputs),

max(1, self.outputs, self.inputs),

maxindex + 1),

dtype=float)

denorder = zeros((self.inputs,), dtype=int)

havenonproper = False

for j in range(self.inputs):

if not len(poles[j]):

# no poles matching this input; only one or more gains

den[j, 0] = 1.0

for i in range(self.outputs):

num[i, j, 0] = poleset[i][j][2]

else:

# create the denominator matching this input

# coefficients should be padded on right, ending at maxindex

maxindex = len(poles[j])

den[j, :maxindex+1] = poly(poles[j])

denorder[j] = maxindex

# now create the numerator, also padded on the right

for i in range(self.outputs):

# start with the current set of zeros for this output

nwzeros = list(poleset[i][j][0])

# add all poles not found in the original denominator,

# and the ones later added from other denominators

for ip in chain(poleset[i][j][3],

range(poleset[i][j][4], maxindex)):

nwzeros.append(poles[j][ip])

numpoly = poleset[i][j][2] * np.atleast_1d(poly(nwzeros))

# td04ad expects a proper transfer function. If the

# numerater has a higher order than the denominator, the

# padding will fail

if len(numpoly) > maxindex + 1:

if allow_nonproper:

havenonproper = True

break

raise ValueError(

self.__str__() +

"is not a proper transfer function. "

"The degree of the numerators must not exceed "

"the degree of the denominators.")

# numerator polynomial should be padded on left and right

# ending at maxindex to line up with what td04ad expects.

num[i, j, maxindex+1-len(numpoly):maxindex+1] = numpoly

# print(num[i, j])

if havenonproper:

num = None

return num, den, denorder

def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):

"""Convert a continuous-time system to discrete time