Merge pull request #62 from cwrowley/matlab

Move matlab compatibility routines to subpackage.

Reviewed all changes and confirmed that these should give the same functionality, modulo the changes @cwrowley noted in the PR.

control/__init__.py

@@ -45,28 +45,26 @@
 """
 
 # Import functions from within the control system library
-# Should probably only import the exact functions we use...
-from .bdalg import series, parallel, negate, feedback
-from .delay import pade
-from .dtime import sample_system
-from .freqplot import bode_plot, nyquist_plot, gangof4_plot
-from .freqplot import bode, nyquist, gangof4
-from .lti import issiso, timebase, timebaseEqual, isdtime, isctime
-from .margins import stability_margins, phase_crossover_frequencies
-from .mateqn import lyap, dlyap, care, dare
-from .modelsimp import hsvd, modred, balred, era, markov, minreal
-from .nichols import nichols_plot, nichols
-from .phaseplot import phase_plot, box_grid
-from .pzmap import pzmap
-from .rlocus import root_locus
-from .statefbk import place, lqr, ctrb, obsv, gram, acker
-from .statesp import StateSpace
-from .timeresp import forced_response, initial_response, step_response, \
-    impulse_response
-from .xferfcn import TransferFunction
+# Note: the functions we use are specified as __all__ variables in the modules
+from .bdalg import *
+from .delay import *
+from .dtime import *
+from .freqplot import *
+from .lti import *
+from .margins import *
+from .mateqn import *
+from .modelsimp import *
+from .nichols import *
+from .phaseplot import *
+from .pzmap import *
+from .rlocus import *
+from .statefbk import *
+from .statesp import *
+from .timeresp import *
+from .xferfcn import *
 from .ctrlutil import *
-from .frdata import FRD
-from .canonical import canonical_form, reachable_form
+from .frdata import *
+from .canonical import *
 
 # Exceptions
 from .exception import *
@@ -77,22 +75,6 @@
 except ImportError:
     __version__ = "dev"
 
-# Import some of the more common (and benign) MATLAB shortcuts
-# By default, don't import conflicting commands here
-#! TODO (RMM, 4 Nov 2012): remove MATLAB dependencies from __init__.py
-#!
-#! Eventually, all functionality should be in modules *other* than matlab.
-#! This will allow inclusion of the matlab module to set up a different set
-#! of defaults from the main package.  At that point, the matlab module will
-#! allow provide compatibility with MATLAB but no package functionality.
-#!
-from .matlab import ss, tf, ss2tf, tf2ss, drss
-from .matlab import pole, zero, evalfr, freqresp, dcgain
-from .matlab import nichols, rlocus, margin
-        # bode and nyquist come directly from freqplot.py
-from .matlab import step, impulse, initial, lsim
-from .matlab import ssdata, tfdata
-
 # The following is to use Numpy's testing framework
 # Tests go under directory tests/, benchmarks under directory benchmarks/
 from numpy.testing import Tester

control/bdalg.py

@@ -58,6 +58,8 @@
 from . import statesp as ss
 from . import frdata as frd
 
+__all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect']
+
 def series(sys1, sys2):
     """Return the series connection sys2 * sys1 for --> sys1 --> sys2 -->.
 

control/canonical.py

@@ -9,6 +9,7 @@
 from numpy import zeros, shape, poly
 from numpy.linalg import inv
 
+__all__ = ['canonical_form', 'reachable_form']
 
 def canonical_form(xsys, form='reachable'):
     """Convert a system into canonical form

control/delay.py

@@ -3,7 +3,7 @@
 #
 # Author: Sawyer Fuller
 # Date: 26 Aug 2010
-# 
+#
 # This file contains functions for implementing time delays (currently
 # only the pade() function).
 #
@@ -16,16 +16,16 @@
 #
 # 1. Redistributions of source code must retain the above copyright
 #    notice, this list of conditions and the following disclaimer.
-# 
+#
 # 2. Redistributions in binary form must reproduce the above copyright
 #    notice, this list of conditions and the following disclaimer in the
 #    documentation and/or other materials provided with the distribution.
-# 
+#
 # 3. Neither the name of the California Institute of Technology nor
 #    the names of its contributors may be used to endorse or promote
 #    products derived from this software without specific prior
 #    written permission.
-# 
+#
 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
@@ -38,30 +38,32 @@
 # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
 # OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 # SUCH DAMAGE.
-# 
+#
 # $Id$
 
 # Python 3 compatability (needs to go here)
 from __future__ import print_function
 
+__all__ = ['pade']
+
 def pade(T, n=1):
-    """ 
+    """
     Create a linear system that approximates a delay.
-    
+
     Return the numerator and denominator coefficients of the Pade approximation.
-    
+
     Parameters
     ----------
     T : number
         time delay
     n : integer
         order of approximation
-        
+
     Returns
-    -------       
+    -------
     num, den : array
         Polynomial coefficients of the delay model, in descending powers of s.
-    
+
     Notes
     -----
     Based on an algorithm in Golub and van Loan, "Matrix Computation" 3rd.
@@ -79,7 +81,7 @@ def pade(T, n=1):
         for k in range(1, n+1):
             c = T * c * (n - k + 1)/(2 * n - k + 1)/k
             num[n - k] = c * (-1)**k
-            den[n - k] = c 
+            den[n - k] = c
         num = [coeff/den[0] for coeff in num]
         den = [coeff/den[0] for coeff in den]
-    return num, den 
+    return num, den

control/dtime.py

@@ -47,6 +47,9 @@
 """
 
 from .lti import isctime
+from .statesp import StateSpace, _convertToStateSpace
+
+__all__ = ['sample_system', 'c2d']
 
 # Sample a continuous time system
 def sample_system(sysc, Ts, method='zoh', alpha=None):
@@ -85,3 +88,30 @@ def sample_system(sysc, Ts, method='zoh', alpha=None):
         raise ValueError("First argument must be continuous time system")
 
     return sysc.sample(Ts, method, alpha)
+
+
+def c2d(sysc, Ts, method='zoh'):
+    '''
+    Return a discrete-time system
+
+    Parameters
+    ----------
+    sysc: LTI (StateSpace or TransferFunction), continuous
+        System to be converted
+
+    Ts: number
+        Sample time for the conversion
+
+    method: string, optional
+        Method to be applied,
+        'zoh'        Zero-order hold on the inputs (default)
+        'foh'        First-order hold, currently not implemented
+        'impulse'    Impulse-invariant discretization, currently not implemented
+        'tustin'     Bilinear (Tustin) approximation, only SISO
+        'matched'    Matched pole-zero method, only SISO
+    '''
+    #  Call the sample_system() function to do the work
+    sysd = sample_system(sysc, Ts, method)
+    if isinstance(sysc, StateSpace) and not isinstance(sysd, StateSpace):
+        return _convertToStateSpace(sysd)
+    return sysd

control/frdata.py

@@ -1,35 +1,9 @@
 from __future__ import division
-"""frdata.py
-
+"""
 Frequency response data representation and functions.
 
-This file contains the FRD class and also functions that operate on
+This module contains the FRD class and also functions that operate on
 FRD data.
-
-Routines in this module:
-
-FRD.__init__
-FRD.copy
-FRD.__str__
-FRD.__neg__
-FRD.__add__
-FRD.__radd__
-FRD.__sub__
-FRD.__rsub__
-FRD.__mul__
-FRD.__rmul__
-FRD.__div__
-FRD.__rdiv__
-FRD.__truediv__
-FRD.__rtruediv__
-FRD.evalfr
-FRD.freqresp
-FRD.pole
-FRD.zero
-FRD.feedback
-FRD._common_den
-_convertToFRD
-
 """
 
 """Copyright (c) 2010 by California Institute of Technology
@@ -80,6 +54,8 @@
 from scipy.interpolate import splprep, splev
 from .lti import LTI
 
+__all__ = ['FRD', 'frd']
+
 class FRD(LTI):
     """A class for models defined by Frequency Response Data (FRD)
 
@@ -480,3 +456,39 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
 
     raise TypeError('''Can't convert given type "%s" to FRD system.''' %
                     sys.__class__)
+
+def frd(*args):
+    '''
+    Construct a Frequency Response Data model, or convert a system
+
+    frd models store the (measured) frequency response of a system.
+
+    This function can be called in different ways:
+
+    ``frd(response, freqs)``
+        Create an frd model with the given response data, in the form of
+        complex response vector, at matching frequency freqs [in rad/s]
+
+    ``frd(sys, freqs)``
+        Convert an LTI system into an frd model with data at frequencies
+        freqs.
+
+    Parameters
+    ----------
+    response: array_like, or list
+        complex vector with the system response
+    freq: array_lik or lis
+        vector with frequencies
+    sys: LTI (StateSpace or TransferFunction)
+        A linear system
+
+    Returns
+    -------
+    sys: FRD
+        New frequency response system
+
+    See Also
+    --------
+    ss, tf
+    '''
+    return FRD(*args)

control/freqplot.py

@@ -49,6 +49,9 @@
 from .bdalg import feedback
 from .lti import isdtime, timebaseEqual
 
+__all__ = ['bode_plot', 'nyquist_plot', 'gangof4_plot',
+           'bode', 'nyquist', 'gangof4']
+
 #
 # Main plotting functions
 #
@@ -74,7 +77,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
     Hz : boolean
         If True, plot frequency in Hz (omega must be provided in rad/sec)
     deg : boolean
-        If True, return phase in degrees (else radians)
+        If True, plot phase in degrees (else radians)
     Plot : boolean
         If True, plot magnitude and phase
     *args, **kwargs:
@@ -85,9 +88,9 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
     mag : array (list if len(syslist) > 1)
         magnitude
     phase : array (list if len(syslist) > 1)
-        phase
+        phase in radians
     omega : array (list if len(syslist) > 1)
-        frequency
+        frequency in rad/sec
 
     Notes
     -----
@@ -126,16 +129,15 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
                 omega = default_frequency_range(syslist)
 
             # Get the magnitude and phase of the system
+            omega = np.array(omega)
             mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega)
             mag = np.atleast_1d(np.squeeze(mag_tmp))
             phase = np.atleast_1d(np.squeeze(phase_tmp))
             phase = unwrap(phase)
             if Hz:
-                omega_plot = omega_sys/(2*sp.pi)
+                omega_plot = omega_sys / (2 * np.pi)
             else:
                 omega_plot = omega_sys
-            if dB: mag = 20*sp.log10(mag)
-            if deg: phase = phase * 180 / sp.pi
 
             mags.append(mag)
             phases.append(phase)
@@ -147,7 +149,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
                 # Magnitude plot
                 plt.subplot(211);
                 if dB:
-                    plt.semilogx(omega_plot, mag, *args, **kwargs)
+                    plt.semilogx(omega_plot, 20 * np.log10(mag), *args, **kwargs)
                 else:
                     plt.loglog(omega_plot, mag, *args, **kwargs)
                 plt.hold(True);
@@ -159,7 +161,11 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
 
                 # Phase plot
                 plt.subplot(212);
-                plt.semilogx(omega_plot, phase, *args, **kwargs)
+                if deg:
+                    phase_plot = phase * 180 / np.pi
+                else:
+                    phase_plot = phase
+                plt.semilogx(omega_plot, phase_plot, *args, **kwargs)
                 plt.hold(True);
 
                 # Add a grid to the plot + labeling