diff --git a/scipy/signal/cont2discrete.py b/scipy/signal/cont2discrete.py
index d17f8f2749db..1f95884d9b89 100644
--- a/scipy/signal/cont2discrete.py
+++ b/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)
\ No newline at end of file
diff --git a/scipy/signal/dltisys.py b/scipy/signal/dltisys.py
index 59d1530acc77..8af12aa897c6 100644
--- a/scipy/signal/dltisys.py
+++ b/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,6 +379,16 @@ 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
@@ -178,7 +396,7 @@ def dimpulse(system, x0=None, t=None, n=None):
     # 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
diff --git a/scipy/signal/tests/test_cont2discrete.py b/scipy/signal/tests/test_cont2discrete.py
index d34b241f8c09..97ad94b3c630 100644
--- a/scipy/signal/tests/test_cont2discrete.py
+++ b/scipy/signal/tests/test_cont2discrete.py
@@ -24,13 +24,14 @@ def test_zoh(self):
         # c and d in discrete should be equal to their continuous counterparts
         dt_requested = 0.5
 
-        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested, method='zoh')
+        # ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested, method='zoh')
+        sys = c2d((ac, bc, cc, dc), dt_requested, method='zoh')
 
-        assert_array_almost_equal(ad_truth, ad)
-        assert_array_almost_equal(bd_truth, bd)
-        assert_array_almost_equal(cc, cd)
-        assert_array_almost_equal(dc, dd)
-        assert_almost_equal(dt_requested, dt)
+        assert_array_almost_equal(ad_truth, sys.A)
+        assert_array_almost_equal(bd_truth, sys.B)
+        assert_array_almost_equal(cc, sys.C)
+        assert_array_almost_equal(dc, sys.D)
+        assert_almost_equal(dt_requested, sys.dt)
 
     def test_gbt(self):
         ac = np.eye(2)
@@ -50,13 +51,15 @@ def test_gbt(self):
                              [0.2],
                              [-0.205]])
 
-        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
-                                 method='gbt', alpha=alpha)
+        #ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
+        #                         method='gbt', alpha=alpha)
+        sys = c2d((ac, bc, cc, dc), dt_requested, method='gbt', alpha=alpha)
 
-        assert_array_almost_equal(ad_truth, ad)
-        assert_array_almost_equal(bd_truth, bd)
-        assert_array_almost_equal(cd_truth, cd)
-        assert_array_almost_equal(dd_truth, dd)
+
+        assert_array_almost_equal(ad_truth, sys.A)
+        assert_array_almost_equal(bd_truth, sys.B)
+        assert_array_almost_equal(cd_truth, sys.C)
+        assert_array_almost_equal(dd_truth, sys.D)
 
     def test_euler(self):
         ac = np.eye(2)
@@ -73,14 +76,13 @@ def test_euler(self):
                              [1.0, 0.25]])
         dd_truth = dc
 
-        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
-                                 method='euler')
+        sys = c2d((ac, bc, cc, dc), dt_requested, method='euler')
 
-        assert_array_almost_equal(ad_truth, ad)
-        assert_array_almost_equal(bd_truth, bd)
-        assert_array_almost_equal(cd_truth, cd)
-        assert_array_almost_equal(dd_truth, dd)
-        assert_almost_equal(dt_requested, dt)
+        assert_array_almost_equal(ad_truth, sys.A)
+        assert_array_almost_equal(bd_truth, sys.B)
+        assert_array_almost_equal(cd_truth, sys.C)
+        assert_array_almost_equal(dd_truth, sys.D)
+        assert_almost_equal(dt_requested, sys.dt)
 
     def test_backward_diff(self):
         ac = np.eye(2)
@@ -99,13 +101,12 @@ def test_backward_diff(self):
                              [1.0],
                              [0.295]])
 
-        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
-                                 method='backward_diff')
+        sys = c2d((ac, bc, cc, dc), dt_requested, method='backward_diff')
 
-        assert_array_almost_equal(ad_truth, ad)
-        assert_array_almost_equal(bd_truth, bd)
-        assert_array_almost_equal(cd_truth, cd)
-        assert_array_almost_equal(dd_truth, dd)
+        assert_array_almost_equal(ad_truth, sys.A)
+        assert_array_almost_equal(bd_truth, sys.B)
+        assert_array_almost_equal(cd_truth, sys.C)
+        assert_array_almost_equal(dd_truth, sys.D)
 
     def test_bilinear(self):
         ac = np.eye(2)
@@ -124,14 +125,13 @@ def test_bilinear(self):
                              [1.0 / 3.0],
                              [-0.121666666666667]])
 
-        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
-                                 method='bilinear')
+        sys = c2d((ac, bc, cc, dc), dt_requested, method='bilinear')
 
-        assert_array_almost_equal(ad_truth, ad)
-        assert_array_almost_equal(bd_truth, bd)
-        assert_array_almost_equal(cd_truth, cd)
-        assert_array_almost_equal(dd_truth, dd)
-        assert_almost_equal(dt_requested, dt)
+        assert_array_almost_equal(ad_truth, sys.A)
+        assert_array_almost_equal(bd_truth, sys.B)
+        assert_array_almost_equal(cd_truth, sys.C)
+        assert_array_almost_equal(dd_truth, sys.D)
+        assert_almost_equal(dt_requested, sys.dt)
 
         # Same continuous system again, but change sampling rate
 
@@ -142,14 +142,14 @@ def test_bilinear(self):
 
         dt_requested = 1.0 / 3.0
 
-        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
+        sys = c2d((ac, bc, cc, dc), dt_requested,
                                  method='bilinear')
 
-        assert_array_almost_equal(ad_truth, ad)
-        assert_array_almost_equal(bd_truth, bd)
-        assert_array_almost_equal(cd_truth, cd)
-        assert_array_almost_equal(dd_truth, dd)
-        assert_almost_equal(dt_requested, dt)
+        assert_array_almost_equal(ad_truth, sys.A)
+        assert_array_almost_equal(bd_truth, sys.B)
+        assert_array_almost_equal(cd_truth, sys.C)
+        assert_array_almost_equal(dd_truth, sys.D)
+        assert_almost_equal(dt_requested, sys.dt)
 
     def test_transferfunction(self):
         numc = np.array([0.25, 0.25, 0.5])
@@ -160,11 +160,11 @@ def test_transferfunction(self):
 
         dt_requested = 0.5
 
-        num, den, dt = c2d((numc, denc), dt_requested, method='zoh')
+        sys = c2d((numc, denc), dt_requested, method='zoh')
 
-        assert_array_almost_equal(numd, num)
-        assert_array_almost_equal(dend, den)
-        assert_almost_equal(dt_requested, dt)
+        assert_array_almost_equal(numd, sys.num)
+        assert_array_almost_equal(dend, sys.den)
+        assert_almost_equal(dt_requested, sys.dt)
 
     def test_zerospolesgain(self):
         zeros_c = np.array([0.5, -0.5])
@@ -178,13 +178,13 @@ def test_zerospolesgain(self):
 
         dt_requested = 0.5
 
-        zeros, poles, k, dt = c2d((zeros_c, poles_c, k_c), dt_requested,
+        sys = 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)
+        assert_array_almost_equal(zeros_d, sys.zeros)
+        assert_array_almost_equal(polls_d, sys.poles)
+        assert_almost_equal(k_d, sys.gain)
+        assert_almost_equal(dt_requested, sys.dt)
 
     def test_gbt_with_sio_tf_and_zpk(self):
         """Test method='gbt' with alpha=0.25 for tf and zpk cases."""
@@ -213,20 +213,21 @@ def test_gbt_with_sio_tf_and_zpk(self):
         dnum, dden = ss2tf(Ad, Bd, Cd, Dd)
 
         # Compute the discrete tf using cont2discrete.
-        c2dnum, c2dden, dt = c2d((cnum, cden), h, method='gbt', alpha=alpha)
+        sys = c2d((cnum, cden), h, method='gbt', alpha=alpha)
 
-        assert_allclose(dnum, c2dnum)
-        assert_allclose(dden, c2dden)
+        assert_allclose(dnum, sys.num)
+        assert_allclose(dden, sys.den)
 
         # Convert explicit solution to zpk.
         dz, dp, dk = ss2zpk(Ad, Bd, Cd, Dd)
 
         # Compute the discrete zpk using cont2discrete.
-        c2dz, c2dp, c2dk, dt = c2d((cz, cp, ck), h, method='gbt', alpha=alpha)
+        #c2dz, c2dp, c2dk, dt = c2d((cz, cp, ck), h, method='gbt', alpha=alpha)
+        sys = c2d((cz, cp, ck), h, method='gbt', alpha=alpha)
 
-        assert_allclose(dz, c2dz)
-        assert_allclose(dp, c2dp)
-        assert_allclose(dk, c2dk)
+        assert_allclose(dz, sys.zeros)
+        assert_allclose(dp, sys.poles)
+        assert_allclose(dk, sys.gain)
 
     def test_discrete_approx(self):
         """