BUG: Fix the automatic frequency selection of scipy.signal.bode.

The automatic selection of the frequencies could miss zeros or resonance peaks, so use the existing function scipy.signal.findfreqs instead of _default_response_frequencies.

Specific changes:
* Modifed the function scipy.signal.bode to use scipy.signal.freqs.
* Added a docstring to scipy.signal.findfreqs.

scipy/signal/filter_design.py

@@ -28,6 +28,34 @@ class BadCoefficients(UserWarning):
 
 
 def findfreqs(num, den, N):
+    """
+    Find an array of frequencies for computing the response of a filter.
+
+    Parameters
+    ----------
+    num, den : array_like, 1-D
+        The polynomial coefficients of the numerator and denominator of the
+        transfer function of the filter or LTI system.  The coefficients are
+        ordered from highest to lowest degree.
+    N : int
+        The length of the array to be computed.
+
+    Returns
+    -------
+    w : (N,) ndarray
+        A 1-D array of frequencies, logarithmically spaced.
+
+    Examples
+    --------
+    Find a set of nine frequencies that span the "interesting part" of the
+    frequency response for the filter with the transfer function
+        H(s) = s / (s^2 + 8s + 25)
+
+    >>> findfreqs([1, 0], [1, 8, 25], N=9)
+    array([  1.00000000e-02,   3.16227766e-02,   1.00000000e-01,
+             3.16227766e-01,   1.00000000e+00,   3.16227766e+00,
+             1.00000000e+01,   3.16227766e+01,   1.00000000e+02])
+    """
     ep = atleast_1d(roots(den)) + 0j
     tz = atleast_1d(roots(num)) + 0j
 

scipy/signal/ltisys.py

@@ -11,7 +11,7 @@
 #   Rewrote lsim2 and added impulse2.
 #
 
-from .filter_design import tf2zpk, zpk2tf, normalize
+from .filter_design import tf2zpk, zpk2tf, normalize, freqs
 import numpy
 from numpy import product, zeros, array, dot, transpose, ones, \
     nan_to_num, zeros_like, linspace
@@ -600,43 +600,6 @@ def _default_response_times(A, n):
     return t
 
 
-def _default_response_frequencies(A, n):
-    """Compute a reasonable set of frequency points for bode plot.
-
-    This function is used by `bode` to compute the frequency points (in rad/s)
-    when the `w` argument to the function is None.
-
-    Parameters
-    ----------
-    A : ndarray
-        The system matrix, which is square.
-    n : int
-        The number of time samples to generate.
-
-    Returns
-    -------
-    w : ndarray
-        The 1-D array of length `n` of frequency samples (in rad/s) at which
-        the response is to be computed.
-    """
-    vals = linalg.eigvals(A)
-    # Remove poles at 0 because they don't help us determine an interesting
-    # frequency range. (And if we pass a 0 to log10() below we will crash.)
-    poles = [pole for pole in vals if pole != 0]
-    # If there are no non-zero poles, just hardcode something.
-    if len(poles) == 0:
-        minpole = 1
-        maxpole = 1
-    else:
-        minpole = min(abs(real(poles)))
-        maxpole = max(abs(real(poles)))
-    # A reasonable frequency range is two orders of magnitude before the
-    # minimum pole (slowest) and two orders of magnitude after the maximum pole
-    # (fastest).
-    w = numpy.logspace(numpy.log10(minpole) - 2, numpy.log10(maxpole) + 2, n)
-    return w
-
-
 def impulse(system, X0=None, T=None, N=None):
     """Impulse response of continuous-time system.
 
@@ -921,12 +884,12 @@ def bode(system, w=None, n=100):
         sys = lti(*system)
 
     if w is None:
-        w = _default_response_frequencies(sys.A, n)
+        worN = n
     else:
-        w = numpy.asarray(w)
+        worN = w
+    w, y = freqs(sys.num, sys.den, worN=worN)
 
-    jw = w * 1j
-    y = numpy.polyval(sys.num, jw) / numpy.polyval(sys.den, jw)
     mag = 20.0 * numpy.log10(abs(y))
     phase = numpy.arctan2(y.imag, y.real) * 180.0 / numpy.pi
+
     return w, mag, phase

scipy/signal/tests/test_ltisys.py

@@ -284,19 +284,15 @@ def test_04(self):
         assert_almost_equal(phase, expected_phase)
 
     def test_05(self):
-        """Test that bode() finds a reasonable frequency range.
-
-        A reasonable frequency range is two orders of magnitude before the
-        minimum (slowest) pole and two orders of magnitude after the maximum
-        (fastest) pole.
-        """
+        """Test that bode() finds a reasonable frequency range."""
         # 1st order low-pass filter: H(s) = 1 / (s + 1)
         system = lti([1], [1, 1])
         vals = linalg.eigvals(system.A)
         minpole = min(abs(np.real(vals)))
         maxpole = max(abs(np.real(vals)))
         n = 10;
-        expected_w = np.logspace(np.log10(minpole) - 2, np.log10(maxpole) + 2, n)
+        # Expected range is from 0.01 to 10.
+        expected_w = np.logspace(-2, 1, n)
         w, mag, phase = bode(system, n=n)
         assert_almost_equal(w, expected_w)
 
@@ -307,5 +303,12 @@ def test_06(self):
         w, mag, phase = bode(system, n=2)
         assert_equal(w[0], 0.01)  # a fail would give not-a-number
 
+    def test_07(self):
+        """bode() should not fail on a system with pure imaginary poles."""
+        # The test passes if bode doesn't raise an exception.
+        system = lti([1], [1, 0, 100])
+        w, mag, phase = bode(system, n=2)
+
+
 if __name__ == "__main__":
     run_module_suite()