Skip to content

Commit 1a6d806

Browse files
sawyerbfullersawyerbfuller
committed
eliminate _evalfr in lti system classes, replaced with __call__, which evaluates at many frequencies at once
1 parent d3142ff commit 1a6d806

16 files changed

Lines changed: 333 additions & 504 deletions

‎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: 67 additions & 101 deletions
Original file line numberDiff line numberDiff line change
@@ -47,8 +47,8 @@
4747
# External function declarations
4848
from warnings import warn
4949
import numpy as np
50-
from numpy import angle, array, empty, ones, \
51-
real, imag, absolute, eye, linalg, where, dot
50+
from numpy import angle, array, empty, ones, isin, \
51+
real, imag, absolute, eye, linalg, where, dot, sort
5252
from scipy.interpolate import splprep, splev
5353
from .lti import LTI
5454

@@ -107,17 +107,10 @@ def __init__(self, *args, **kwargs):
107107
# not an FRD, but still a system, second argument should be
108108
# the frequency range
109109
otherlti = args[0]
110-
self.omega = array(args[1], dtype=float)
111-
self.omega.sort()
110+
self.omega = sort(np.asarray(args[1], dtype=float))
112111
numfreq = len(self.omega)
113-
114112
# 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-
113+
self.fresp = otherlti(1j * self.omega, squeeze=False)
121114
else:
122115
# The user provided a response and a freq vector
123116
self.fresp = array(args[0], dtype=complex)
@@ -141,7 +134,7 @@ def __init__(self, *args, **kwargs):
141134
self.omega = args[0].omega
142135
self.fresp = args[0].fresp
143136
else:
144-
raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
137+
raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
145138

146139
# create interpolation functions
147140
if smooth:
@@ -163,7 +156,7 @@ def __str__(self):
163156
mimo = self.inputs > 1 or self.outputs > 1
164157
outstr = ['Frequency response data']
165158

166-
mt, pt, wt = self.freqresp(self.omega)
159+
#mt, pt, wt = self.frequency_response(self.omega)
167160
for i in range(self.inputs):
168161
for j in range(self.outputs):
169162
if mimo:
@@ -173,7 +166,7 @@ def __str__(self):
173166
outstr.extend(
174167
['%12.3f %10.4g%+10.4gj' % (w, m, p)
175168
for m, p, w in zip(real(self.fresp[j, i, :]),
176-
imag(self.fresp[j, i, :]), wt)])
169+
imag(self.fresp[j, i, :]), self.omega)])
177170

178171
return '\n'.join(outstr)
179172

@@ -342,110 +335,84 @@ def __pow__(self, other):
342335
return (FRD(ones(self.fresp.shape), self.omega) / self) * \
343336
(self**(other+1))
344337

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-
361338
# Define the `eval` function to evaluate an FRD at a given (real)
362339
# frequency. Note that we choose to use `eval` instead of `evalfr` to
363340
# avoid confusion with :func:`evalfr`, which takes a complex number as its
364341
# argument. Similarly, we don't use `__call__` to avoid confusion between
365342
# G(s) for a transfer function and G(omega) for an FRD object.
366-
def eval(self, omega):
343+
# update Sawyer B. Fuller 2020.08.14: __call__ added to provide a uniform
344+
# interface to systems in general and the lti.frequency_response method
345+
def eval(self, omega, squeeze=True):
367346
"""Evaluate a transfer function at a single angular frequency.
368347
369348
self.evalfr(omega) returns the value of the frequency response
370349
at frequency omega.
371350
372351
Note that a "normal" FRD only returns values for which there is an
373352
entry in the omega vector. An interpolating FRD can return
374-
intermediate values.
353+
intermediate values.
354+
355+
squeeze: bool, optional (default=True)
356+
If True and sys is single input, single output (SISO), return a
357+
1D array or scalar depending on omega's length.
375358
376359
"""
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)
360+
omega_array = np.array(omega, ndmin=1) # array-like version of omega
361+
if any(omega_array.imag > 0):
362+
raise ValueError("FRD.eval can only accept real-valued omega")
387363

388364
if self.ifunc is None:
389-
try:
390-
out = self.fresp[:, :, where(self.omega == omega)[0][0]]
391-
except Exception:
365+
elements = isin(self.omega, omega)
366+
if sum(elements) < len(omega_array):
392367
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]
368+
"not all frequencies omega are in frequency list of FRD "
369+
"system. Try an interpolating FRD for additional points.")
402370
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-
371+
out = self.fresp[:, :, elements]
372+
else:
373+
out = empty((self.outputs, self.inputs, len(omega_array)))
374+
for i in range(self.outputs):
375+
for j in range(self.inputs):
376+
for k, w in enumerate(omega_array):
377+
frraw = splev(w, self.ifunc[i, j], der=0)
378+
out[i, j, k] = frraw[0] + 1.0j * frraw[1]
379+
if not hasattr(omega, '__len__'):
380+
# omega is a scalar, squeeze down array along last dim
381+
out = np.squeeze(out, axis=2)
382+
if squeeze and self.issiso():
383+
out = out[0][0]
408384
return out
409385

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.
413-
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.
417-
418-
Parameters
419-
----------
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.
424-
425-
Returns
426-
-------
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))
440-
441-
omega.sort()
386+
def __call__(self, s, squeeze=True):
387+
"""Evaluate the system's transfer function at complex frequencies s.
442388
443-
for k, w in enumerate(omega):
444-
fresp = self._evalfr(w)
445-
mag[:, :, k] = abs(fresp)
446-
phase[:, :, k] = angle(fresp)
389+
For a SISO system, returns the complex value of the
390+
transfer function. For a MIMO transfer fuction, returns a
391+
matrix of values.
392+
393+
Raises an error if s is not purely imaginary.
394+
395+
squeeze: bool, optional (default=True)
396+
If True and sys is single input, single output (SISO), return a
397+
1D array or scalar depending on omega's length.
447398
448-
return mag, phase, omega
399+
"""
400+
if any(abs(np.array(s, ndmin=1).real) > 0):
401+
raise ValueError("__call__: FRD systems can only accept"
402+
"purely imaginary frequencies")
403+
# need to preserve array or scalar status
404+
if hasattr(s, '__len__'):
405+
return self.eval(np.asarray(s).imag, squeeze=squeeze)
406+
else:
407+
return self.eval(complex(s).imag, squeeze=squeeze)
408+
409+
def freqresp(self, omega):
410+
warn("FrequencyResponseData.freqresp(omega) will be deprecated in a "
411+
"future release of python-control; use "
412+
"FrequencyResponseData.frequency_response(sys, omega), or "
413+
"freqresp(sys, omega) in the MATLAB compatibility module "
414+
"instead", PendingDeprecationWarning)
415+
return self.frequency_response(omega)
449416

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

517484
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-
485+
omega = np.sort(omega)
486+
fresp = sys(1j * omega)
487+
if len(fresp.shape) == 1:
488+
fresp = fresp[np.newaxis, np.newaxis, :]
523489
return FRD(fresp, omega, smooth=True)
524490

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

‎control/freqplot.py‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -136,7 +136,7 @@ def bode_plot(syslist, omega=None,
136136
Notes
137137
-----
138138
1. Alternatively, you may use the lower-level method (mag, phase, freq)
139-
= sys.freqresp(freq) to generate the frequency response for a system,
139+
= sys.frequency_response(freq) to generate the frequency response for a system,
140140
but it returns a MIMO response.
141141
142142
2. If a discrete time model is given, the frequency response is plotted
@@ -207,7 +207,7 @@ def bode_plot(syslist, omega=None,
207207
else:
208208
nyquistfrq = None
209209
# Get the magnitude and phase of the system
210-
mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega_sys)
210+
mag_tmp, phase_tmp, omega_sys = sys.frequency_response(omega_sys)
211211
mag = np.atleast_1d(np.squeeze(mag_tmp))
212212
phase = np.atleast_1d(np.squeeze(phase_tmp))
213213
phase = unwrap(phase)
@@ -524,7 +524,7 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
524524
"Nyquist is currently only implemented for SISO systems.")
525525
else:
526526
# Get the magnitude and phase of the system
527-
mag_tmp, phase_tmp, omega = sys.freqresp(omega)
527+
mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
528528
mag = np.squeeze(mag_tmp)
529529
phase = np.squeeze(phase_tmp)
530530

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

656656
# TODO: Need to add in the mag = 1 lines
657-
mag_tmp, phase_tmp, omega = S.freqresp(omega)
657+
mag_tmp, phase_tmp, omega = S.frequency_response(omega)
658658
mag = np.squeeze(mag_tmp)
659659
if dB:
660660
plot_axes['s'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -664,7 +664,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
664664
plot_axes['s'].tick_params(labelbottom=False)
665665
plot_axes['s'].grid(grid, which='both')
666666

667-
mag_tmp, phase_tmp, omega = (P * S).freqresp(omega)
667+
mag_tmp, phase_tmp, omega = (P * S).frequency_response(omega)
668668
mag = np.squeeze(mag_tmp)
669669
if dB:
670670
plot_axes['ps'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -674,7 +674,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
674674
plot_axes['ps'].set_ylabel("$|PS|$" + " (dB)" if dB else "")
675675
plot_axes['ps'].grid(grid, which='both')
676676

677-
mag_tmp, phase_tmp, omega = (C * S).freqresp(omega)
677+
mag_tmp, phase_tmp, omega = (C * S).frequency_response(omega)
678678
mag = np.squeeze(mag_tmp)
679679
if dB:
680680
plot_axes['cs'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
@@ -685,7 +685,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
685685
plot_axes['cs'].set_ylabel("$|CS|$" + " (dB)" if dB else "")
686686
plot_axes['cs'].grid(grid, which='both')
687687

688-
mag_tmp, phase_tmp, omega = T.freqresp(omega)
688+
mag_tmp, phase_tmp, omega = T.frequency_response(omega)
689689
mag = np.squeeze(mag_tmp)
690690
if dB:
691691
plot_axes['t'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)

0 commit comments

Comments
 (0)