[python-control-discuss] SF.net SVN: python-control:[33] trunk/src

[<mur...@us...>]

Revision: 33
          http://python-control.svn.sourceforge.net/python-control/?rev=33&view=rev
Author:   murrayrm
Date:     2010-11-26 21:59:57 +0000 (Fri, 26 Nov 2010)

Log Message:
-----------
* Added Nichols plot, contributed by Allan McInnes
* Allan also split out code to determine default frequency range

Modified Paths:
--------------
    trunk/src/freqplot.py
    trunk/src/matlab.py

Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py	2010-11-26 21:56:26 UTC (rev 32)
+++ trunk/src/freqplot.py	2010-11-26 21:59:57 UTC (rev 33)
@@ -4,7 +4,7 @@
 # Date: 24 May 09
 # 
 # This file contains some standard control system plots: Bode plots,
-# Nyquist plots and pole-zero diagrams
+# Nyquist plots, Nichols plots and pole-zero diagrams
 #
 # Copyright (c) 2010 by California Institute of Technology
 # All rights reserved.
@@ -46,6 +46,13 @@
 from ctrlutil import unwrap
 from bdalg import feedback
 
+#
+# Main plotting functions
+#
+# This section of the code contains the functions for generating
+# frequency domain plots
+#
+   
 # Bode plot
 def bode(syslist, omega=None, dB=False, Hz=False):
     """Bode plot for a system
@@ -81,36 +88,10 @@
     if (not getattr(syslist, '__iter__', False)):
         syslist = (syslist,)
 
-    #
     # Select a default range if none is provided
-    #
-    # This code looks at the poles and zeros of all of the systems that
-    # we are plotting and sets the frequency range to be one decade above
-    # and below the min and max feature frequencies, rounded to the nearest
-    # integer.  It excludes poles and zeros at the origin.  If no features
-    # are found, it turns logspace(-1, 1)
-    #
     if (omega == None):
-        # Find the list of all poles and zeros in the systems
-        features = np.array(())
-        for sys in syslist:
-            # Add new features to the list
-            features = np.concatenate((features, np.abs(sys.poles)))
-            features = np.concatenate((features, np.abs(sys.zeros)))
+        omega = default_frequency_range(syslist)
 
-        # Get rid of poles and zeros at the origin
-        features = features[features != 0];
-
-        # Make sure there is at least one point in the range
-        if (features.shape[0] == 0): features = [1];
-
-        # Take the log of the features
-        features = np.log10(features)
-
-        # Set the to be an order of magnitude beyond any features
-        omega = sp.logspace(np.floor(np.min(features))-1, 
-                            np.ceil(np.max(features))+1)
-
     for sys in syslist:
         # Get the magnitude and phase of the system
         mag, phase, omega = sys.freqresp(omega)
@@ -154,7 +135,7 @@
     return (211, 212)
 
 # Nyquist plot
-def nyquist(sys, omega=None):
+def nyquist(syslist, omega=None):
     """Nyquist plot for a system
 
     Usage
@@ -165,8 +146,8 @@
 
     Parameters
     ----------
-    sys : linsys
-        Linear input/output system
+    syslist : linsys
+        List of linear input/output systems (single system is OK)
     omega : freq_range
         Range of frequencies (list or bounds) in rad/sec
 
@@ -174,27 +155,82 @@
     -------------
     None
     """
-
+    # If argument was a singleton, turn it into a list
+    if (not getattr(syslist, '__iter__', False)):
+        syslist = (syslist,)
+        
     # Select a default range if none is provided
-    #! TODO: This needs to be made more intelligent
     if (omega == None):
-        omega = sp.logspace(-2, 2);
+        omega = default_frequency_range(syslist)
 
-    # Get the magnitude and phase of the system
-    mag, phase, omega = sys.freqresp(omega)
+    for sys in syslist:
+        # Get the magnitude and phase of the system
+        mag, phase, omega = sys.freqresp(omega)
+    
+        # Compute the primary curve
+        x = sp.multiply(mag, sp.cos(phase));
+        y = sp.multiply(mag, sp.sin(phase));
+    
+        # Plot the primary curve and mirror image
+        plt.plot(x, y, '-');
+        plt.plot(x, -y, '--');
 
-    # Compute the primary curve
-    x = sp.multiply(mag, sp.cos(phase));
-    y = sp.multiply(mag, sp.sin(phase));
+    # Mark the -1 point
+    plt.plot([-1], [0], 'r+')
 
-    # Plot the primary curve and mirror image
-    plt.plot(x, y, '-');
-    plt.plot(x, -y, '--');
+# Nichols plot
+# Contributed by Allan McInnes <All...@ca...>
+#! TODO: need unit test code
+def nichols(syslist, omega=None):
+    """Nichols plot for a system
 
-    # Mark the -1 point
-    plt.plot([-1], [0], '+k')
+    Usage
+    =====
+    magh = nichols(sys, omega=None)
 
+    Plots a Nichols plot for the system over a (optional) frequency range.
+
+    Parameters
+    ----------
+    syslist : linsys
+        List of linear input/output systems (single system is OK)
+    omega : freq_range
+        Range of frequencies (list or bounds) in rad/sec
+
+    Return values
+    -------------
+    None
+    """
+
+    # If argument was a singleton, turn it into a list
+    if (not getattr(syslist, '__iter__', False)):
+        syslist = (syslist,)
+
+    # Select a default range if none is provided
+    if (omega == None):
+        omega = default_frequency_range(syslist)
+
+    for sys in syslist:
+        # Get the magnitude and phase of the system
+        mag, phase, omega = sys.freqresp(omega)
+    
+        # Convert to Nichols-plot format (phase in degrees, 
+        # and magnitude in dB)
+        x = unwrap(sp.degrees(phase), 360)
+        y = 20*sp.log10(mag)
+    
+        # Generate the plot
+        plt.plot(x, y)
+        
+    plt.xlabel('Phase (deg)')
+    plt.ylabel('Magnitude (dB)')
+    plt.title('Nichols Plot')
+
+    # Mark the -180 point
+    plt.plot([-180], [0], 'r+')
+    
 # Gang of Four
+#! TODO: think about how (and whether) to handle lists of systems
 def gangof4(P, C, omega=None):
     """Plot the "Gang of 4" transfer functions for a system
 
@@ -220,7 +256,7 @@
     # Select a default range if none is provided
     #! TODO: This needs to be made more intelligent
     if (omega == None):
-        omega = sp.logspace(-2, 2);
+        omega = default_frequency_range((P,C))
 
     # Compute the senstivity functions
     L = P*C;
@@ -240,3 +276,60 @@
 
     mag, phase, omega = S.freqresp(omega);
     plt.subplot(224); plt.loglog(omega, mag);
+
+#
+# Utility functions
+#
+# This section of the code contains some utility functions for
+# generating frequency domain plots
+#
+   
+# Compute reasonable defaults for axes
+def default_frequency_range(syslist):
+    """Compute a reasonable default frequency range for frequency
+    domain plots.
+
+    Usage
+    =====
+    omega = default_frequency_range(syslist)
+
+    Finds a reasonable default frequency range by examining the features
+    (poles and zeros) of the systems in syslist.
+
+    Parameters
+    ----------
+    syslist : linsys
+        List of linear input/output systems (single system is OK)
+
+    Return values
+    -------------
+    omega : freq_range
+        Range of frequencies in rad/sec
+    """
+    # This code looks at the poles and zeros of all of the systems that
+    # we are plotting and sets the frequency range to be one decade above
+    # and below the min and max feature frequencies, rounded to the nearest
+    # integer.  It excludes poles and zeros at the origin.  If no features
+    # are found, it turns logspace(-1, 1)
+    
+    # Find the list of all poles and zeros in the systems
+    features = np.array(())
+    for sys in syslist:
+        # Add new features to the list
+        features = np.concatenate((features, np.abs(sys.poles)))
+        features = np.concatenate((features, np.abs(sys.zeros)))
+
+    # Get rid of poles and zeros at the origin
+    features = features[features != 0];
+
+    # Make sure there is at least one point in the range
+    if (features.shape[0] == 0): features = [1];
+
+    # Take the log of the features
+    features = np.log10(features)
+
+    # Set the range to be an order of magnitude beyond any features
+    omega = sp.logspace(np.floor(np.min(features))-1, 
+                        np.ceil(np.max(features))+1)   
+                        
+    return omega

Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py	2010-11-26 21:56:26 UTC (rev 32)
+++ trunk/src/matlab.py	2010-11-26 21:59:57 UTC (rev 33)
@@ -66,7 +66,7 @@
 
 # Import MATLAB-like functions that can be used as-is
 from ctrlutil import unwrap
-from freqplot import nyquist, gangof4
+from freqplot import nyquist, nichols, gangof4
 from bdalg import series, parallel, negate, feedback
 from pzmap import pzmap
 from statefbk import ctrb, obsv, place, lqr
@@ -155,7 +155,7 @@
    lti/bodemag    - Bode magnitude diagram only
    sigma          - singular value frequency plot
 *  nyquist        - Nyquist plot
-   nichols        - Nichols plot
+*  nichols        - Nichols plot
    margin         - gain and phase margins
    lti/allmargin  - all crossover frequencies and related gain/phase margins
    lti/freqresp   - frequency response over a frequency grid


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.