Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

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

Merged
merged 1 commit into from

4 participants

@WarrenWeckesser
Collaborator

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.

More details

The function _default_response_frequencies in ltisys.py, used by signal.bode, attempts to compute a good range of frequencies for showing the frequency response of an LTI system. This actually duplicates the functionality provided by the function findfreqs in filter_design.py.

Because _default_response_frequencies doesn't use the zeros or the imaginary parts of the poles of the system, it can generate bad frequency ranges.

For example, the follow either fail or miss important features of the frequency response:

>>> w, mag, phase = bode(([1],[1, 0, 100]))    # Raises exception
>>> w, mag, phase = bode(([1],[1, 1e-2, 100])) # Misses the resonance peak near 10
>>> w, mag, phase = bode(([0.01j, -0.01j], [2,2], [0.5])) # Misses the zero at 0.01.

Compare the plots created by semilogx(w, mag) to the corresponding plots with data generated by scipy.signal.freqs, e.g

>>> system = lti([1], [1, 1e-2, 100])
>>> fw, fh = freqs(system.num, system.den, worN=1000)
>>> semilogx(fw, 20*numpy.log10(abs(fh)))

So in this pull request, I have changed the bode function to use scipy.signal.freqs instead of _default_response_frequencies. This improves the behavior of bode and removes redundant functionality.

I also added a docstring to findfreqs in filter_design.py. This is the function used by freqs to find the range of frequencies when the frequencies are not specified. It takes into account the zeros and the poles (real and imaginary parts) to choose the frequencies.

In the cases where _default_response_frequencies worked, it generally produced a larger range of values than findfreqs does, and this often looked better because it went a little further into the asymptotic ranges. Whether findfreqs should be tweaked to generate a wider range is something that could be discussed here or in another pull request.

@WarrenWeckesser WarrenWeckesser 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.
aac6bd0
@bjornfor

You are right, _default_response_frequencies() has bugs. But I don't yet understand all of what findfreqs() does. So I'm currently trying to figure that out (and testing it and comparing with Matlab). And pondering whether adding clarifying comments to findfreqs() or fixing _default_response_frequencies() is the better solution.

Taking zeros into account in _default... is not difficult. It's when the poles/zeros are in the positive half-plane that I am unsure of what to do.

@WarrenWeckesser
Collaborator

We want a function that can automatically create a good range of frequencies for an LTI system, but we don't need two of them. findfreqs is already there and works pretty well, which is why this pull request drops _default_response_frequencies in favor of findfreqs. Rather than fix _default_response_frequencies, why not improve findfreqs? Note that improve might eventually mean completely replace the code. (But that would be done in a separate pull request.)

The code in findfreqs appears to use some numerical heuristics for finding a good frequency range. It also includes the flourish of rounding the ends of the interval to powers of 10. It is easy enough to read the code to know what it does, but unfortunately, the reasoning behind the heuristics is not explained in the code, so why it does it is a bit opaque. It also makes it tough to add tests. Updating findfreqs with well-documented code and tests would be a welcome contribution. (I don't think this pull request has to wait for such a change, though.)

@bjornfor

Good points. Ok, you get an ACK from me. Thanks for fixing bugs in bode() :-)

I'll add "document/refactor findfreqs" to my TODO. But I wouldn't mind if someone beat me to it :-)

@pv
Owner
pv commented

I guess if both of you agree, this can be merged?

@bjornfor

Yes, I agree. Let's merge this to fix bode and remove duplicated code.

@pv pv merged commit aac6bd0 into scipy:master
@pv
Owner
pv commented

Thanks, merged.

@hazelnusse

@pv @bjornfor @WarrenWeckesser I filed a bug on a previous version of ltisys.py:

http://projects.scipy.org/scipy/ticket/1863

It is still present with what was just merged into master. I believe it is

https://github.com/scipy/scipy/blob/master/scipy/signal/ltisys.py#L904

I haven't explored this completely, but I think the issue is when sys.num has dimension 2 rather than dimension 1. This happens, for example with the following bit of code:

import numpy as np
from scipy.signal import bode

# Basic double integrator system
A = np.array([[0, 1], [0, 0]])
B = np.array([[0], [1]])
C = np.array([[1, 0]])
D = np.array([[0]])

w, mag, phase = bode((A, B, C, D))
@hazelnusse

Unless I'm using bode incorrectly, the example I gave should probably be added as some sort of test to make sure the bode() function supports the (A,B,C,D) interface as advertised.

@ClemensFMN ClemensFMN referenced this pull request from a commit
Commit has since been removed from the repository and is no longer available.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Feb 12, 2013
  1. @WarrenWeckesser

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

    WarrenWeckesser authored
    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.
This page is out of date. Refresh to see the latest.
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
47 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
View
17 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()
Something went wrong with that request. Please try again.