Speed up freqresp for transfer functions by an order of magnitude

cwrowley · cwrowley · commit 242d2b43252a · 2014-08-13T09:43:13.000-04:00
When determining the frequency response, the list of frequencies
is now stored as a numpy array, instead of a list, and evaluated
all at once, instead of using `map`.  This is much faster.

The improved performance can be measured using a script
`bench/time_freqresp.py`.  On my machine, I get a speedup of
more than a factor of 20 (from 3.9 seconds to 0.18 seconds) for
a 10th-order transfer function.
diff --git a/control/bench/time_freqresp.py b/control/bench/time_freqresp.py
@@ -0,0 +1,14 @@
+from control import tf
+from control.matlab import rss
+from numpy import logspace
+from timeit import timeit
+
+nstates = 10
+sys = rss(nstates)
+sys_tf = tf(sys)
+w = logspace(-1,1,50)
+ntimes = 1000
+time_ss = timeit("sys.freqresp(w)", setup="from __main__ import sys, w", number=ntimes)
+time_tf = timeit("sys_tf.freqresp(w)", setup="from __main__ import sys_tf, w", number=ntimes)
+print("State-space model on %d states: %f" % (nstates, time_ss))
+print("Transfer-function model on %d states: %f" % (nstates, time_tf))
diff --git a/control/xferfcn.py b/control/xferfcn.py
@@ -242,18 +242,18 @@ def __init__(self, *args):
         self._truncatecoeff()
 
     def __call__(self, s):
-        """Evaluate the system's transfer function for a complex vairable
+        """Evaluate the system's transfer function for a complex variable
 
         For a SISO transfer function, returns the value of the
         transfer function.  For a MIMO transfer fuction, returns a
         matrix of values evaluated at complex variable s."""
 
-        if (self.inputs > 1 or self.outputs > 1):
-            # MIMO transfer function, return a matrix
-            return self.horner(s)
-        else:
-            # SISO transfer function, return a scalar
+        if self.issiso():
+            # return a scalar
             return self.horner(s)[0][0]
+        else:
+            # return a matrix
+            return self.horner(s)
 
     def _truncatecoeff(self):
         """Remove extraneous zero coefficients from num and den.
@@ -615,15 +615,13 @@ def freqresp(self, omega):
             if (max(omega) * dt > pi):
                 warn("evalfr: frequency evaluation above Nyquist frequency")
         else:
-            slist = map(lambda w: 1.j * w, omega)
+            slist = np.array([1j * w for w in omega])
 
         # Compute frequency response for each input/output pair
         for i in range(self.outputs):
             for j in range(self.inputs):
-                fresp = map(lambda s: (polyval(self.num[i][j], s) /
-                            polyval(self.den[i][j], s)), slist)
-                fresp = array(list(fresp))
-
+                fresp = (polyval(self.num[i][j], slist) /
+                         polyval(self.den[i][j], slist))
                 mag[i, j, :] = abs(fresp)
                 phase[i, j, :] = angle(fresp)
 
