Skip to content

Commit 8ebb73b

Browse files
murrayrmmurrayrm
committed
discrete time frequency plots + more consistent timebase methods
1 parent e1e250c commit 8ebb73b

7 files changed

Lines changed: 146 additions & 88 deletions

File tree

‎ChangeLog‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,36 @@
1+
2012-10-18 Richard Murray <murray@dn0a1594a2.sunet>
2+
3+
* src/freqplot.py (gangof4_plot, nyquist_plot): discrete time
4+
systems supported (through fresp)
5+
6+
* src/statesp.py (StateSpace.evalfr): added warning if discrete time
7+
evaluation is above Nyquist frequency
8+
9+
* src/xferfcn.py (TransferFunction.evalfr, freqresp): : added
10+
warning if discrete time evaluation is above Nyquist frequency
11+
12+
* tests/discrete_test.py (TestDiscrete.test_discrete_bode):
13+
rudimentary unit test for discrete bode plot
14+
15+
* src/freqplot.py (bode_plot): fixed bug where input frequency list
16+
was being overwritten by default frequency range
17+
18+
* src/xferfcn.py (TransferFunction.evalfr, freqresp): handles dtime
19+
20+
* src/statesp.py (StateSpace.evalfr): handles dtime
21+
22+
* src/lti.py (timebaseEqual): converted function to take systems as
23+
inputs instead of dt as input (also affected functions in statesp.py
24+
and xferfcn.py)
25+
26+
* src/lti.py (timebase): timebase now returns dt instead of
27+
'ctime'/'dtime'
28+
29+
2012-10-13 Richard Murray <murray@dn0a1594a2.sunet>
30+
31+
* src/lti.py (timebase): new function to return timebase. Doesn't
32+
yet handle the case where dt == None.
33+
134
2012-10-07 Richard Murray <murray@altura.local>
235

336
* src/matlab.py (c2d): MATLAB compatible function (just calls

‎src/ctrlutil.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -96,6 +96,7 @@ def unwrap(angle, period=2*pi):
9696
# Determine if an object is a system
9797
def issys(object):
9898
# Check for a member of one of the classes that we define here
99+
#! TODO: this should probably look for an LTI object instead??
99100
if (isinstance(object, (statesp.StateSpace, xferfcn.TransferFunction))):
100101
return True
101102

‎src/freqplot.py‎

Lines changed: 16 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -44,9 +44,10 @@
4444
import matplotlib.pyplot as plt
4545
import scipy as sp
4646
import numpy as np
47+
from warnings import warn
4748
from ctrlutil import unwrap
4849
from bdalg import feedback
49-
from lti import isdtime
50+
from lti import isdtime, timebaseEqual
5051

5152
#
5253
# Main plotting functions
@@ -90,9 +91,14 @@ def bode_plot(syslist, omega=None, dB=False, Hz=False, deg=True,
9091
9192
Notes
9293
-----
93-
1. Alternatively, you may use the lower-level method
94-
(mag, phase, freq) = sys.freqresp(freq) to generate the frequency
95-
response for a system, but it returns a MIMO response.
94+
1. Alternatively, you may use the lower-level method (mag, phase, freq)
95+
= sys.freqresp(freq) to generate the frequency response for a system,
96+
but it returns a MIMO response.
97+
98+
2. If a discrete time model is given, the frequency response is plotted
99+
along the upper branch of the unit circle, using the mapping z = exp(j
100+
\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
101+
time base. If not timebase is specified (dt = True), dt is set to 1.
96102
97103
Examples
98104
--------
@@ -105,16 +111,12 @@ def bode_plot(syslist, omega=None, dB=False, Hz=False, deg=True,
105111

106112
mags, phases, omegas = [], [], []
107113
for sys in syslist:
108-
# TODO: implement for discrete time systems
109-
if (isdtime(sys, strict=True)):
110-
raise(NotImplementedError("Function not implemented in discrete time"))
111-
112114
if (sys.inputs > 1 or sys.outputs > 1):
113115
#TODO: Add MIMO bode plots.
114116
raise NotImplementedError("Bode is currently only implemented for SISO systems.")
115117
else:
116-
# Select a default range if none is provided
117118
if (omega == None):
119+
# Select a default range if none is provided
118120
omega = default_frequency_range(syslist)
119121

120122
# Get the magnitude and phase of the system
@@ -204,6 +206,7 @@ def nyquist_plot(syslist, omega=None, Plot=True, color='b',
204206

205207
# Select a default range if none is provided
206208
if (omega == None):
209+
#! TODO: think about doing something smarter for discrete
207210
omega = default_frequency_range(syslist)
208211

209212
# Interpolate between wmin and wmax if a tuple or list are provided
@@ -214,10 +217,6 @@ def nyquist_plot(syslist, omega=None, Plot=True, color='b',
214217
omega = np.logspace(np.log10(omega[0]), np.log10(omega[1]),
215218
num=50, endpoint=True, base=10.0)
216219
for sys in syslist:
217-
# TODO: implement for discrete time systems
218-
if (isdtime(sys, strict=True)):
219-
raise(NotImplementedError("Function not implemented in discrete time"))
220-
221220
if (sys.inputs > 1 or sys.outputs > 1):
222221
#TODO: Add MIMO nyquist plots.
223222
raise NotImplementedError("Nyquist is currently only implemented for SISO systems.")
@@ -281,10 +280,6 @@ def gangof4_plot(P, C, omega=None):
281280
-------
282281
None
283282
"""
284-
# TODO: implement for discrete time systems
285-
if (isdtime(P, strict=True) or isdtime(C, strict=True)):
286-
raise(NotImplementedError("Function not implemented in discrete time"))
287-
288283
if (P.inputs > 1 or P.outputs > 1 or C.inputs > 1 or C.outputs >1):
289284
#TODO: Add MIMO go4 plots.
290285
raise NotImplementedError("Gang of four is currently only implemented for SISO systems.")
@@ -365,11 +360,8 @@ def default_frequency_range(syslist):
365360
# detect if single sys passed by checking if it is sequence-like
366361
if (not getattr(syslist, '__iter__', False)):
367362
syslist = (syslist,)
368-
for sys in syslist:
369-
# TODO: implement for discrete time systems
370-
if (isdtime(sys, strict=True)):
371-
raise(NotImplementedError("Function not implemented in discrete time"))
372363

364+
for sys in syslist:
373365
# Add new features to the list
374366
features = np.concatenate((features, np.abs(sys.pole())))
375367
features = np.concatenate((features, np.abs(sys.zero())))
@@ -383,10 +375,12 @@ def default_frequency_range(syslist):
383375
# Take the log of the features
384376
features = np.log10(features)
385377

378+
#! TODO: Add a check in discrete case to make sure we don't get aliasing
379+
386380
# Set the range to be an order of magnitude beyond any features
387381
omega = sp.logspace(np.floor(np.min(features))-1,
388382
np.ceil(np.max(features))+1)
389-
383+
390384
return omega
391385

392386
#

‎src/lti.py‎

Lines changed: 25 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,8 @@
1212
timebaseEqual()
1313
"""
1414

15+
from types import NoneType
16+
1517
class Lti:
1618

1719
"""Lti is a parent class to linear time invariant control (LTI) objects.
@@ -65,36 +67,40 @@ def __init__(self, inputs=1, outputs=1, dt=None):
6567
self.outputs = outputs
6668
self.dt = dt
6769

68-
# Return the timebase of a system
69-
def timebase(sys):
70-
# TODO: add docstring
71-
# If we get passed a constant, timebase is None
70+
# Return the timebase (with conversion if unspecified)
71+
def timebase(sys, strict=True):
72+
"""Return the timebase for an Lti system
73+
74+
dt = timebase(sys)
75+
76+
returns the timebase for a system 'sys'. If the strict option is
77+
set to False, dt = True will be returned as 1.
78+
"""
79+
# System needs to be either a constant or an Lti system
7280
if isinstance(sys, (int, float, long, complex)):
7381
return None
82+
elif not isinstance(sys, Lti):
83+
raise ValueError("Timebase not defined")
7484

75-
# Check for a transfer fucntion or state space object
76-
if isinstance(sys, Lti):
77-
if sys.dt > 0 or sys.dt == True:
78-
return 'dtime';
79-
elif sys.dt == 0:
80-
return 'ctime';
81-
elif sys.dt == None:
82-
return None
85+
# Return the dample time, with converstion to float if strict is false
86+
if (sys.dt == None):
87+
return None
88+
elif (strict):
89+
return float(sys.dt)
8390

84-
# Got pased something we don't recognize or bad timebase
85-
return False;
91+
return sys.dt
8692

8793
# Check to see if two timebases are equal
88-
def timebaseEqual(dt1, dt2):
94+
def timebaseEqual(sys1, sys2):
8995
# TODO: add docstring
90-
if (type(dt1) == bool or type(dt2) == bool):
96+
if (type(sys1.dt) == bool or type(sys2.dt) == bool):
9197
# Make sure both are unspecified discrete timebases
92-
return type(dt1) == type(dt2) and dt1 == dt2
93-
elif (type(dt1) == None or type(dt2) == None):
98+
return type(sys1.dt) == type(sys2.dt) and sys1.dt == sys2.dt
99+
elif (type(sys1.dt) == NoneType or type(sys2.dt) == NoneType):
94100
# One or the other is unspecified => the other can be anything
95101
return True
96102
else:
97-
return dt1 == dt2
103+
return sys1.dt == sys2.dt
98104

99105
# Check to see if a system is a discrete time system
100106
def isdtime(sys, strict=False):

‎src/statesp.py‎

Lines changed: 14 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -81,7 +81,7 @@
8181
from numpy.linalg.linalg import LinAlgError
8282
from scipy.signal import lti
8383
import warnings
84-
from lti import Lti, timebaseEqual
84+
from lti import Lti, timebaseEqual, isdtime
8585
import xferfcn
8686

8787
class StateSpace(Lti):
@@ -249,7 +249,7 @@ def __add__(self, other):
249249
if (self.dt == None and other.dt != None):
250250
dt = other.dt # use dt from second argument
251251
elif (other.dt == None and self.dt != None) or \
252-
(timebaseEqual(self.dt, other.dt)):
252+
(timebaseEqual(self, other)):
253253
dt = self.dt # use dt from first argument
254254
else:
255255
raise ValueError, "Systems have different sampling times"
@@ -307,7 +307,7 @@ def __mul__(self, other):
307307
if (self.dt == None and other.dt != None):
308308
dt = other.dt # use dt from second argument
309309
elif (other.dt == None and self.dt != None) or \
310-
(timebaseEqual(self.dt, other.dt)):
310+
(timebaseEqual(self, other)):
311311
dt = self.dt # use dt from first argument
312312
else:
313313
raise ValueError, "Systems have different sampling times"
@@ -359,12 +359,16 @@ def evalfr(self, omega):
359359
input value s = i * omega.
360360
361361
"""
362-
363-
# TODO: implement for discrete time systems
364-
if (self.dt != 0 and self.dt != None):
365-
raise(NotImplementedError("Function not implemented in discrete time"))
362+
# Figure out the point to evaluate the transfer function
363+
if isdtime(self, strict=True):
364+
dt = timebase(self)
365+
s = exp(1.j * omega * dt)
366+
if (omega * dt > pi):
367+
warn("evalfr: frequency evaluation above Nyquist frequency")
368+
else:
369+
s = omega * 1.j
366370

367-
fresp = self.C * solve(omega * 1.j * eye(self.states) - self.A,
371+
fresp = self.C * solve(s * eye(self.states) - self.A,
368372
self.B) + self.D
369373

370374
return array(fresp)
@@ -382,11 +386,6 @@ def freqresp(self, omega):
382386
input omega.
383387
384388
"""
385-
386-
# TODO: implement for discrete time systems
387-
if (self.dt != 0 and self.dt != None):
388-
raise(NotImplementedError("Function not implemented in discrete time"))
389-
390389
# Preallocate outputs.
391390
numfreq = len(omega)
392391
mag = empty((self.outputs, self.inputs, numfreq))
@@ -395,6 +394,7 @@ def freqresp(self, omega):
395394

396395
omega.sort()
397396

397+
# Evaluate response at each frequency
398398
for k in range(numfreq):
399399
fresp[:, :, k] = self.evalfr(omega[k])
400400

@@ -439,7 +439,7 @@ def feedback(self, other, sign=-1):
439439
if (self.dt == None and other.dt != None):
440440
dt = other.dt # use dt from second argument
441441
elif (other.dt == None and self.dt != None) or \
442-
timebaseEqual(self.dt, other.dt):
442+
timebaseEqual(self, other):
443443
dt = self.dt # use dt from first argument
444444
else:
445445
raise ValueError, "Systems have different sampling times"

‎src/xferfcn.py‎

Lines changed: 26 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -76,10 +76,11 @@
7676

7777
# External function declarations
7878
from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \
79-
polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze
79+
polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze, exp, pi
8080
from scipy.signal import lti
8181
from copy import deepcopy
82-
from lti import Lti, timebaseEqual
82+
from lti import Lti, timebaseEqual, timebase, isdtime
83+
from warnings import warn
8384
import statesp
8485

8586
class TransferFunction(Lti):
@@ -325,7 +326,7 @@ def __add__(self, other):
325326
if (self.dt == None and other.dt != None):
326327
dt = other.dt # use dt from second argument
327328
elif (other.dt == None and self.dt != None) or \
328-
(timebaseEqual(self.dt, other.dt)):
329+
(timebaseEqual(self, other)):
329330
dt = self.dt # use dt from first argument
330331
else:
331332
raise ValueError, "Systems have different sampling times"
@@ -470,16 +471,22 @@ def evalfr(self, omega):
470471
"""
471472

472473
# TODO: implement for discrete time systems
473-
if (self.dt != 0 and self.dt != None):
474-
raise(NotImplementedError("Function not implemented in discrete time"))
474+
if isdtime(self, strict=True):
475+
# Convert the frequency to discrete time
476+
dt = timebase(self)
477+
s = exp(1.j * omega * dt)
478+
if (omega * dt > pi):
479+
warn("evalfr: frequency evaluation above Nyquist frequency")
480+
else:
481+
s = 1.j * omega
475482

476483
# Preallocate the output.
477484
out = empty((self.outputs, self.inputs), dtype=complex)
478485

479486
for i in range(self.outputs):
480487
for j in range(self.inputs):
481-
out[i][j] = (polyval(self.num[i][j], omega * 1.j) /
482-
polyval(self.den[i][j], omega * 1.j))
488+
out[i][j] = (polyval(self.num[i][j], s) /
489+
polyval(self.den[i][j], s))
483490

484491
return out
485492

@@ -495,21 +502,26 @@ def freqresp(self, omega):
495502
496503
"""
497504

498-
# TODO: implement for discrete time systems
499-
if (self.dt != 0 and self.dt != None):
500-
raise(NotImplementedError("Function not implemented in discrete time"))
501-
502505
# Preallocate outputs.
503506
numfreq = len(omega)
504507
mag = empty((self.outputs, self.inputs, numfreq))
505508
phase = empty((self.outputs, self.inputs, numfreq))
506509

507-
omega.sort()
510+
# Figure out the frequencies
511+
omega.sort();
512+
if isdtime(self, strict=True):
513+
dt = timebase(self)
514+
slist = map(lambda w: exp(1.j * w * dt), omega)
515+
if (max(omega) * dt > pi):
516+
warn("evalfr: frequency evaluation above Nyquist frequency")
517+
else:
518+
slist = map(lambda w: 1.j * w, omega)
508519

520+
# Compute frequency response for each input/output pair
509521
for i in range(self.outputs):
510522
for j in range(self.inputs):
511-
fresp = map(lambda w: (polyval(self.num[i][j], w * 1.j) /
512-
polyval(self.den[i][j], w * 1.j)), omega)
523+
fresp = map(lambda s: (polyval(self.num[i][j], s) /
524+
polyval(self.den[i][j], s)), slist)
513525
fresp = array(fresp)
514526

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

0 commit comments

Comments
 (0)