Dltisys enhancement

scipy/signal/cont2discrete.py

@@ -9,7 +9,8 @@
 import numpy as np
 from scipy import linalg
 
-from .ltisys import tf2ss, ss2tf, zpk2ss, ss2zpk
+from .ltisys import tf2ss, ss2tf, zpk2ss, ss2zpk, lti
+from .dltisys import dlti
 
 __all__ = ['cont2discrete']
 
@@ -48,6 +49,7 @@ def cont2discrete(sys, dt, method="zoh", alpha=None):
     sysd : tuple containing the discrete system
         Based on the input type, the output will be of the form
 
+        * (lti instance,dt)   for lti class object input
         * (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
@@ -74,20 +76,13 @@ def cont2discrete(sys, dt, method="zoh", alpha=None):
         (http://www.ece.ualberta.ca/~gfzhang/research/ZCC07_preprint.pdf)
 
     """
-    if len(sys) == 2:
-        sysd = cont2discrete(tf2ss(sys[0], sys[1]), dt, method=method,
-                             alpha=alpha)
-        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,
-                             alpha=alpha)
-        return ss2zpk(sysd[0], sysd[1], sysd[2], sysd[3]) + (dt,)
-    elif len(sys) == 4:
-        a, b, c, d = sys
-    else:
-        raise ValueError("First argument must either be a tuple of 2 (tf), "
-                         "3 (zpk), or 4 (ss) arrays.")
 
+
+    if isinstance(sys, lti):
+        sys = sys
+    else:
+        sys = lti(*sys)
+    a,b,c,d = sys.A,sys.B,sys.C,sys.D
     if method == 'gbt':
         if alpha is None:
             raise ValueError("Alpha parameter must be specified for the "
@@ -139,4 +134,4 @@ def cont2discrete(sys, dt, method="zoh", alpha=None):
     else:
         raise ValueError("Unknown transformation method '%s'" % method)
 
-    return ad, bd, cd, dd, dt
+    return dlti(ad, bd, cd, dd, dt)

scipy/signal/dltisys.py

@@ -6,13 +6,210 @@
 # April 4, 2011
 from __future__ import division, print_function, absolute_import
 
+from .filter_design import tf2zpk, zpk2tf, normalize, freqs
 import numpy as np
 from scipy.interpolate import interp1d
-from .ltisys import tf2ss, zpk2ss
+from .ltisys import tf2ss, ss2tf, zpk2ss, ss2zpk, abcd_normalize
 
-__all__ = ['dlsim', 'dstep', 'dimpulse']
+__all__ = ['dlsim', 'dstep', 'dimpulse', 'dlti']
 
 
+class dlti(object):
+    """Discrete Linear Time Invariant class which simplifies representation.
+
+    Parameters
+    ----------
+    args : arguments
+        The `dlti` class can be instantiated with either 2, 3 or 4 arguments.
+        The following gives the number of elements in the tuple and the
+        interpretation:
+
+            * 2: (numerator, denominator)
+            * 3: (zeros, poles, gain)
+            * 4: (A, B, C, D)
+
+        Each argument can be an array or sequence.
+
+    Notes
+    -----
+    `dlti` instances have all types of representations available; for example
+    after creating an instance z with ``(zeros, poles, gain)`` the transfer
+    function representation (numerator, denominator) can be accessed as
+    ``z.num`` and ``z.den``.
+
+    """
+
+    def __init__(self, *args, **kwords):
+        """
+        Initialize the DLTI system using either:
+
+            - (numerator, denominator, dt)
+            - (zeros, poles, gain, dt)
+            - (A, B, C, D, dt) : state-space.
+
+        """
+
+        N = len(args)
+        if N == 3:  # Numerator denominator dt transfer function input
+            self._num, self._den, = normalize(args[0],args[1])
+            self._dt = args[2]
+            self._update(N)
+            self.inputs = 1
+            if len(self.num.shape) > 1:
+                self.outputs = self.num.shape[0]
+            else:
+                self.outputs = 1
+        elif N == 4:      # Zero-pole-gain dt form
+            self._zeros, self._poles, self._gain, self._dt = args
+            self._update(N)
+            # make sure we have numpy arrays
+            self.zeros = np.asarray(self.zeros)
+            self.poles = np.asarray(self.poles)
+            self.inputs = 1
+            if len(self.zeros.shape) > 1:
+                self.outputs = self.zeros.shape[0]
+            else:
+                self.outputs = 1
+        elif N == 5:       # State-space with dt form
+            self._A, self._B, self._C, self._D, = abcd_normalize(args[0],args[1],args[2],args[3])
+            self._dt = args[4]
+            self._update(N)
+            self.inputs = self.B.shape[-1]
+            self.outputs = self.C.shape[0]
+        else:
+            raise ValueError("Needs 3, 4, or 5 arguments.")
+
+    def __repr__(self):
+        """
+        Canonical representation using state-space to preserve numerical
+        precision and any MIMO information
+        """
+        return '{0}(\n{1},\n{2},\n{3},\n{4}\n,{5})'.format(
+            self.__class__.__name__,
+            repr(self.A),
+            repr(self.B),
+            repr(self.C),
+            repr(self.D),
+            repr(self.dt),
+            )
+
+    @property
+    def num(self):
+        return self._num
+
+    @num.setter
+    def num(self, value):
+        self._num = value
+        self._update(3)
+
+    @property
+    def den(self):
+        return self._den
+
+    @den.setter
+    def den(self, value):
+        self._den = value
+        self._update(3)
+
+    @property
+    def zeros(self):
+        return self._zeros
+
+    @zeros.setter
+    def zeros(self, value):
+        self._zeros = value
+        self._update(4)
+
+    @property
+    def poles(self):
+        return self._poles
+
+    @poles.setter
+    def poles(self, value):
+        self._poles = value
+        self._update(4)
+
+    @property
+    def gain(self):
+        return self._gain
+
+    @gain.setter
+    def gain(self, value):
+        self._gain = value
+        self._update(4)
+
+    @property
+    def A(self):
+        return self._A
+
+    @A.setter
+    def A(self, value):
+        self._A = value
+        self._update(5)
+
+    @property
+    def B(self):
+        return self._B
+
+    @B.setter
+    def B(self, value):
+        self._B = value
+        self._update(5)
+
+    @property
+    def C(self):
+        return self._C
+
+    @C.setter
+    def C(self, value):
+        self._C = value
+        self._update(5)
+
+    @property
+    def D(self):
+        return self._D
+
+    @D.setter
+    def D(self, value):
+        self._D = value
+        self._update(5)
+
+    @property
+    def dt(self):
+        return self._dt
+
+    @dt.setter
+    def dt(self,value):
+        self._dt = value
+
+    def _update(self, N):
+        if N == 3:
+            self._zeros, self._poles, self._gain = tf2zpk(self.num, self.den)
+            self._A, self._B, self._C, self._D = tf2ss(self.num, self.den)
+        if N == 4:
+            self._num, self._den = zpk2tf(self.zeros, self.poles, self.gain)
+            self._A, self._B, self._C, self._D = zpk2ss(self.zeros,
+                                                        self.poles, self.gain)
+        if N == 5:
+            self._num, self._den = ss2tf(self.A, self.B, self.C, self.D)
+            self._zeros, self._poles, self._gain = ss2zpk(self.A, self.B,
+                                                          self.C, self.D)
+
+
+    def dimpulse(self, x0=None, t=None, n=None):
+        """
+        Return the impulse response of a discrete-time system.
+
+        """
+        return dimpulse(self, x0=x0, t=t, n=n)
+
+    def dstep(self, x0=None, t=None, n=None):
+        """
+        Return the step response of a continuous-time system.
+
+        """
+        return dstep(self, x0=x0, t=t, n=n)
+
 def dlsim(system, u, t=None, x0=None):
     """
     Simulate output of a discrete-time linear system.
@@ -66,6 +263,8 @@ def dlsim(system, u, t=None, x0=None):
     >>> y
     array([ 0.,  0.,  0.,  1.])
 
+    """
+    ### OBRAS #####
     """
     if len(system) == 3:
         a, b, c, d = tf2ss(system[0], system[1])
@@ -79,6 +278,13 @@ def dlsim(system, u, t=None, x0=None):
         raise ValueError("System argument should be a discrete transfer " +
                          "function, zeros-poles-gain specification, or " +
                          "state-space system")
+    """
+    if isinstance(system, dlti):
+        sys = system
+    else:
+        sys = dlti(*system)
+
+    a, b, c, d, dt = sys.A, sys.B, sys.C, sys.D, sys.dt
 
     if t is None:
         out_samples = max(u.shape)
@@ -109,12 +315,12 @@ def dlsim(system, u, t=None, x0=None):
         u_dt = u_dt_interp(tout).transpose()
 
     # Simulate the system
-    for i in range(0, out_samples - 1):
-        xout[i+1,:] = np.dot(a, xout[i,:]) + np.dot(b, u_dt[i,:])
-        yout[i,:] = np.dot(c, xout[i,:]) + np.dot(d, u_dt[i,:])
+    for i in range(0, out_samples-1):
+        xout[i + 1,:] = np.dot(a, xout[i,:]) + np.dot(b, u_dt[i,:])
+        yout[i, :] = np.dot(c, xout[i,:]) + np.dot(d, u_dt[i,:])
 
     # Last point
-    yout[out_samples-1,:] = np.dot(c, xout[out_samples-1,:]) + \
+    yout[out_samples - 1,:] = np.dot(c, xout[out_samples-1,:]) + \
                             np.dot(d, u_dt[out_samples-1,:])
 
     if len(system) == 5:
@@ -155,6 +361,8 @@ def dimpulse(system, x0=None, t=None, n=None):
     --------
     impulse, dstep, dlsim, cont2discrete
 
+    """
+    # Obras ###
     """
     # Determine the system type and set number of inputs and time steps
     if len(system) == 3:
@@ -171,14 +379,24 @@ def dimpulse(system, x0=None, t=None, n=None):
                          "function, zeros-poles-gain specification, or " +
                          "state-space system")
 
+    """
+
+    if isinstance(system, dlti):
+        system = system
+    else:
+        system = dlti(*system)
+
+    n_inputs = system.A.shape[1]
+
+
     # Default to 100 samples if unspecified
     if n is None:
         n = 100
 
     # If time is not specified, use the number of samples
     # and system dt
     if t is None:
-        t = np.linspace(0, n * dt, n, endpoint=False)
+        t = np.linspace(0, n * system.dt, n, endpoint=False)
 
     # For each input, implement a step change
     yout = None