Skip to content

Commit 90da4fb

Browse files
murrayrmmurrayrm
committed
update squeeze processing for freq response + unit tests, doc updates
1 parent 451d6d2 commit 90da4fb

File tree

7 files changed

+230
-105
lines changed

7 files changed

+230
-105
lines changed

‎control/config.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616
# Package level default values
1717
_control_defaults = {
1818
'control.default_dt': 0,
19-
'control.squeeze_frequency_response': True
19+
'control.squeeze_frequency_response': None
2020
}
2121
defaults = dict(_control_defaults)
2222

‎control/frdata.py‎

Lines changed: 38 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -353,25 +353,33 @@ def eval(self, omega, squeeze=None):
353353
354354
Parameters
355355
----------
356-
omega : float or array_like
356+
omega : float or 1D array_like
357357
Frequencies in radians per second
358-
squeeze : bool, optional (default=True)
359-
If True and the system is single-input single-output (SISO),
360-
return a 1D array rather than a 3D array. Default value (True)
361-
set by config.defaults['control.squeeze_frequency_response'].
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'].
362364
363365
Returns
364366
-------
365-
fresp : (self.outputs, self.inputs, len(x)) or (len(x), ) complex ndarray
366-
The frequency response of the system. Array is ``len(x)``
367-
if and only 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.
368374
369375
"""
370-
# Set value of squeeze argument if not set
371-
if squeeze is None:
372-
squeeze = config.defaults['control.squeeze_frequency_response']
373-
374376
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
375383
if any(omega_array.imag > 0):
376384
raise ValueError("FRD.eval can only accept real-valued omega")
377385

@@ -396,7 +404,7 @@ def eval(self, omega, squeeze=None):
396404

397405
def __call__(self, s, squeeze=None):
398406
"""Evaluate system's transfer function at complex frequencies.
399-
407+
400408
Returns the complex frequency response `sys(s)` of system `sys` with
401409
`m = sys.inputs` number of inputs and `p = sys.outputs` number of
402410
outputs.
@@ -406,19 +414,24 @@ def __call__(self, s, squeeze=None):
406414
407415
Parameters
408416
----------
409-
s : complex scalar or array_like
417+
s : complex scalar or 1D array_like
410418
Complex frequencies
411419
squeeze : bool, optional (default=True)
412-
If True and the system is single-input single-output (SISO),
413-
return a 1D array rather than a 3D array. Default value (True)
414-
set by config.defaults['control.squeeze_frequency_response'].
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'].
415425
416426
Returns
417427
-------
418-
fresp : (p, m, len(s)) complex ndarray or (len(s),) complex ndarray
419-
The frequency response of the system. Array is ``(len(s), )`` if
420-
and only if system is SISO and ``squeeze=True``.
421-
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.
422435
423436
Raises
424437
------
@@ -427,13 +440,14 @@ def __call__(self, s, squeeze=None):
427440
:class:`FrequencyDomainData` systems are only defined at imaginary
428441
frequency values.
429442
"""
430-
# Set value of squeeze argument if not set
431-
if squeeze is None:
432-
squeeze = config.defaults['control.squeeze_frequency_response']
443+
# Make sure that we are operating on a simple list
444+
if len(np.array(s, ndmin=1).shape) > 1:
445+
raise ValueError("input list must be 1D")
433446

434447
if any(abs(np.array(s, ndmin=1).real) > 0):
435448
raise ValueError("__call__: FRD systems can only accept "
436449
"purely imaginary frequencies")
450+
437451
# need to preserve array or scalar status
438452
if hasattr(s, '__len__'):
439453
return self.eval(np.asarray(s).imag, squeeze=squeeze)

‎control/lti.py‎

Lines changed: 59 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -131,21 +131,26 @@ def frequency_response(self, omega, squeeze=None):
131131
132132
Parameters
133133
----------
134-
omega : float or array_like
134+
omega : float or 1D array_like
135135
A list, tuple, array, or scalar value of frequencies in
136136
radians/sec at which the system will be evaluated.
137137
squeeze : bool, optional
138-
If True and the system is single-input single-output (SISO),
139-
return a 1D array rather than a 3D array. Default value (True)
140-
set by config.defaults['control.squeeze_frequency_response'].
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'].
141143
142144
Returns
143145
-------
144-
mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray
146+
mag : ndarray
145147
The magnitude (absolute value, not dB or log10) of the system
146-
frequency response. Array is ``(len(omega), )`` if
147-
and only if system is SISO and ``squeeze=True``.
148-
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
149154
The wrapped phase in radians of the system frequency response.
150155
omega : ndarray
151156
The (sorted) frequencies at which the response was evaluated.
@@ -482,18 +487,24 @@ def evalfr(sys, x, squeeze=None):
482487
----------
483488
sys: StateSpace or TransferFunction
484489
Linear system
485-
x : complex scalar or array_like
490+
x : complex scalar or 1D array_like
486491
Complex frequency(s)
487492
squeeze : bool, optional (default=True)
488-
If True and the system is single-input single-output (SISO), return a
489-
1D array rather than a 3D array. Default value (True) set by
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
490497
config.defaults['control.squeeze_frequency_response'].
491498
492499
Returns
493500
-------
494-
fresp : (p, m, len(x)) complex ndarray or (len(x),) complex ndarray
495-
The frequency response of the system. Array is ``(len(x), )`` if
496-
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.
497508
498509
See Also
499510
--------
@@ -519,7 +530,7 @@ def evalfr(sys, x, squeeze=None):
519530

520531
def freqresp(sys, omega, squeeze=None):
521532
"""Frequency response of an LTI system at multiple angular frequencies.
522-
533+
523534
In general the system may be multiple input, multiple output (MIMO), where
524535
`m = sys.inputs` number of inputs and `p = sys.outputs` number of
525536
outputs.
@@ -528,23 +539,27 @@ def freqresp(sys, omega, squeeze=None):
528539
----------
529540
sys: StateSpace or TransferFunction
530541
Linear system
531-
omega : float or array_like
542+
omega : float or 1D array_like
532543
A list of frequencies in radians/sec at which the system should be
533544
evaluated. The list can be either a python list or a numpy array
534545
and will be sorted before evaluation.
535-
squeeze : bool, optional (default=True)
536-
If True and the system is single-input single-output (SISO), return a
537-
1D array rather than a 3D array. Default value (True) set by
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
538551
config.defaults['control.squeeze_frequency_response'].
539552
540553
Returns
541554
-------
542-
mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray
555+
mag : ndarray
543556
The magnitude (absolute value, not dB or log10) of the system
544-
frequency response. Array is ``(len(omega), )`` if and only if system
545-
is SISO and ``squeeze=True``.
546-
547-
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
548563
The wrapped phase in radians of the system frequency response.
549564
omega : ndarray
550565
The list of sorted frequencies at which the response was
@@ -601,14 +616,30 @@ def dcgain(sys):
601616

602617
# Process frequency responses in a uniform way
603618
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+
604623
if not hasattr(omega, '__len__'):
605624
# received a scalar x, squeeze down the array along last dim
606625
out = np.squeeze(out, axis=2)
607626

627+
#
608628
# Get rid of unneeded dimensions
609-
if squeeze is None:
610-
squeeze = config.defaults['control.squeeze_frequency_response']
611-
if squeeze and sys.issiso():
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
612641
return out[0][0]
613-
else:
642+
elif squeeze is False or squeeze is None:
614643
return out
644+
else:
645+
raise ValueError("unknown squeeze value")

‎control/margins.py‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -294,25 +294,25 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
294294
# frequency for gain margin: phase crosses -180 degrees
295295
w_180 = _poly_iw_real_crossing(num_iw, den_iw, epsw)
296296
with np.errstate(all='ignore'): # den=0 is okay
297-
w180_resp = evalfr(sys, 1J * w_180, squeeze=True)
297+
w180_resp = evalfr(sys, 1J * w_180)
298298

299299
# frequency for phase margin : gain crosses magnitude 1
300300
wc = _poly_iw_mag1_crossing(num_iw, den_iw, epsw)
301-
wc_resp = evalfr(sys, 1J * wc, squeeze=True)
301+
wc_resp = evalfr(sys, 1J * wc)
302302

303303
# stability margin
304304
wstab = _poly_iw_wstab(num_iw, den_iw, epsw)
305-
ws_resp = evalfr(sys, 1J * wstab, squeeze=True)
305+
ws_resp = evalfr(sys, 1J * wstab)
306306

307307
else: # Discrete Time
308308
zargs = _poly_z_invz(sys)
309309
# gain margin
310310
z, w_180 = _poly_z_real_crossing(*zargs, epsw=epsw)
311-
w180_resp = evalfr(sys, z, squeeze=True)
311+
w180_resp = evalfr(sys, z)
312312

313313
# phase margin
314314
z, wc = _poly_z_mag1_crossing(*zargs, epsw=epsw)
315-
wc_resp = evalfr(sys, z, squeeze=True)
315+
wc_resp = evalfr(sys, z)
316316

317317
# stability margin
318318
z, wstab = _poly_z_wstab(*zargs, epsw=epsw)
@@ -437,11 +437,11 @@ def phase_crossover_frequencies(sys):
437437
omega = _poly_iw_real_crossing(num_iw, den_iw, 0.)
438438

439439
# using real() to avoid rounding errors and results like 1+0j
440-
gain = np.real(evalfr(sys, 1J * omega, squeeze=True))
440+
gain = np.real(evalfr(sys, 1J * omega))
441441
else:
442442
zargs = _poly_z_invz(sys)
443443
z, omega = _poly_z_real_crossing(*zargs, epsw=0.)
444-
gain = np.real(evalfr(sys, z, squeeze=True))
444+
gain = np.real(evalfr(sys, z))
445445

446446
return omega, gain
447447

‎control/statesp.py‎

Lines changed: 24 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -645,7 +645,7 @@ def __call__(self, x, squeeze=None):
645645
646646
Returns the complex frequency response `sys(x)` where `x` is `s` for
647647
continuous-time systems and `z` for discrete-time systems.
648-
648+
649649
In general the system may be multiple input, multiple output
650650
(MIMO), where `m = self.inputs` number of inputs and `p =
651651
self.outputs` number of outputs.
@@ -657,18 +657,24 @@ def __call__(self, x, squeeze=None):
657657
658658
Parameters
659659
----------
660-
x : complex or complex array_like
660+
x : complex or complex 1D array_like
661661
Complex frequencies
662662
squeeze : bool, optional
663-
If True and the system is single-input single-output (SISO),
664-
return a 1D array rather than a 3D array. Default value (True)
665-
set by config.defaults['control.squeeze_frequency_response'].
663+
If squeeze=True, remove single-dimensional entries from the shape
664+
of the output even if the system is not SISO. If squeeze=False,
665+
keep all indices (output, input and, if omega is array_like,
666+
frequency) even if the system is SISO. The default value can be
667+
set using config.defaults['control.squeeze_frequency_response'].
666668
667669
Returns
668670
-------
669-
fresp : (p, m, len(x)) complex ndarray or (len(x),) complex ndarray
670-
The frequency response of the system. Array is ``len(x)`` if and
671-
only if system is SISO and ``squeeze=True``.
671+
fresp : complex ndarray
672+
The frequency response of the system. If the system is SISO and
673+
squeeze is not True, the shape of the array matches the shape of
674+
omega. If the system is not SISO or squeeze is False, the first
675+
two dimensions of the array are indices for the output and input
676+
and the remaining dimensions match omega. If ``squeeze`` is True
677+
then single-dimensional axes are removed.
672678
673679
"""
674680
# Use Slycot if available
@@ -693,9 +699,13 @@ def slycot_laub(self, x):
693699
Frequency response
694700
"""
695701
from slycot import tb05ad
702+
x_arr = np.atleast_1d(x) # array-like version of x
703+
704+
# Make sure that we are operating on a simple list
705+
if len(x_arr.shape) > 1:
706+
raise ValueError("input list must be 1D")
696707

697708
# preallocate
698-
x_arr = np.atleast_1d(x) # array-like version of x
699709
n = self.states
700710
m = self.inputs
701711
p = self.outputs
@@ -755,6 +765,11 @@ def horner(self, x):
755765
# Fall back because either Slycot unavailable or cannot handle
756766
# certain cases.
757767
x_arr = np.atleast_1d(x) # force to be an array
768+
769+
# Make sure that we are operating on a simple list
770+
if len(x_arr.shape) > 1:
771+
raise ValueError("input list must be 1D")
772+
758773
# Preallocate
759774
out = empty((self.outputs, self.inputs, len(x_arr)), dtype=complex)
760775

0 commit comments

Comments
 (0)