Skip to content

Commit de78c03

Browse files
repaghrepagh
committed
looks to be succesfull
1 parent d42454c commit de78c03

2 files changed

Lines changed: 74 additions & 55 deletions

File tree

‎src/frdata.py‎

Lines changed: 41 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -75,7 +75,7 @@
7575
# External function declarations
7676
from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \
7777
polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze, inner, \
78-
real, imag, matrix, absolute, eye, linalg, pi
78+
real, imag, matrix, absolute, eye, linalg, pi, where
7979
from scipy.interpolate import splprep, splev
8080
from copy import deepcopy
8181
from lti import Lti
@@ -106,7 +106,7 @@ class FRD(Lti):
106106

107107
epsw = 1e-8
108108

109-
def __init__(self, *args):
109+
def __init__(self, *args, **kwargs):
110110
"""Construct a transfer function.
111111
112112
The default constructor is FRD(w, d), where w is an iterable of
@@ -122,6 +122,7 @@ def __init__(self, *args):
122122
object, other than an FRD, call FRD(sys, omega)
123123
124124
"""
125+
smooth = kwargs.get('smooth', False)
125126

126127
if len(args) == 2:
127128
if not isinstance(args[0], FRD) and isinstance(args[0], Lti):
@@ -165,15 +166,17 @@ def __init__(self, *args):
165166
raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
166167

167168
# create interpolation functions
168-
self.ifunc = empty((self.fresp.shape[1], self.fresp.shape[0]),
169-
dtype=tuple)
170-
for i in range(self.fresp.shape[1]):
171-
for j in range(self.fresp.shape[0]):
172-
self.ifunc[i,j],u = splprep(
173-
u=self.omega, x=[real(self.fresp[i, j, :]),
174-
imag(self.fresp[i, j, :])],
175-
w=1.0/(absolute(self.fresp[i, j, :])+0.001), s=0.0)
176-
169+
if smooth:
170+
self.ifunc = empty((self.fresp.shape[1], self.fresp.shape[0]),
171+
dtype=tuple)
172+
for i in range(self.fresp.shape[1]):
173+
for j in range(self.fresp.shape[0]):
174+
self.ifunc[i,j],u = splprep(
175+
u=self.omega, x=[real(self.fresp[i, j, :]),
176+
imag(self.fresp[i, j, :])],
177+
w=1.0/(absolute(self.fresp[i, j, :])+0.001), s=0.0)
178+
else:
179+
self.ifunc = None
177180
Lti.__init__(self, self.fresp.shape[1], self.fresp.shape[0])
178181

179182
def __str__(self):
@@ -187,14 +190,11 @@ def __str__(self):
187190
for j in range(self.outputs):
188191
if mimo:
189192
outstr.append("Input %i to output %i:" % (i + 1, j + 1))
190-
outstr.append('Freq [rad/s] Magnitude Phase')
191-
outstr.append('------------ ----------- -----------')
192-
# outstr.extend(
193-
# [ '%12.3f %11.3e %11.2f' % (w, m, p*180.0/pi)
194-
# for m, p, w in zip(mt[i][j], pt[i][j], wt) ])
193+
outstr.append('Freq [rad/s] Response ')
194+
outstr.append('------------ ------------------------')
195195
outstr.extend(
196196
[ '%12.3f %10.4g + %10.4g' % (w, m, p)
197-
for m, p, w in zip(real(self.fresp[i,j,:]), imag(self.fresp[i,j,:]), wt) ])
197+
for m, p, w in zip(real(self.fresp[j,i,:]), imag(self.fresp[j,i,:]), wt) ])
198198

199199

200200
return '\n'.join(outstr)
@@ -316,10 +316,12 @@ def __rdiv__(self, other):
316316

317317
if (self.inputs > 1 or self.outputs > 1 or
318318
other.inputs > 1 or other.outputs > 1):
319-
raise NotImplementedError("FRD.__rdiv__ is currently \
320-
implemented only for SISO systems.")
319+
raise NotImplementedError(
320+
"FRD.__rdiv__ is currently implemented only for SISO systems.")
321321

322322
return other / self
323+
324+
323325
def __pow__(self,other):
324326
if not type(other) == int:
325327
raise ValueError("Exponent must be an integer")
@@ -342,10 +344,18 @@ def evalfr(self, omega):
342344
# Preallocate the output.
343345
out = empty((self.outputs, self.inputs), dtype=complex)
344346

345-
for i in range(self.outputs):
346-
for j in range(self.inputs):
347-
frraw = splev(omega, self.ifunc[i,j], der=0)
348-
out[i,j] = frraw[0] + 1.0j*frraw[1]
347+
if self.ifunc is None:
348+
try:
349+
out = self.fresp[:, :, where(self.omega == omega)[0][0]]
350+
except:
351+
raise ValueError(
352+
"Frequency %f not in frequency list, try an interpolating"
353+
" FRD if you want additional points")
354+
else:
355+
for i in range(self.outputs):
356+
for j in range(self.inputs):
357+
frraw = splev(omega, self.ifunc[i,j], der=0)
358+
out[i,j] = frraw[0] + 1.0j*frraw[1]
349359

350360
return out
351361

@@ -389,19 +399,13 @@ def feedback(self, other, sign=-1):
389399
# TODO: vectorize this
390400
# TODO: handle omega re-mapping
391401
for k, w in enumerate(other.omega):
392-
fresp[:, :, k] = linalg.solve(
402+
fresp[:, :, k] = self.fresp[:, :, k].view(type=matrix)* \
403+
linalg.solve(
393404
eye(self.inputs) +
394-
self.fresp[:, :, k].view(type=matrix) *
395-
other.fresp[:, :, k].view(type=matrix),
396-
eye(self.inputs))*self.fresp[:, :, k].view(type=matrix)
405+
other.fresp[:, :, k].view(type=matrix) *
406+
self.fresp[:, :, k].view(type=matrix),
407+
eye(self.inputs))
397408

398-
# for i in range(self.inputs):
399-
# for j in range(self.outputs):
400-
# fresp[i, j, k] = \
401-
# self.fresp[i, j, k] / \
402-
# (1.0-sign*inner(self.fresp[:, j, k],
403-
# other.fresp[i, :, k]))
404-
405409
return FRD(fresp, other.omega)
406410

407411
def _convertToFRD(sys, omega, inputs=1, outputs=1):
@@ -429,26 +433,9 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
429433
# frequencies match, and system was already frd; simply use
430434
return sys
431435

432-
# omega becomes lowest common range
433-
omega = omega[omega >= min(sys.omega)]
434-
omega = omega[omega <= max(sys.omega)]
435-
if not omega:
436-
raise ValueError("Frequency ranges of FRD do not overlap")
437-
438-
# if there would be data beyond the extremes, add omega points at the
439-
# edges
440-
if omega[0] - sys.omega[0] > FRD.epsw:
441-
omega.insert(sys.omega[0], 0)
442-
if sys.omega[-1] - omega[-1] > FRD.epsw:
443-
omega.append(sys.omega[-1])
444-
print "Adjusting frequency range in FRD"
445-
446-
fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
447-
for k, w in enumerate(omega):
448-
fresp[:, :, k] = sys.evalfr(w)
449-
450-
return FRD(fresp, omega)
451-
436+
raise NotImplementedError(
437+
"Frequency ranges of FRD do not match, conversion not implemented")
438+
452439
elif isinstance(sys, Lti):
453440
omega.sort()
454441
fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)

‎tests/frd_test.py‎

Lines changed: 33 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -101,7 +101,7 @@ def testAuto(self):
101101
def testNyquist(self):
102102
h1 = TransferFunction([1], [1, 2, 2])
103103
omega = np.logspace(-1, 2, 40)
104-
f1 = FRD(h1, omega)
104+
f1 = FRD(h1, omega, smooth=True)
105105
control.freqplot.nyquist(f1, np.logspace(-1, 2, 100))
106106
plt.savefig('/dev/null', format='svg')
107107
plt.figure(2)
@@ -136,7 +136,39 @@ def testMIMOfb(self):
136136
np.testing.assert_array_almost_equal(
137137
f1.freqresp([0.1, 1.0, 10])[1],
138138
f2.freqresp([0.1, 1.0, 10])[1])
139+
140+
def testMIMOfb2(self):
141+
sys = StateSpace(np.matrix('-2.0 0 0; 0 -1 1; 0 0 -3'),
142+
np.matrix('1.0 0; 0 0; 0 1'),
143+
np.eye(3), np.zeros((3,2)))
144+
omega = np.logspace(-1, 2, 10)
145+
K = np.matrix('1 0.3 0; 0.1 0 0')
146+
f1 = FRD(sys, omega).feedback(K)
147+
f2 = FRD(sys.feedback(K), omega)
148+
np.testing.assert_array_almost_equal(
149+
f1.freqresp([0.1, 1.0, 10])[0],
150+
f2.freqresp([0.1, 1.0, 10])[0])
151+
np.testing.assert_array_almost_equal(
152+
f1.freqresp([0.1, 1.0, 10])[1],
153+
f2.freqresp([0.1, 1.0, 10])[1])
154+
155+
def testAgainstOctave(self):
156+
157+
# with data from octave:
158+
#sys = ss([-2 0 0; 0 -1 1; 0 0 -3], [1 0; 0 0; 0 1], eye(3), zeros(3,2))
159+
#bfr = frd(bsys, [1])
160+
sys = StateSpace(np.matrix('-2.0 0 0; 0 -1 1; 0 0 -3'),
161+
np.matrix('1.0 0; 0 0; 0 1'),
162+
np.eye(3), np.zeros((3,2)))
163+
omega = np.logspace(-1, 2, 10)
164+
f1 = FRD(sys, omega)
165+
np.testing.assert_array_almost_equal(
166+
(f1.freqresp([1.0])[0] *
167+
np.exp(1j*f1.freqresp([1.0])[1])).reshape(3,2),
168+
np.matrix('0.4-0.2j 0; 0 0.1-0.2j; 0 0.3-0.1j'))
169+
139170
if __name__ == "__main__":
171+
140172
unittest.main()
141173
sys.exit(0)
142174

0 commit comments

Comments
 (0)