Skip to content

Commit 0eb4396

Browse files
murrayrmmurrayrm
authored
Merge pull request #449 from sawyerbfuller/call-method
use __call__ instead of evalfr in lti system classes
2 parents c2b00c3 + a631098 commit 0eb4396

File tree

18 files changed

+628
-546
lines changed

18 files changed

+628
-546
lines changed

‎control/bench/time_freqresp.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
sys_tf = tf(sys)
99
w = logspace(-1,1,50)
1010
ntimes = 1000
11-
time_ss = timeit("sys.freqresp(w)", setup="from __main__ import sys, w", number=ntimes)
12-
time_tf = timeit("sys_tf.freqresp(w)", setup="from __main__ import sys_tf, w", number=ntimes)
11+
time_ss = timeit("sys.freqquency_response(w)", setup="from __main__ import sys, w", number=ntimes)
12+
time_tf = timeit("sys_tf.frequency_response(w)", setup="from __main__ import sys_tf, w", number=ntimes)
1313
print("State-space model on %d states: %f" % (nstates, time_ss))
1414
print("Transfer-function model on %d states: %f" % (nstates, time_tf))

‎control/frdata.py‎

Lines changed: 96 additions & 97 deletions
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,7 @@
4848
from warnings import warn
4949
import numpy as np
5050
from numpy import angle, array, empty, ones, \
51-
real, imag, absolute, eye, linalg, where, dot
51+
real, imag, absolute, eye, linalg, where, dot, sort
5252
from scipy.interpolate import splprep, splev
5353
from .lti import LTI
5454

@@ -100,24 +100,18 @@ def __init__(self, *args, **kwargs):
100100
object, other than an FRD, call FRD(sys, omega)
101101
102102
"""
103+
# TODO: discrete-time FRD systems?
103104
smooth = kwargs.get('smooth', False)
104105

105106
if len(args) == 2:
106107
if not isinstance(args[0], FRD) and isinstance(args[0], LTI):
107108
# not an FRD, but still a system, second argument should be
108109
# the frequency range
109110
otherlti = args[0]
110-
self.omega = array(args[1], dtype=float)
111-
self.omega.sort()
111+
self.omega = sort(np.asarray(args[1], dtype=float))
112112
numfreq = len(self.omega)
113-
114113
# calculate frequency response at my points
115-
self.fresp = empty(
116-
(otherlti.outputs, otherlti.inputs, numfreq),
117-
dtype=complex)
118-
for k, w in enumerate(self.omega):
119-
self.fresp[:, :, k] = otherlti._evalfr(w)
120-
114+
self.fresp = otherlti(1j * self.omega, squeeze=False)
121115
else:
122116
# The user provided a response and a freq vector
123117
self.fresp = array(args[0], dtype=complex)
@@ -141,7 +135,7 @@ def __init__(self, *args, **kwargs):
141135
self.omega = args[0].omega
142136
self.fresp = args[0].fresp
143137
else:
144-
raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
138+
raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
145139

146140
# create interpolation functions
147141
if smooth:
@@ -163,17 +157,17 @@ def __str__(self):
163157
mimo = self.inputs > 1 or self.outputs > 1
164158
outstr = ['Frequency response data']
165159

166-
mt, pt, wt = self.freqresp(self.omega)
167160
for i in range(self.inputs):
168161
for j in range(self.outputs):
169162
if mimo:
170163
outstr.append("Input %i to output %i:" % (i + 1, j + 1))
171164
outstr.append('Freq [rad/s] Response')
172165
outstr.append('------------ ---------------------')
173166
outstr.extend(
174-
['%12.3f %10.4g%+10.4gj' % (w, m, p)
175-
for m, p, w in zip(real(self.fresp[j, i, :]),
176-
imag(self.fresp[j, i, :]), wt)])
167+
['%12.3f %10.4g%+10.4gj' % (w, re, im)
168+
for w, re, im in zip(self.omega,
169+
real(self.fresp[j, i, :]),
170+
imag(self.fresp[j, i, :]))])
177171

178172
return '\n'.join(outstr)
179173

@@ -342,110 +336,116 @@ def __pow__(self, other):
342336
return (FRD(ones(self.fresp.shape), self.omega) / self) * \
343337
(self**(other+1))
344338

345-
def evalfr(self, omega):
346-
"""Evaluate a transfer function at a single angular frequency.
347-
348-
self._evalfr(omega) returns the value of the frequency response
349-
at frequency omega.
350-
351-
Note that a "normal" FRD only returns values for which there is an
352-
entry in the omega vector. An interpolating FRD can return
353-
intermediate values.
354-
355-
"""
356-
warn("FRD.evalfr(omega) will be deprecated in a future release "
357-
"of python-control; use sys.eval(omega) instead",
358-
PendingDeprecationWarning) # pragma: no coverage
359-
return self._evalfr(omega)
360-
361339
# Define the `eval` function to evaluate an FRD at a given (real)
362340
# frequency. Note that we choose to use `eval` instead of `evalfr` to
363341
# avoid confusion with :func:`evalfr`, which takes a complex number as its
364342
# argument. Similarly, we don't use `__call__` to avoid confusion between
365343
# G(s) for a transfer function and G(omega) for an FRD object.
366-
def eval(self, omega):
367-
"""Evaluate a transfer function at a single angular frequency.
368-
369-
self.evalfr(omega) returns the value of the frequency response
370-
at frequency omega.
344+
# update Sawyer B. Fuller 2020.08.14: __call__ added to provide a uniform
345+
# interface to systems in general and the lti.frequency_response method
346+
def eval(self, omega, squeeze=True):
347+
"""Evaluate a transfer function at angular frequency omega.
371348
372349
Note that a "normal" FRD only returns values for which there is an
373350
entry in the omega vector. An interpolating FRD can return
374351
intermediate values.
375352
353+
Parameters
354+
----------
355+
omega : float or array_like
356+
Frequencies in radians per second
357+
squeeze : bool, optional (default=True)
358+
If True and `sys` is single input single output (SISO), returns a
359+
1D array rather than a 3D array.
360+
361+
Returns
362+
-------
363+
fresp : (self.outputs, self.inputs, len(x)) or (len(x), ) complex ndarray
364+
The frequency response of the system. Array is ``len(x)`` if and only
365+
if system is SISO and ``squeeze=True``.
376366
"""
377-
return self._evalfr(omega)
378-
379-
# Internal function to evaluate the frequency responses
380-
def _evalfr(self, omega):
381-
"""Evaluate a transfer function at a single angular frequency."""
382-
# Preallocate the output.
383-
if getattr(omega, '__iter__', False):
384-
out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
385-
else:
386-
out = empty((self.outputs, self.inputs), dtype=complex)
367+
omega_array = np.array(omega, ndmin=1) # array-like version of omega
368+
if any(omega_array.imag > 0):
369+
raise ValueError("FRD.eval can only accept real-valued omega")
387370

388371
if self.ifunc is None:
389-
try:
390-
out = self.fresp[:, :, where(self.omega == omega)[0][0]]
391-
except Exception:
372+
elements = np.isin(self.omega, omega) # binary array
373+
if sum(elements) < len(omega_array):
392374
raise ValueError(
393-
"Frequency %f not in frequency list, try an interpolating"
394-
" FRD if you want additional points" % omega)
395-
else:
396-
if getattr(omega, '__iter__', False):
397-
for i in range(self.outputs):
398-
for j in range(self.inputs):
399-
for k, w in enumerate(omega):
400-
frraw = splev(w, self.ifunc[i, j], der=0)
401-
out[i, j, k] = frraw[0] + 1.0j * frraw[1]
375+
"not all frequencies omega are in frequency list of FRD "
376+
"system. Try an interpolating FRD for additional points.")
402377
else:
403-
for i in range(self.outputs):
404-
for j in range(self.inputs):
405-
frraw = splev(omega, self.ifunc[i, j], der=0)
406-
out[i, j] = frraw[0] + 1.0j * frraw[1]
407-
378+
out = self.fresp[:, :, elements]
379+
else:
380+
out = empty((self.outputs, self.inputs, len(omega_array)),
381+
dtype=complex)
382+
for i in range(self.outputs):
383+
for j in range(self.inputs):
384+
for k, w in enumerate(omega_array):
385+
frraw = splev(w, self.ifunc[i, j], der=0)
386+
out[i, j, k] = frraw[0] + 1.0j * frraw[1]
387+
if not hasattr(omega, '__len__'):
388+
# omega is a scalar, squeeze down array along last dim
389+
out = np.squeeze(out, axis=2)
390+
if squeeze and self.issiso():
391+
out = out[0][0]
408392
return out
409393

410-
# Method for generating the frequency response of the system
411-
def freqresp(self, omega):
412-
"""Evaluate the frequency response at a list of angular frequencies.
394+
def __call__(self, s, squeeze=True):
395+
"""Evaluate system's transfer function at complex frequencies.
396+
397+
Returns the complex frequency response `sys(s)` of system `sys` with
398+
`m = sys.inputs` number of inputs and `p = sys.outputs` number of
399+
outputs.
413400
414-
Reports the value of the magnitude, phase, and angular frequency of
415-
the requency response evaluated at omega, where omega is a list of
416-
angular frequencies, and is a sorted version of the input omega.
401+
To evaluate at a frequency omega in radians per second, enter
402+
``s = omega * 1j`` or use ``sys.eval(omega)``
417403
418404
Parameters
419405
----------
420-
omega : array_like
421-
A list of frequencies in radians/sec at which the system should be
422-
evaluated. The list can be either a python list or a numpy array
423-
and will be sorted before evaluation.
406+
s : complex scalar or array_like
407+
Complex frequencies
408+
squeeze : bool, optional (default=True)
409+
If True and `sys` is single input single output (SISO), i.e. `m=1`,
410+
`p=1`, return a 1D array rather than a 3D array.
424411
425412
Returns
426413
-------
427-
mag : (self.outputs, self.inputs, len(omega)) ndarray
428-
The magnitude (absolute value, not dB or log10) of the system
429-
frequency response.
430-
phase : (self.outputs, self.inputs, len(omega)) ndarray
431-
The wrapped phase in radians of the system frequency response.
432-
omega : ndarray or list or tuple
433-
The list of sorted frequencies at which the response was
434-
evaluated.
435-
"""
436-
# Preallocate outputs.
437-
numfreq = len(omega)
438-
mag = empty((self.outputs, self.inputs, numfreq))
439-
phase = empty((self.outputs, self.inputs, numfreq))
414+
fresp : (p, m, len(s)) complex ndarray or (len(s),) complex ndarray
415+
The frequency response of the system. Array is ``(len(s), )`` if
416+
and only if system is SISO and ``squeeze=True``.
440417
441-
omega.sort()
442418
443-
for k, w in enumerate(omega):
444-
fresp = self._evalfr(w)
445-
mag[:, :, k] = abs(fresp)
446-
phase[:, :, k] = angle(fresp)
419+
Raises
420+
------
421+
ValueError
422+
If `s` is not purely imaginary, because
423+
:class:`FrequencyDomainData` systems are only defined at imaginary
424+
frequency values.
425+
"""
426+
if any(abs(np.array(s, ndmin=1).real) > 0):
427+
raise ValueError("__call__: FRD systems can only accept "
428+
"purely imaginary frequencies")
429+
# need to preserve array or scalar status
430+
if hasattr(s, '__len__'):
431+
return self.eval(np.asarray(s).imag, squeeze=squeeze)
432+
else:
433+
return self.eval(complex(s).imag, squeeze=squeeze)
447434

448-
return mag, phase, omega
435+
def freqresp(self, omega):
436+
"""(deprecated) Evaluate transfer function at complex frequencies.
437+
438+
.. deprecated::0.9.0
439+
Method has been given the more pythonic name
440+
:meth:`FrequencyResponseData.frequency_response`. Or use
441+
:func:`freqresp` in the MATLAB compatibility module.
442+
"""
443+
warn("FrequencyResponseData.freqresp(omega) will be removed in a "
444+
"future release of python-control; use "
445+
"FrequencyResponseData.frequency_response(omega), or "
446+
"freqresp(sys, omega) in the MATLAB compatibility module "
447+
"instead", DeprecationWarning)
448+
return self.frequency_response(omega)
449449

450450
def feedback(self, other=1, sign=-1):
451451
"""Feedback interconnection between two FRD objects."""
@@ -515,11 +515,10 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
515515
"Frequency ranges of FRD do not match, conversion not implemented")
516516

517517
elif isinstance(sys, LTI):
518-
omega.sort()
519-
fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
520-
for k, w in enumerate(omega):
521-
fresp[:, :, k] = sys._evalfr(w)
522-
518+
omega = np.sort(omega)
519+
fresp = sys(1j * omega)
520+
if len(fresp.shape) == 1:
521+
fresp = fresp[np.newaxis, np.newaxis, :]
523522
return FRD(fresp, omega, smooth=True)
524523

525524
elif isinstance(sys, (int, float, complex, np.number)):

‎control/freqplot.py‎

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -151,14 +151,14 @@ def bode_plot(syslist, omega=None,
151151
152152
Notes
153153
-----
154-
1. Alternatively, you may use the lower-level method
155-
``(mag, phase, freq) = sys.freqresp(freq)`` to generate the frequency
156-
response for a system, but it returns a MIMO response.
154+
1. Alternatively, you may use the lower-level methods
155+
:meth:`LTI.frequency_response` or ``sys(s)`` or ``sys(z)`` or to
156+
generate the frequency response for a single system.
157157
158158
2. If a discrete time model is given, the frequency response is plotted
159-
along the upper branch of the unit circle, using the mapping z = exp(j
160-
\\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
161-
timebase. If not timebase is specified (dt = True), dt is set to 1.
159+
along the upper branch of the unit circle, using the mapping ``z = exp(1j
160+
* omega * dt)`` where `omega` ranges from 0 to `pi/dt` and `dt` is the discrete
161+
timebase. If timebase not specified (``dt=True``), `dt` is set to 1.
162162
163163
Examples
164164
--------
@@ -228,7 +228,7 @@ def bode_plot(syslist, omega=None,
228228
nyquistfrq = None
229229

230230
# Get the magnitude and phase of the system
231-
mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega_sys)
231+
mag_tmp, phase_tmp, omega_sys = sys.frequency_response(omega_sys)
232232
mag = np.atleast_1d(np.squeeze(mag_tmp))
233233
phase = np.atleast_1d(np.squeeze(phase_tmp))
234234

@@ -588,7 +588,7 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
588588
"Nyquist is currently only implemented for SISO systems.")
589589
else:
590590
# Get the magnitude and phase of the system
591-
mag_tmp, phase_tmp, omega = sys.freqresp(omega)
591+
mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
592592
mag = np.squeeze(mag_tmp)
593593
phase = np.squeeze(phase_tmp)
594594

@@ -718,7 +718,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
718718
omega_plot = omega / (2. * math.pi) if Hz else omega
719719

720720
# TODO: Need to add in the mag = 1 lines
721-
mag_tmp, phase_tmp, omega = S.freqresp(omega)
721+
mag_tmp, phase_tmp, omega = S.frequency_response(omega)
722722
mag = np.squeeze(mag_tmp)
723723
if dB:
724724
plot_axes['s'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -728,7 +728,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
728728
plot_axes['s'].tick_params(labelbottom=False)
729729
plot_axes['s'].grid(grid, which='both')
730730

731-
mag_tmp, phase_tmp, omega = (P * S).freqresp(omega)
731+
mag_tmp, phase_tmp, omega = (P * S).frequency_response(omega)
732732
mag = np.squeeze(mag_tmp)
733733
if dB:
734734
plot_axes['ps'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -738,7 +738,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
738738
plot_axes['ps'].set_ylabel("$|PS|$" + " (dB)" if dB else "")
739739
plot_axes['ps'].grid(grid, which='both')
740740

741-
mag_tmp, phase_tmp, omega = (C * S).freqresp(omega)
741+
mag_tmp, phase_tmp, omega = (C * S).frequency_response(omega)
742742
mag = np.squeeze(mag_tmp)
743743
if dB:
744744
plot_axes['cs'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -749,7 +749,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
749749
plot_axes['cs'].set_ylabel("$|CS|$" + " (dB)" if dB else "")
750750
plot_axes['cs'].grid(grid, which='both')
751751

752-
mag_tmp, phase_tmp, omega = T.freqresp(omega)
752+
mag_tmp, phase_tmp, omega = T.frequency_response(omega)
753753
mag = np.squeeze(mag_tmp)
754754
if dB:
755755
plot_axes['t'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)

0 commit comments

Comments
 (0)