Skip to content

Commit 6f674cf

Browse files
murrayrmmurrayrm
authored
Merge pull request #507 from murrayrm/standardize_squeeze
Standardize squeeze processing in frequency response functions
2 parents b47ed08 + b338e32 commit 6f674cf

File tree

6 files changed

+290
-106
lines changed

6 files changed

+290
-106
lines changed

‎control/config.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,8 @@
1515

1616
# Package level default values
1717
_control_defaults = {
18-
'control.default_dt':0
18+
'control.default_dt': 0,
19+
'control.squeeze_frequency_response': None
1920
}
2021
defaults = dict(_control_defaults)
2122

‎control/frdata.py‎

Lines changed: 48 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,8 @@
5050
from numpy import angle, array, empty, ones, \
5151
real, imag, absolute, eye, linalg, where, dot, sort
5252
from scipy.interpolate import splprep, splev
53-
from .lti import LTI
53+
from .lti import LTI, _process_frequency_response
54+
from . import config
5455

5556
__all__ = ['FrequencyResponseData', 'FRD', 'frd']
5657

@@ -343,7 +344,7 @@ def __pow__(self, other):
343344
# G(s) for a transfer function and G(omega) for an FRD object.
344345
# update Sawyer B. Fuller 2020.08.14: __call__ added to provide a uniform
345346
# interface to systems in general and the lti.frequency_response method
346-
def eval(self, omega, squeeze=True):
347+
def eval(self, omega, squeeze=None):
347348
"""Evaluate a transfer function at angular frequency omega.
348349
349350
Note that a "normal" FRD only returns values for which there is an
@@ -352,19 +353,33 @@ def eval(self, omega, squeeze=True):
352353
353354
Parameters
354355
----------
355-
omega : float or array_like
356+
omega : float or 1D array_like
356357
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.
358+
squeeze : bool, optional
359+
If squeeze=True, remove single-dimensional entries from the shape
360+
of the output even if the system is not SISO. If squeeze=False,
361+
keep all indices (output, input and, if omega is array_like,
362+
frequency) even if the system is SISO. The default value can be
363+
set using config.defaults['control.squeeze_frequency_response'].
360364
361365
Returns
362366
-------
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``.
367+
fresp : complex ndarray
368+
The frequency response of the system. If the system is SISO and
369+
squeeze is not True, the shape of the array matches the shape of
370+
omega. If the system is not SISO or squeeze is False, the first
371+
two dimensions of the array are indices for the output and input
372+
and the remaining dimensions match omega. If ``squeeze`` is True
373+
then single-dimensional axes are removed.
374+
366375
"""
367376
omega_array = np.array(omega, ndmin=1) # array-like version of omega
377+
378+
# Make sure that we are operating on a simple list
379+
if len(omega_array.shape) > 1:
380+
raise ValueError("input list must be 1D")
381+
382+
# Make sure that frequencies are all real-valued
368383
if any(omega_array.imag > 0):
369384
raise ValueError("FRD.eval can only accept real-valued omega")
370385

@@ -384,16 +399,12 @@ def eval(self, omega, squeeze=True):
384399
for k, w in enumerate(omega_array):
385400
frraw = splev(w, self.ifunc[i, j], der=0)
386401
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]
392-
return out
393-
394-
def __call__(self, s, squeeze=True):
402+
403+
return _process_frequency_response(self, omega, out, squeeze=squeeze)
404+
405+
def __call__(self, s, squeeze=None):
395406
"""Evaluate system's transfer function at complex frequencies.
396-
407+
397408
Returns the complex frequency response `sys(s)` of system `sys` with
398409
`m = sys.inputs` number of inputs and `p = sys.outputs` number of
399410
outputs.
@@ -403,18 +414,24 @@ def __call__(self, s, squeeze=True):
403414
404415
Parameters
405416
----------
406-
s : complex scalar or array_like
417+
s : complex scalar or 1D array_like
407418
Complex frequencies
408419
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.
420+
If squeeze=True, remove single-dimensional entries from the shape
421+
of the output even if the system is not SISO. If squeeze=False,
422+
keep all indices (output, input and, if omega is array_like,
423+
frequency) even if the system is SISO. The default value can be
424+
set using config.defaults['control.squeeze_frequency_response'].
411425
412426
Returns
413427
-------
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``.
417-
428+
fresp : complex ndarray
429+
The frequency response of the system. If the system is SISO and
430+
squeeze is not True, the shape of the array matches the shape of
431+
omega. If the system is not SISO or squeeze is False, the first
432+
two dimensions of the array are indices for the output and input
433+
and the remaining dimensions match omega. If ``squeeze`` is True
434+
then single-dimensional axes are removed.
418435
419436
Raises
420437
------
@@ -423,9 +440,14 @@ def __call__(self, s, squeeze=True):
423440
:class:`FrequencyDomainData` systems are only defined at imaginary
424441
frequency values.
425442
"""
426-
if any(abs(np.array(s, ndmin=1).real) > 0):
443+
# Make sure that we are operating on a simple list
444+
if len(np.atleast_1d(s).shape) > 1:
445+
raise ValueError("input list must be 1D")
446+
447+
if any(abs(np.atleast_1d(s).real) > 0):
427448
raise ValueError("__call__: FRD systems can only accept "
428449
"purely imaginary frequencies")
450+
429451
# need to preserve array or scalar status
430452
if hasattr(s, '__len__'):
431453
return self.eval(np.asarray(s).imag, squeeze=squeeze)

‎control/lti.py‎

Lines changed: 85 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
import numpy as np
1616
from numpy import absolute, real, angle, abs
1717
from warnings import warn
18+
from . import config
1819

1920
__all__ = ['issiso', 'timebase', 'common_timebase', 'timebaseEqual',
2021
'isdtime', 'isctime', 'pole', 'zero', 'damp', 'evalfr',
@@ -111,7 +112,7 @@ def damp(self):
111112
Z = -real(splane_poles)/wn
112113
return wn, Z, poles
113114

114-
def frequency_response(self, omega, squeeze=True):
115+
def frequency_response(self, omega, squeeze=None):
115116
"""Evaluate the linear time-invariant system at an array of angular
116117
frequencies.
117118
@@ -124,30 +125,36 @@ def frequency_response(self, omega, squeeze=True):
124125
125126
G(exp(j*omega*dt)) = mag*exp(j*phase).
126127
127-
In general the system may be multiple input, multiple output (MIMO), where
128-
`m = self.inputs` number of inputs and `p = self.outputs` number of
129-
outputs.
128+
In general the system may be multiple input, multiple output (MIMO),
129+
where `m = self.inputs` number of inputs and `p = self.outputs` number
130+
of outputs.
130131
131132
Parameters
132133
----------
133-
omega : float or array_like
134+
omega : float or 1D array_like
134135
A list, tuple, array, or scalar value of frequencies in
135136
radians/sec at which the system will be evaluated.
136-
squeeze : bool, optional (default=True)
137-
If True and the system is single input single output (SISO), i.e. `m=1`,
138-
`p=1`, return a 1D array rather than a 3D array.
137+
squeeze : bool, optional
138+
If squeeze=True, remove single-dimensional entries from the shape
139+
of the output even if the system is not SISO. If squeeze=False,
140+
keep all indices (output, input and, if omega is array_like,
141+
frequency) even if the system is SISO. The default value can be
142+
set using config.defaults['control.squeeze_frequency_response'].
139143
140144
Returns
141145
-------
142-
mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray
146+
mag : ndarray
143147
The magnitude (absolute value, not dB or log10) of the system
144-
frequency response. Array is ``(len(omega), )`` if
145-
and only if system is SISO and ``squeeze=True``.
146-
phase : (p, m, len(omega)) ndarray or (len(omega),) ndarray
148+
frequency response. If the system is SISO and squeeze is not
149+
True, the array is 1D, indexed by frequency. If the system is not
150+
SISO or squeeze is False, the array is 3D, indexed by the output,
151+
input, and frequency. If ``squeeze`` is True then
152+
single-dimensional axes are removed.
153+
phase : ndarray
147154
The wrapped phase in radians of the system frequency response.
148155
omega : ndarray
149156
The (sorted) frequencies at which the response was evaluated.
150-
157+
151158
"""
152159
omega = np.sort(np.array(omega, ndmin=1))
153160
if isdtime(self, strict=True):
@@ -463,9 +470,8 @@ def damp(sys, doprint=True):
463470
(p.real, p.imag, d, w))
464471
return wn, damping, poles
465472

466-
def evalfr(sys, x, squeeze=True):
467-
"""
468-
Evaluate the transfer function of an LTI system for complex frequency x.
473+
def evalfr(sys, x, squeeze=None):
474+
"""Evaluate the transfer function of an LTI system for complex frequency x.
469475
470476
Returns the complex frequency response `sys(x)` where `x` is `s` for
471477
continuous-time systems and `z` for discrete-time systems, with
@@ -481,17 +487,24 @@ def evalfr(sys, x, squeeze=True):
481487
----------
482488
sys: StateSpace or TransferFunction
483489
Linear system
484-
x : complex scalar or array_like
490+
x : complex scalar or 1D array_like
485491
Complex frequency(s)
486492
squeeze : bool, optional (default=True)
487-
If True and `sys` is single input single output (SISO), i.e. `m=1`,
488-
`p=1`, return a 1D array rather than a 3D array.
493+
If squeeze=True, remove single-dimensional entries from the shape of
494+
the output even if the system is not SISO. If squeeze=False, keep all
495+
indices (output, input and, if omega is array_like, frequency) even if
496+
the system is SISO. The default value can be set using
497+
config.defaults['control.squeeze_frequency_response'].
489498
490499
Returns
491500
-------
492-
fresp : (p, m, len(x)) complex ndarray or (len(x),) complex ndarray
493-
The frequency response of the system. Array is ``(len(x), )`` if
494-
and only if system is SISO and ``squeeze=True``.
501+
fresp : complex ndarray
502+
The frequency response of the system. If the system is SISO and
503+
squeeze is not True, the shape of the array matches the shape of
504+
omega. If the system is not SISO or squeeze is False, the first two
505+
dimensions of the array are indices for the output and input and the
506+
remaining dimensions match omega. If ``squeeze`` is True then
507+
single-dimensional axes are removed.
495508
496509
See Also
497510
--------
@@ -511,13 +524,13 @@ def evalfr(sys, x, squeeze=True):
511524
>>> # This is the transfer function matrix evaluated at s = i.
512525
513526
.. todo:: Add example with MIMO system
527+
514528
"""
515529
return sys.__call__(x, squeeze=squeeze)
516530

517-
def freqresp(sys, omega, squeeze=True):
518-
"""
519-
Frequency response of an LTI system at multiple angular frequencies.
520-
531+
def freqresp(sys, omega, squeeze=None):
532+
"""Frequency response of an LTI system at multiple angular frequencies.
533+
521534
In general the system may be multiple input, multiple output (MIMO), where
522535
`m = sys.inputs` number of inputs and `p = sys.outputs` number of
523536
outputs.
@@ -526,22 +539,27 @@ def freqresp(sys, omega, squeeze=True):
526539
----------
527540
sys: StateSpace or TransferFunction
528541
Linear system
529-
omega : float or array_like
542+
omega : float or 1D array_like
530543
A list of frequencies in radians/sec at which the system should be
531544
evaluated. The list can be either a python list or a numpy array
532545
and will be sorted before evaluation.
533-
squeeze : bool, optional (default=True)
534-
If True and `sys` is single input, single output (SISO), returns
535-
1D array rather than a 3D array.
546+
squeeze : bool, optional
547+
If squeeze=True, remove single-dimensional entries from the shape of
548+
the output even if the system is not SISO. If squeeze=False, keep all
549+
indices (output, input and, if omega is array_like, frequency) even if
550+
the system is SISO. The default value can be set using
551+
config.defaults['control.squeeze_frequency_response'].
536552
537553
Returns
538554
-------
539-
mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray
555+
mag : ndarray
540556
The magnitude (absolute value, not dB or log10) of the system
541-
frequency response. Array is ``(len(omega), )`` if and only if system
542-
is SISO and ``squeeze=True``.
543-
544-
phase : (p, m, len(omega)) ndarray or (len(omega),) ndarray
557+
frequency response. If the system is SISO and squeeze is not True,
558+
the array is 1D, indexed by frequency. If the system is not SISO or
559+
squeeze is False, the array is 3D, indexed by the output, input, and
560+
frequency. If ``squeeze`` is True then single-dimensional axes are
561+
removed.
562+
phase : ndarray
545563
The wrapped phase in radians of the system frequency response.
546564
omega : ndarray
547565
The list of sorted frequencies at which the response was
@@ -579,6 +597,7 @@ def freqresp(sys, omega, squeeze=True):
579597
#>>> # input to the 1st output, and the phase (in radians) of the
580598
#>>> # frequency response from the 1st input to the 2nd output, for
581599
#>>> # s = 0.1i, i, 10i.
600+
582601
"""
583602
return sys.frequency_response(omega, squeeze=squeeze)
584603

@@ -593,3 +612,34 @@ def dcgain(sys):
593612
at the origin
594613
"""
595614
return sys.dcgain()
615+
616+
617+
# Process frequency responses in a uniform way
618+
def _process_frequency_response(sys, omega, out, squeeze=None):
619+
# Set value of squeeze argument if not set
620+
if squeeze is None:
621+
squeeze = config.defaults['control.squeeze_frequency_response']
622+
623+
if not hasattr(omega, '__len__'):
624+
# received a scalar x, squeeze down the array along last dim
625+
out = np.squeeze(out, axis=2)
626+
627+
#
628+
# Get rid of unneeded dimensions
629+
#
630+
# There are three possible values for the squeeze keyword at this point:
631+
#
632+
# squeeze=None: squeeze input/output axes iff SISO
633+
# squeeze=True: squeeze all single dimensional axes (ala numpy)
634+
# squeeze-False: don't squeeze any axes
635+
#
636+
if squeeze is True:
637+
# Squeeze everything that we can if that's what the user wants
638+
return np.squeeze(out)
639+
elif squeeze is None and sys.issiso():
640+
# SISO system output squeezed unless explicitly specified otherwise
641+
return out[0][0]
642+
elif squeeze is False or squeeze is None:
643+
return out
644+
else:
645+
raise ValueError("unknown squeeze value")

0 commit comments

Comments
 (0)