Skip to content

Commit 28adf70

Browse files
murrayrmmurrayrm
committed
improvements to nyquist and bode plots
1 parent 6abf6aa commit 28adf70

2 files changed

Lines changed: 98 additions & 37 deletions

File tree

‎src/__init__.py‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,3 +72,11 @@
7272
from timeresp import forced_response, initial_response, step_response, \
7373
impulse_response
7474
from xferfcn import TransferFunction
75+
76+
# Import some of the more common (and benign) MATLAB shortcuts
77+
from matlab import ss, tf, ss2tf, tf2ss, drss
78+
from matlab import pole, zero, evalfr, freqresp, dcgain
79+
from matlab import nichols, rlocus, margin
80+
# bode and nyquist come directly from freqplot.py
81+
from matlab import step, impulse, initial, lsim
82+

‎src/freqplot.py‎

Lines changed: 90 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@
5656

5757
# Bode plot
5858
def bode_plot(syslist, omega=None, dB=False, Hz=False, deg=True,
59-
color=None, Plot=True):
59+
Plot=True, *args, **kwargs):
6060
"""Bode plot for a system
6161
6262
Plots a Bode plot for the system over a (optional) frequency range.
@@ -71,12 +71,12 @@ def bode_plot(syslist, omega=None, dB=False, Hz=False, deg=True,
7171
If True, plot result in dB
7272
Hz : boolean
7373
If True, plot frequency in Hz (omega must be provided in rad/sec)
74-
color : matplotlib color
75-
Color of line in bode plot
7674
deg : boolean
7775
If True, return phase in degrees (else radians)
7876
Plot : boolean
7977
If True, plot magnitude and phase
78+
*args, **kwargs:
79+
Additional options to matplotlib (color, linestyle, etc)
8080
8181
Returns
8282
-------
@@ -95,7 +95,6 @@ def bode_plot(syslist, omega=None, dB=False, Hz=False, deg=True,
9595
9696
Examples
9797
--------
98-
>>> from matlab import ss
9998
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
10099
>>> mag, phase, omega = bode(sys)
101100
"""
@@ -132,53 +131,37 @@ def bode_plot(syslist, omega=None, dB=False, Hz=False, deg=True,
132131
# Magnitude plot
133132
plt.subplot(211);
134133
if dB:
135-
if color==None:
136-
plt.semilogx(omega, mag)
137-
else:
138-
plt.semilogx(omega, mag, color=color)
139-
plt.ylabel("Magnitude (dB)")
134+
plt.semilogx(omega, mag, *args, **kwargs)
140135
else:
141-
if color==None:
142-
plt.loglog(omega, mag)
143-
else:
144-
plt.loglog(omega, mag, color=color)
145-
plt.ylabel("Magnitude")
136+
plt.loglog(omega, mag, *args, **kwargs)
137+
plt.hold(True);
146138

147-
# Add a grid to the plot
139+
# Add a grid to the plot + labeling
148140
plt.grid(True)
149141
plt.grid(True, which='minor')
150-
plt.hold(True);
142+
plt.ylabel("Magnitude (dB)" if dB else "Magnitude")
151143

152144
# Phase plot
153145
plt.subplot(212);
154-
if deg:
155-
phase_deg = phase
156-
else:
157-
phase_deg = phase * 180 / sp.pi
158-
if color==None:
159-
plt.semilogx(omega, phase_deg)
160-
else:
161-
plt.semilogx(omega, phase_deg, color=color)
162-
plt.hold(True)
146+
plt.semilogx(omega, phase, *args, **kwargs)
147+
plt.hold(True);
163148

164-
# Add a grid to the plot
149+
# Add a grid to the plot + labeling
165150
plt.grid(True)
166151
plt.grid(True, which='minor')
167-
plt.ylabel("Phase (deg)")
152+
plt.ylabel("Phase (deg)" if deg else "Phase (rad)")
168153

169154
# Label the frequency axis
170-
if Hz:
171-
plt.xlabel("Frequency (Hz)")
172-
else:
173-
plt.xlabel("Frequency (rad/sec)")
155+
plt.xlabel("Frequency (Hz)" if Hz else "Frequency (rad/sec)")
174156

175157
if len(syslist) == 1:
176158
return mags[0], phases[0], omegas[0]
177159
else:
178160
return mags, phases, omegas
179161

180162
# Nyquist plot
181-
def nyquist_plot(syslist, omega=None, Plot=True):
163+
def nyquist_plot(syslist, omega=None, Plot=True, color='b',
164+
labelFreq=0, *args, **kwargs):
182165
"""Nyquist plot for a system
183166
184167
Plots a Nyquist plot for the system over a (optional) frequency range.
@@ -190,7 +173,11 @@ def nyquist_plot(syslist, omega=None, Plot=True):
190173
omega : freq_range
191174
Range of frequencies (list or bounds) in rad/sec
192175
Plot : boolean
193-
if True, plot magnitude
176+
If True, plot magnitude
177+
labelFreq : int
178+
Label every nth frequency on the plot
179+
*args, **kwargs:
180+
Additional options to matplotlib (color, linestyle, etc)
194181
195182
Returns
196183
-------
@@ -203,7 +190,6 @@ def nyquist_plot(syslist, omega=None, Plot=True):
203190
204191
Examples
205192
--------
206-
>>> from matlab import ss
207193
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
208194
>>> real, imag, freq = nyquist(sys)
209195
"""
@@ -214,12 +200,14 @@ def nyquist_plot(syslist, omega=None, Plot=True):
214200
# Select a default range if none is provided
215201
if (omega == None):
216202
omega = default_frequency_range(syslist)
203+
217204
# Interpolate between wmin and wmax if a tuple or list are provided
218205
elif (isinstance(omega,list) | isinstance(omega,tuple)):
219206
# Only accept tuple or list of length 2
220207
if (len(omega) != 2):
221208
raise ValueError("Supported frequency arguments are (wmin,wmax) tuple or list, or frequency vector. ")
222-
omega = np.logspace(np.log10(omega[0]),np.log10(omega[1]),num=50,endpoint=True,base=10.0)
209+
omega = np.logspace(np.log10(omega[0]), np.log10(omega[1]),
210+
num=50, endpoint=True, base=10.0)
223211
for sys in syslist:
224212
if (sys.inputs > 1 or sys.outputs > 1):
225213
#TODO: Add MIMO nyquist plots.
@@ -236,11 +224,33 @@ def nyquist_plot(syslist, omega=None, Plot=True):
236224

237225
if (Plot):
238226
# Plot the primary curve and mirror image
239-
plt.plot(x, y, '-');
240-
plt.plot(x, -y, '--');
227+
plt.plot(x, y, '-', color=color, *args, **kwargs);
228+
plt.plot(x, -y, '--', color=color, *args, **kwargs);
241229
# Mark the -1 point
242230
plt.plot([-1], [0], 'r+')
243231

232+
# Label the frequencies of the points
233+
if (labelFreq):
234+
for xpt, ypt, omegapt in zip(x, y, omega)[::labelFreq]:
235+
# Convert to Hz
236+
f = omegapt/(2*sp.pi)
237+
238+
# Factor out multiples of 1000 and limit the
239+
# result to the range [-8, 8].
240+
pow1000 = max(min(get_pow1000(f),8),-8)
241+
242+
# Get the SI prefix.
243+
prefix = gen_prefix(pow1000)
244+
245+
# Apply the text. (Use a space before the text to
246+
# prevent overlap with the data.)
247+
#
248+
# np.round() is used because 0.99... appears
249+
# instead of 1.0, and this would otherwise be
250+
# truncated to 0.
251+
plt.text(xpt, ypt,
252+
' ' + str(int(np.round(f/1000**pow1000, 0))) +
253+
' ' + prefix + 'Hz')
244254
return x, y, omega
245255

246256
# Gang of Four
@@ -362,6 +372,49 @@ def default_frequency_range(syslist):
362372

363373
return omega
364374

375+
#
376+
# KLD 5/23/11: Two functions to create nice looking labels
377+
#
378+
def get_pow1000(num):
379+
'''Determine the exponent for which the significand of a number is within the
380+
range [1, 1000).
381+
'''
382+
# Based on algorithm from http://www.mail-archive.com/matplotlib-users@lists.sourceforge.net/msg14433.html, accessed 2010/11/7
383+
# by Jason Heeris 2009/11/18
384+
from decimal import Decimal
385+
from math import floor
386+
dnum = Decimal(str(num))
387+
if dnum == 0:
388+
return 0
389+
elif dnum < 0:
390+
dnum = -dnum
391+
return int(floor(dnum.log10()/3))
392+
393+
def gen_prefix(pow1000):
394+
'''Return the SI prefix for a power of 1000.
395+
'''
396+
# Prefixes according to Table 5 of [BIPM 2006] (excluding hecto,
397+
# deca, deci, and centi).
398+
if pow1000 < -8 or pow1000 > 8:
399+
raise ValueError("Value is out of the range covered by the SI prefixes.")
400+
return ['Y', # yotta (10^24)
401+
'Z', # zetta (10^21)
402+
'E', # exa (10^18)
403+
'P', # peta (10^15)
404+
'T', # tera (10^12)
405+
'G', # giga (10^9)
406+
'M', # mega (10^6)
407+
'k', # kilo (10^3)
408+
'', # (10^0)
409+
'm', # milli (10^-3)
410+
r'$\mu$', # micro (10^-6)
411+
'n', # nano (10^-9)
412+
'p', # pico (10^-12)
413+
'f', # femto (10^-15)
414+
'a', # atto (10^-18)
415+
'z', # zepto (10^-21)
416+
'y'][8 - pow1000] # yocto (10^-24)
417+
365418
# Function aliases
366419
bode = bode_plot
367420
nyquist = nyquist_plot

0 commit comments

Comments
 (0)