Skip to content
Browse files

Merge pull request #432 from WarrenWeckesser/bode-freqs

BUG: signal: Fix the automatic frequency selection of scipy.signal.bode
  • Loading branch information...
2 parents e5b5de4 + aac6bd0 commit 608b418241931629c19642895e2176dce61c3bb0 @pv pv committed
Showing with 42 additions and 49 deletions.
  1. +28 −0 scipy/signal/filter_design.py
  2. +4 −42 scipy/signal/ltisys.py
  3. +10 −7 scipy/signal/tests/test_ltisys.py
View
28 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
View
46 scipy/signal/ltisys.py
@@ -613,44 +613,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 frequency 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.
@@ -936,14 +898,14 @@ 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
View
17 scipy/signal/tests/test_ltisys.py
@@ -285,19 +285,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)
@@ -308,6 +304,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)
+
class Test_freqresp(object):
@@ -377,5 +379,6 @@ def test_pole_zero(self):
w, H = freqresp(system, n=2)
assert_equal(w[0], 0.01) # a fail would give not-a-number
+
if __name__ == "__main__":
run_module_suite()

0 comments on commit 608b418

Please sign in to comment.
Something went wrong with that request. Please try again.