Skip to content

Commit 27b9e9c

Browse files
repaghrepagh
committed
some progress made
1 parent dc4727d commit 27b9e9c

File tree

2 files changed

+125
-129
lines changed

2 files changed

+125
-129
lines changed

‎src/frdata.py

Lines changed: 91 additions & 129 deletions
Original file line numberDiff line numberDiff line change
@@ -3,13 +3,12 @@
33
Frequency response data representation and functions.
44
55
This file contains the FRD class and also functions
6-
that operate on transfer functions. This is the primary representation
6+
that operate on FRD data. This is the primary representation
77
for the python-control library.
88
99
Routines in this module:
1010
1111
FRD.__init__
12-
FRD._truncatecoeff
1312
FRD.copy
1413
FRD.__str__
1514
FRD.__neg__
@@ -26,16 +25,13 @@
2625
FRD.pole
2726
FRD.zero
2827
FRD.feedback
29-
FRD.returnScipySignalLti
3028
FRD._common_den
31-
_tfpolyToString
32-
_addSISO
3329
_convertToFRD
3430
3531
"""
3632

3733
"""Copyright (c) 2010 by California Institute of Technology
38-
and 2012 Delft University of Technology
34+
Copyright (c) 2012 by Delft University of Technology
3935
All rights reserved.
4036
4137
Redistribution and use in source and binary forms, with or without
@@ -49,7 +45,8 @@
4945
notice, this list of conditions and the following disclaimer in the
5046
documentation and/or other materials provided with the distribution.
5147
52-
3. Neither the name of the California Institute of Technology nor
48+
3. Neither the names of the California Institute of Technology nor
49+
the Delft University of Technology nor
5350
the names of its contributors may be used to endorse or promote
5451
products derived from this software without specific prior
5552
written permission.
@@ -67,7 +64,7 @@
6764
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
6865
SUCH DAMAGE.
6966
70-
Author: M.M. (Rene) van Paassen
67+
Author: M.M. (Rene) van Paassen (using xferfcn.py as basis)
7168
Date: 02 Oct 12
7269
Revised:
7370
@@ -77,8 +74,8 @@
7774

7875
# External function declarations
7976
from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \
80-
polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze
81-
from scipy.signal import lti
77+
polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze, inner
78+
from scipy.interpolate import interp1d
8279
from copy import deepcopy
8380
from lti import Lti
8481
import statesp
@@ -87,7 +84,7 @@ class FRD(Lti):
8784

8885
"""The FRD class represents (measured?) frequency response
8986
TF instances and functions.
90-
87+
4
9188
The FRD class is derived from the Lti parent class. It is used
9289
throughout the python-control library to represent systems in frequency
9390
response data form.
@@ -99,11 +96,15 @@ class FRD(Lti):
9996
10097
>>> frdata[2][5] = numpy.array([1., 0.8-0.2j, 0.2-0.8j])
10198
102-
means that the frequency response from the 6th input to the
103-
3rd output at the frequencies defined in omega is set the array above.
99+
means that the frequency response from the 6th input to the 3rd
100+
output at the frequencies defined in omega is set to the array
101+
above, i.e. the rows represent the outputs and the columns
102+
represent the inputs.
104103
105104
"""
106105

106+
epsw = 1e-8
107+
107108
def __init__(self, *args):
108109
"""Construct a transfer function.
109110
@@ -116,21 +117,26 @@ def __init__(self, *args):
116117
To call the copy constructor, call FRD(sys), where sys is a
117118
FRD object.
118119
120+
To construct frequency response data for an existing Lti
121+
object, other than an FRD, call FRD(sys, omega)
122+
119123
"""
120124

121125
if len(args) == 2:
122126
if not isinstance(args[0], FRD) and isinstance(args[0], Lti):
123127
# not an FRD, but still a system, second argument should be
124128
# the frequency range
125129
otherlti = args[0]
126-
127-
self.omega = array(args[1].sort(), dtype=float)
130+
print args[1]
131+
self.omega = array(args[1], dtype=float)
132+
self.omega.sort()
133+
numfreq = len(self.omega)
128134

129135
# calculate frequency response at my points
130136
self.fresp = empty(
131137
(otherlti.outputs, otherlti.inputs, numfreq),
132138
dtype=complex)
133-
for k, w in enumerate(omega):
139+
for k, w in enumerate(self.omega):
134140
self.fresp[:, :, k] = otherlti.evalfr(w)
135141

136142
else:
@@ -158,6 +164,13 @@ def __init__(self, *args):
158164
else:
159165
raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
160166

167+
# create interpolation functions
168+
self.ifunc = empty((self.fresp.shape[1], self.fresp.shape[0]),
169+
dtype=interp1d)
170+
for i in range(self.fresp.shape[1]):
171+
for j in range(self.fresp.shape[0]):
172+
self.ifunc[i,j] = interp1d(self.omega, self.fresp[i, j, :])
173+
161174
Lti.__init__(self, self.fresp.shape[1], self.fresp.shape[0])
162175

163176
def __str__(self):
@@ -180,7 +193,7 @@ def __str__(self):
180193
def __neg__(self):
181194
"""Negate a transfer function."""
182195

183-
return FRD(self.omega, -self.fresp)
196+
return FRD(-self.fresp, self.omega)
184197

185198
def __add__(self, other):
186199
"""Add two LTI objects (parallel connection)."""
@@ -203,7 +216,7 @@ def __add__(self, other):
203216
raise ValueError("The first summand has %i output(s), but the \
204217
second has %i." % (self.outputs, other.outputs))
205218

206-
return FRD(other.omega, self.frd + other.frd)
219+
return FRD(self.fresp + other.fresp, other.omega)
207220

208221
def __radd__(self, other):
209222
"""Right add two LTI objects (parallel connection)."""
@@ -226,9 +239,9 @@ def __mul__(self, other):
226239
# Convert the second argument to a transfer function.
227240
if isinstance(other, (int, float, long, complex)):
228241
other = _convertToFRD(other, inputs=self.inputs,
229-
outputs=self.inputs)
242+
outputs=self.inputs, omega=self.omega)
230243
else:
231-
other = _convertToFRD(other)
244+
other = _convertToFRD(other, omega=self.omega)
232245

233246
# Check that the input-output sizes are consistent.
234247
if self.inputs != other.outputs:
@@ -268,9 +281,9 @@ def __div__(self, other):
268281

269282
if isinstance(other, (int, float, long, complex)):
270283
other = _convertToFRD(other, inputs=self.inputs,
271-
outputs=self.inputs)
284+
outputs=self.inputs, omega=self.omega)
272285
else:
273-
other = _convertToFRD(other)
286+
other = _convertToFRD(other, omega=self.omega)
274287

275288

276289
if (self.inputs > 1 or self.outputs > 1 or
@@ -288,9 +301,9 @@ def __rdiv__(self, other):
288301
"""Right divide two LTI objects."""
289302
if isinstance(other, (int, float, long, complex)):
290303
other = _convertToFRD(other, inputs=self.inputs,
291-
outputs=self.inputs)
304+
outputs=self.inputs, omega=self.omega)
292305
else:
293-
other = _convertToFRD(other)
306+
other = _convertToFRD(other, omega=self.omega)
294307

295308
if (self.inputs > 1 or self.outputs > 1 or
296309
other.inputs > 1 or other.outputs > 1):
@@ -322,8 +335,7 @@ def evalfr(self, omega):
322335

323336
for i in range(self.outputs):
324337
for j in range(self.inputs):
325-
out[i][j] = (polyval(self.num[i][j], omega * 1.j) /
326-
polyval(self.den[i][j], omega * 1.j))
338+
out[i,j] = self.ifunc[i,j](omega)
327339

328340
return out
329341

@@ -348,71 +360,43 @@ def freqresp(self, omega):
348360

349361
for i in range(self.outputs):
350362
for j in range(self.inputs):
351-
fresp = map(lambda w: (polyval(self.num[i][j], w * 1.j) /
352-
polyval(self.den[i][j], w * 1.j)), omega)
353-
fresp = array(fresp)
354-
363+
fresp = self.ifunc[i,j](omega)
355364
mag[i, j, :] = abs(fresp)
356365
phase[i, j, :] = angle(fresp)
357366

358367
return mag, phase, omega
359368

360369
def feedback(self, other, sign=-1):
361-
"""Feedback interconnection between two LTI objects."""
362-
363-
other = _convertToFRD(other)
364-
365-
if (self.inputs > 1 or self.outputs > 1 or
366-
other.inputs > 1 or other.outputs > 1):
367-
# TODO: MIMO feedback
368-
raise NotImplementedError("FRD.feedback is currently \
369-
only implemented for SISO functions.")
370-
371-
num1 = self.num[0][0]
372-
den1 = self.den[0][0]
373-
num2 = other.num[0][0]
374-
den2 = other.den[0][0]
375-
376-
num = polymul(num1, den2)
377-
den = polyadd(polymul(den2, den1), -sign * polymul(num2, num1))
378-
379-
return FRD(num, den)
380-
381-
# For MIMO or SISO systems, the analytic expression is
382-
# self / (1 - sign * other * self)
383-
# But this does not work correctly because the state size will be too
384-
# large.
385-
386-
def returnScipySignalLti(self):
387-
"""Return a list of a list of scipy.signal.lti objects.
388-
389-
For instance,
370+
"""Feedback interconnection between two FRD objects."""
390371

391-
>>> out = tfobject.returnScipySignalLti()
392-
>>> out[3][5]
393-
394-
is a signal.scipy.lti object corresponding to the transfer function from
395-
the 6th input to the 4th output.
396-
397-
"""
398-
399-
# Preallocate the output.
400-
out = [[[] for j in range(self.inputs)] for i in range(self.outputs)]
372+
other = _convertToFRD(other, omega=self.omega)
401373

402-
for i in range(self.outputs):
403-
for j in range(self.inputs):
404-
out[i][j] = lti(self.num[i][j], self.den[i][j])
405-
406-
return out
374+
if (self.outputs != other.inputs or
375+
self.inputs != other.outputs):
376+
raise ValueError(
377+
"FRD.feedback, inputs/outputs mismatch")
378+
fresp = empty(self.outputs, self.inputs, len(other.omega))
379+
# TODO: vectorize this
380+
# TODO: handle omega re-mapping
381+
for k, w in enumerate(other.omega):
382+
for i in range(0, self.inputs):
383+
for j in range(0, self.outputs):
384+
fresp[i, j, k] = \
385+
self.fresp[i, j, k] / \
386+
(1.0 + inner(self.fresp[i, :, k], other.fresp[:, i, k]))
387+
388+
return FRD(fresp, other.omega)
407389

408390

409-
def _convertToFRD(sys, omega):
410-
"""Convert a system to transfer function form (if needed).
391+
def _convertToFRD(sys, omega, inputs=1, outputs=1):
392+
"""Convert a system to frequency response data form (if needed).
411393
412-
If sys is already a transfer function, then it is returned. If sys is a
413-
state space object, then it is converted to a transfer function and
414-
returned. If sys is a scalar, then the number of inputs and outputs can be
415-
specified manually, as in:
394+
If sys is already an frd, and its frequency range matches or
395+
overlaps the range given in omega then it is returned. If sys is
396+
another Lti object or a transfer function, then it is converted to
397+
a frequency response data at the specified omega. If sys is a
398+
scalar, then the number of inputs and outputs can be specified
399+
manually, as in:
416400
417401
>>> sys = _convertToFRD(3.) # Assumes inputs = outputs = 1
418402
>>> sys = _convertToFRD(1., inputs=3, outputs=2)
@@ -424,7 +408,8 @@ def _convertToFRD(sys, omega):
424408

425409
if isinstance(sys, FRD):
426410

427-
if (abs(omega - sys.omega) < eps).all():
411+
omega.sort()
412+
if (abs(omega - sys.omega) < FRD.epsw).all():
428413
# frequencies match, and system was already frd; simply use
429414
return sys
430415

@@ -433,56 +418,33 @@ def _convertToFRD(sys, omega):
433418
omega = omega[omega <= max(sys.omega)]
434419
if not omega:
435420
raise ValueError("Frequency ranges of FRD do not overlap")
436-
return FRDsys
437-
438-
elif isinstance(sys, statesp.StateSpace):
439-
try:
440-
from slycot import tb04ad
441-
if len(kw):
442-
raise TypeError("If sys is a StateSpace, \
443-
_convertToFRD cannot take keywords.")
444-
445-
# Use Slycot to make the transformation
446-
# Make sure to convert system matrices to numpy arrays
447-
tfout = tb04ad(sys.states, sys.inputs, sys.outputs, array(sys.A),
448-
array(sys.B), array(sys.C), array(sys.D), tol1=0.0)
449-
450-
# Preallocate outputs.
451-
num = [[[] for j in range(sys.inputs)] for i in range(sys.outputs)]
452-
den = [[[] for j in range(sys.inputs)] for i in range(sys.outputs)]
453-
454-
for i in range(sys.outputs):
455-
for j in range(sys.inputs):
456-
num[i][j] = list(tfout[6][i, j, :])
457-
# Each transfer function matrix row has a common denominator.
458-
den[i][j] = list(tfout[5][i, :])
459-
# print num
460-
# print den
461-
except ImportError:
462-
# If slycot is not available, use signal.lti (SISO only)
463-
if (sys.inputs != 1 or sys.outputs != 1):
464-
raise TypeError("No support for MIMO without slycot")
465-
466-
lti_sys = lti(sys.A, sys.B, sys.C, sys.D)
467-
num = squeeze(lti_sys.num)
468-
den = squeeze(lti_sys.den)
469-
print num
470-
print den
421+
422+
# if there would be data beyond the extremes, add omega points at the
423+
# edges
424+
if omega[0] - sys.omega[0] > FRD.epsw:
425+
omega.insert(sys.omega[0], 0)
426+
if sys.omega[-1] - omega[-1] > FRD.epsw:
427+
omega.append(sys.omega[-1])
428+
print "Adjusting frequency range in FRD"
429+
430+
fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
431+
for k, w in enumerate(omega):
432+
fresp[:, :, k] = sys.evalfr(w)
433+
434+
return FRD(fresp, omega)
471435

472-
return FRD(num, den)
473-
elif isinstance(sys, (int, long, float, complex)):
474-
if "inputs" in kw:
475-
inputs = kw["inputs"]
476-
else:
477-
inputs = 1
478-
if "outputs" in kw:
479-
outputs = kw["outputs"]
480-
else:
481-
outputs = 1
436+
elif isinstance(sys, Lti):
437+
omega.sort()
438+
fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
439+
for k, w in enumerate(omega):
440+
fresp[:, :, k] = sys.evalfr(w)
441+
442+
return FRD(fresp, omega)
482443

483-
num = [[[sys] for j in range(inputs)] for i in range(outputs)]
484-
den = [[[1] for j in range(inputs)] for i in range(outputs)]
444+
elif isinstance(sys, (int, long, float, complex)):
485445

486-
return FRD(num, den)
446+
fresp = ones((outputs, inputs, len(omega)), dtype=float)*sys
447+
448+
return FRD(fresp, omega)
487449
else:
488450
raise TypeError("Can't convert given type to FRD system.")

0 commit comments

Comments
 (0)