Added zeros-poles-gain support to discrete LTI functions.

scipy/signal/cont2discrete.py

@@ -9,7 +9,7 @@
 import numpy.linalg
 import scipy.linalg
 
-from ltisys import tf2ss, ss2tf
+from ltisys import tf2ss, ss2tf, zpk2ss, ss2zpk
 
 __all__ = ['cont2discrete']
 
@@ -26,24 +26,24 @@ def cont2discrete(sys, dt, method="zoh"):
     -----------
     sys : a tuple describing the system.
         The following gives the number of elements in the tuple and
-        the interpretation::
-
-          - 2: (num, den)    defining a transfer function
-          - 4: (A, B, C, D)  defining a state-space system
-
+        the interpretation:
+        * 2: (num, den)
+        * 3: (zeros, poles, gain)
+        * 4: (A, B, C, D)
+        
     dt : float
         The discretization time step.
     method : {"bilinear", "zoh"}
         Which method to use, bilinear or zero-order hold ("zoh", the default).
 
     Returns
     -------
-    sysd : tuple
-        containing the discrete system
-        Based on the input type, the output will be of the form::
-
-          (num, den, dt)   for transfer function input
-          (A, B, C, D, dt) for state-space system input
+    sysd : tuple containing the discrete system
+        Based on the input type, the output will be of the form
+        
+        (num, den, dt)   for transfer function input
+        (zeros, poles, gain, dt)   for zeros-poles-gain input
+        (A, B, C, D, dt) for state-space system input
 
     Notes
     -----
@@ -61,6 +61,9 @@ def cont2discrete(sys, dt, method="zoh"):
     if len(sys) == 2:
         sysd = cont2discrete(tf2ss(sys[0], sys[1]), dt, method=method)
         return ss2tf(sysd[0], sysd[1], sysd[2], sysd[3]) + (dt,)
+    elif len(sys) == 3:
+        sysd = cont2discrete(zpk2ss(sys[0], sys[1], sys[2]), dt, method=method)
+        return ss2zpk(sysd[0], sysd[1], sysd[2], sysd[3]) + (dt,)
     elif len(sys) == 4:
         a, b, c, d = sys
     else:

scipy/signal/dltisys.py

@@ -8,11 +8,10 @@
 import math
 import numpy as np
 from scipy.interpolate import interp1d
-from ltisys import tf2ss
+from ltisys import tf2ss, zpk2ss
 
 __all__ = ['dlsim', 'dstep', 'dimpulse']
 
-
 def dlsim(system, u, t=None, x0=None):
     """
     Simulate output of a discrete-time linear system.
@@ -22,11 +21,10 @@ def dlsim(system, u, t=None, x0=None):
     system : class instance or tuple
         An instance of the LTI class, or a tuple describing the system.
         The following gives the number of elements in the tuple and
-        the interpretation::
-
-          - 3: (num, den, dt)
-          - 5: (A, B, C, D, dt)
-
+        the interpretation:
+          * 3: (num, den, dt)
+          * 4: (zeros, poles, gain, dt)
+          * 5: (A, B, C, D, dt)
     u : array_like
         An input array describing the input at each time `t` (interpolation is
         assumed between given times).  If there are multiple inputs, then each
@@ -55,11 +53,15 @@ def dlsim(system, u, t=None, x0=None):
     if len(system) == 3:
         a, b, c, d = tf2ss(system[0], system[1])
         dt = system[2]
+    elif len(system) == 4:
+        a, b, c, d = zpk2ss(system[0], system[1], system[2])
+        dt = system[3]
     elif len(system) == 5:
         a, b, c, d, dt = system
     else:
-        raise ValueError("System argument should be a discrete transfer "
-                         "function or state-space system")
+        raise ValueError("System argument should be a discrete transfer " +
+                         "function, zeros-poles-gain specification, or " +
+                         "state-space system")
 
     if t is None:
         out_samples = max(u.shape)
@@ -94,20 +96,23 @@ def dlsim(system, u, t=None, x0=None):
     # Last point
     yout[out_samples-1,:] = np.dot(c, xout[out_samples-1,:]) + \
                             np.dot(d, u_dt[out_samples-1,:])
-
-    if len(system) == 3:
-        return tout, yout
-    elif len(system) == 5:
+
+    if len(system) == 5:
         return tout, yout, xout
+    else:
+        return tout, yout
 
 def dimpulse(system, x0=None, t=None, n=None):
     """Impulse response of discrete-time system.
 
     Parameters
     ----------
     system : tuple
-        If specified as a tuple, the system is described as
-        ``(num, den)``, ``(zero, pole, gain)``, or ``(A, B, C, D)``.
+        The following gives the number of elements in the tuple and
+        the interpretation.
+          * 3: (num, den, dt)
+          * 4: (zeros, poles, gain, dt)
+          * 5: (A, B, C, D, dt)
     x0 : array_like, optional
         Initial state-vector.  Defaults to zero.
     t : array_like, optional
@@ -132,10 +137,17 @@ def dimpulse(system, x0=None, t=None, n=None):
     if len(system) == 3:
         n_inputs = 1
         dt = system[2]
+    elif len(system) == 4:
+        n_inputs = 1
+        dt = system[3]
     elif len(system) == 5:
         n_inputs = system[1].shape[1]
         dt = system[4]
-
+    else:
+        raise ValueError("System argument should be a discrete transfer " +
+                         "function, zeros-poles-gain specification, or " +
+                         "state-space system")
+
     # Default to 100 samples if unspecified
     if n is None:
         n = 100
@@ -169,11 +181,10 @@ def dstep(system, x0=None, t=None, n=None):
     ----------
     system : a tuple describing the system.
         The following gives the number of elements in the tuple and
-        the interpretation::
-
-          - 3 (num, den, dt)
-          - 5 (A, B, C, D, dt)
-
+        the interpretation.
+          * 3: (num, den, dt)
+          * 4: (zeros, poles, gain, dt)
+          * 5: (A, B, C, D, dt)
     x0 : array_like, optional
         Initial state-vector (default is zero).
     t : array_like, optional
@@ -198,10 +209,17 @@ def dstep(system, x0=None, t=None, n=None):
     if len(system) == 3:
         n_inputs = 1
         dt = system[2]
+    elif len(system) == 4:
+        n_inputs = 1
+        dt = system[3]
     elif len(system) == 5:
         n_inputs = system[1].shape[1]
         dt = system[4]
-
+    else:
+        raise ValueError("System argument should be a discrete transfer " +
+                         "function, zeros-poles-gain specification, or " +
+                         "state-space system")
+
     # Default to 100 samples if unspecified
     if n is None:
         n = 100

scipy/signal/tests/test_cont2discrete.py

@@ -115,6 +115,27 @@ def test_transferfunction(self):
         assert_array_almost_equal(dend, den)
 
         assert_almost_equal(dt_requested, dt)
+
+    def test_zerospolesgain(self):
+
+        zeros_c = np.array([0.5, -0.5])
+        poles_c = np.array([1.j / np.sqrt(2), -1.j / np.sqrt(2)])
+        k_c = 1.0
+
+        zeros_d = [1.23371727305860, 0.735356894461267]
+        polls_d = [0.938148335039729 + 0.346233593780536j,
+                   0.938148335039729 - 0.346233593780536j]
+        k_d = 1.0
+
+        dt_requested = 0.5
+
+        zeros, poles, k, dt = c2d((zeros_c, poles_c, k_c), dt_requested, 
+                                  method='zoh')
+
+        assert_array_almost_equal(zeros_d, zeros)
+        assert_array_almost_equal(polls_d, poles)
+        assert_almost_equal(k_d, k)
+        assert_almost_equal(dt_requested, dt)
 
 if __name__ == "__main__":
     run_module_suite()

scipy/signal/tests/test_dltisys.py

@@ -5,7 +5,7 @@
 import numpy as np
 from numpy.testing import TestCase, run_module_suite, assert_equal, \
                           assert_array_almost_equal, assert_almost_equal
-from scipy.signal import dlsim, dstep, dimpulse
+from scipy.signal import dlsim, dstep, dimpulse, tf2zpk
 
 class TestDLTI(TestCase):
 
@@ -67,6 +67,17 @@ def test_dlsim(self):
 
         assert_array_almost_equal(yout,yout_truth)
         assert_array_almost_equal(t_in, tout)
+
+        # zeros-poles-gain representation
+        zd = np.array([0.5, -0.5])
+        pd = np.array([1.j / np.sqrt(2), -1.j / np.sqrt(2)])
+        k = 1.0
+        yout_truth = np.asmatrix([0.0, 1.0, 2.0, 2.25, 2.5]).transpose()
+
+        tout, yout = dlsim((zd, pd, k, 0.5), u[:,0], t_in)
+
+        assert_array_almost_equal(yout,yout_truth)
+        assert_array_almost_equal(t_in, tout)
 
     def test_dstep(self):
 
@@ -96,6 +107,19 @@ def test_dstep(self):
         for i in range(0, len(yout)):
             assert_equal(yout[i].shape[0], 10)
             assert_array_almost_equal(yout[i].flatten(), yout_step_truth[i])
+
+        # Check that the other two inputs (tf, zpk) will work as well
+        tfin = ([1.0,], [1.0, 1.0], 0.5)
+        yout_tfstep = np.asarray([0.0, 1.0, 0.0])
+        tout, yout = dstep(tfin, n=3)
+        assert_equal(len(yout), 1)
+        assert_array_almost_equal(yout[0].flatten(), yout_tfstep)
+
+        zpkin = tf2zpk(tfin[0], tfin[1]) + (0.5,)
+        tout, yout = dstep(zpkin, n=3)
+        assert_equal(len(yout), 1)
+        assert_array_almost_equal(yout[0].flatten(), yout_tfstep)
+
 
     def test_dimpulse(self):
 
@@ -124,3 +148,16 @@ def test_dimpulse(self):
         for i in range(0, len(yout)):
             assert_equal(yout[i].shape[0], 10)
             assert_array_almost_equal(yout[i].flatten(), yout_imp_truth[i])
+
+        # Check that the other two inputs (tf, zpk) will work as well
+        tfin = ([1.0,], [1.0, 1.0], 0.5)
+        yout_tfimpulse = np.asarray([0.0, 1.0, -1.0])
+        tout, yout = dimpulse(tfin, n=3)
+        assert_equal(len(yout), 1)
+        assert_array_almost_equal(yout[0].flatten(), yout_tfimpulse)
+
+        zpkin = tf2zpk(tfin[0], tfin[1]) + (0.5,)
+        tout, yout = dimpulse(tfin, n=3)
+        assert_equal(len(yout), 1)
+        assert_array_almost_equal(yout[0].flatten(), yout_tfimpulse)
+