Skip to content

BUG: Use zpk in IIR filter design to improve accuracy #3085

Merged
merged 31 commits into from Jan 31, 2014

10 participants

@endolith

Temporary fix for #2443, using private functions until a decision is made about how to arrange the functions for backwards compatibility, etc.

Examples:
17-pole digital Butterworth bandpass
Before:
bad butter 17

After:
digital butter17

21-pole analog Butterworth lowpass at Wn=0.01
Before:
signal butter 21 01 pz

After:
mine butter 21 01 pz

cheby2 bandpass
Before:
sloppy cheby2 freq resp

After:
better cheby2 freq resp

endolith added some commits Nov 14, 2013
@endolith endolith ENH: Add 'lp' and 'hp' to band_dict
(since 'bp' is there and functions like 'lp2hp' imply it)
636be9b
@endolith endolith BUG: Use zpk in iirfilter to improve accuracy of higher-order filters
Previously iirfilter converted to tf representation internally, which is not numerically stable for high filter orders or cutoff frequencies close to 0.  Now it uses zpk format internally to improve accuracy.  Currently this uses private functions like _zpklp2lp, but they should eventually be merged into the existing lp2lp etc functions or made into their own functions, when a decision is made on the best way to do that.
99495fd
@endolith endolith BUG: Make all filter design functions return empty arrays for N=0 b2aa105
@coveralls

Coverage Status

Coverage remained the same when pulling b2aa105 on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@endolith endolith ENH: Clearer error message for incorrect Wn size
Previously it said "IndexError: 0-d arrays can't be indexed"
Combine bandpass and bandstop calculations for DRY
Also add quotes to error messages
557695b
@coveralls

Coverage Status

Coverage remained the same when pulling 557695b on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@Eric89GXL

Ahh @endolith I didn't see this PR, ignore my questions in a (much) older issue.

@Eric89GXL Eric89GXL referenced this pull request in mne-tools/mne-python Nov 21, 2013
Closed

BUG/API with filtering #912

@Eric89GXL Eric89GXL and 1 other commented on an outdated diff Nov 21, 2013
scipy/signal/filter_design.py
@@ -5,9 +5,10 @@
import warnings
import numpy
-from numpy import atleast_1d, poly, polyval, roots, real, asarray, allclose, \
- resize, pi, absolute, logspace, r_, sqrt, tan, log10, arctan, arcsinh, \
- cos, exp, cosh, arccosh, ceil, conjugate, zeros, sinh
+from numpy import (atleast_1d, poly, polyval, roots, real, asarray, allclose,
+ resize, pi, absolute, logspace, r_, sqrt, tan, log10, arctan, arcsinh,
@Eric89GXL
Eric89GXL added a note Nov 21, 2013

With parens you might as well indent to the open paren, PEP8.

@endolith
endolith added a note Nov 22, 2013

Ok, but is that really required? The example in http://www.python.org/dev/peps/pep-0328/ is this way:

from Tkinter import (Tk, Frame, Button, Entry, Canvas, Text,
    LEFT, DISABLED, NORMAL, RIDGE, END)

but http://docs.python.org/release/2.4/whatsnew/node10.html says

from SimpleXMLRPCServer import (SimpleXMLRPCServer,
                                SimpleXMLRPCRequestHandler,
                                CGIXMLRPCRequestHandler,
                                resolve_dotted_attribute)

so maybe it is?

@endolith
endolith added a note Nov 22, 2013

I guess it looks ok since numpy isn't long:

from numpy import (atleast_1d, poly, polyval, roots, real, asarray, allclose,
                   resize, pi, absolute, logspace, r_, sqrt, tan, log10, 
                   arctan, arcsinh, cos, exp, cosh, arccosh, ceil, conjugate, 
                   zeros, sinh, append, concatenate, prod, ones)

Should these be alphabetized or anything?

@endolith
endolith added a note Nov 22, 2013

also I thought it was discouraged to use r_ except as a shortcut while working on the command line

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@Eric89GXL Eric89GXL and 3 others commented on an outdated diff Nov 21, 2013
scipy/signal/filter_design.py
+
+# TODO: merge these into existing functions or make public versions
+
+def _zpkbilinear(z, p, k, fs):
+ """
+ Transform a set of poles and zeros from the analog s-plane to the digital
+ z-plane using Tustin's method, which maintains the shape of the
+ frequency response.
+ """
+ degree = len(p) - len(z)
+ if degree < 0:
+ raise ValueError("Improper transfer function. "
+ "Must have at least as many poles as zeros.")
+
+ # Compensate for gain change
+ k *= real(prod(2*fs - z) / prod(2*fs - p))
@Eric89GXL
Eric89GXL added a note Nov 21, 2013

This (and other places) also violates PEP8, need to do 2 * fs instead.

@endolith
endolith added a note Nov 21, 2013

"If operators with different priorities are used, consider adding whitespace around the operators with the lowest priority(ies)." Recommended:

x = x*2 - 1
hypot2 = x*x + y*y
c = (a+b) * (a-b)

http://www.python.org/dev/peps/pep-0008/#other-recommendations

@Eric89GXL
Eric89GXL added a note Nov 21, 2013

Ahh, your PEP8-fu is better than mine :) Thanks.

@Tillsten
Tillsten added a note Nov 22, 2013

This is a recent and very welcome change 👍

@endolith
endolith added a note Nov 22, 2013

Yeah I just realized this was changed a year ago, but I agree with the change. http://bugs.python.org/issue16239#msg173990

@Eric89GXL
Eric89GXL added a note Nov 22, 2013

Yeah, it looks quite a bit better. Now I get to try to convince the people I work with on other packages...

@rgommers
SciPy member
rgommers added a note Dec 16, 2013

This modifies an input parameter. k = real(k * prod(2*fs - z) / prod(2*fs - p)) would be better.

@rgommers
SciPy member
rgommers added a note Dec 16, 2013

Also, do fs = 2 * fs first, saves repeating 2*fs six times.

@endolith
endolith added a note Dec 17, 2013

Good point, I will fix these when I have time.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@Eric89GXL Eric89GXL and 1 other commented on an outdated diff Nov 21, 2013
scipy/signal/filter_design.py
@@ -1624,6 +1762,8 @@ def cheb1ap(N, rp):
defined as the point at which the gain first drops below -`rp`.
"""
+ if N == 0:
@Eric89GXL
Eric89GXL added a note Nov 21, 2013

Should we add a check for < 0 in case someone does something silly?

@endolith
endolith added a note Nov 22, 2013

Probably. Right now it doesn't even raise an error

In [491]: butter(0, 1)
Out[491]: (array([ 1.]), array([ 1.]))

In [492]: butter(-1, 1)
Out[492]: (array([ 1.]), array([ 1.]))

In [493]: butter(-10, 1)
Out[493]: (array([ 1.]), array([ 1.]))

Also probably should raise an error if analog==False and not 0 <= Wn <= 1? Right now it looks like it aliases, but I don't know why anyone would intentionally want to do this.

@endolith
endolith added a note Nov 22, 2013

It also accepts non-integer orders and spits out invalid things without complaining, etc. How much checking is too much?

@Eric89GXL
Eric89GXL added a note Nov 22, 2013

Could you do something like:

if not np.uint32(N) == N:
    raise ValueError(...)

I wonder if there is some trick like this that would take care of all the checking in one line.

@endolith
endolith added a note Nov 22, 2013

or this?

if abs(int(N)) != N:
    raise ValueError("Filter order must be a positive integer")

or should it round to the nearest integer instead?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@Eric89GXL

This looks nice (aside from my nitpicks), I'm excited about having ZPK implemented natively!

@endolith

still needs second-order sections to make high-order digital filters usable, though.

still needs sos

@Eric89GXL

Makes sense. Are you planning on including that with this PR, in a subsequent PR, or going to leave that to others? If you don't want to do it, I'll start reading up now on what will be required to implement SOS so I can be ready to give it a shot once you're done.

@endolith

Not in this PR.

I can work on it or someone else can. I started to work on it a while ago but didn't get very far, and I'm not even sure that's the correct or best way to do it. numpy.poly uses exact float equality to combine pairs instead of approximate equality, but it probably needs to be approximate for general usage. Maybe better to use some kind of least squared error fit thing to match up the poles, since very nearby poles might get incorrectly matched up if just based on sorting by real part?

@endolith endolith ENH: Add check for negative or fractional filter orders (besselap alr…
…eady does it)

Indent imports

Use string for error or else besselap says "not supported for order 3" when it means "not supported for order 3.5"
d9e1278
@coveralls

Coverage Status

Coverage remained the same when pulling d9e1278 on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@coveralls

Coverage Status

Coverage remained the same when pulling 8fbe568 on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@Eric89GXL

@endolith good question about pole matching. I'll see if I can find anything in Oppenheim & Schafer (Discrete-Time Signal Processing) to see if they have recommendations about SOS implementation.

@endolith

Yeah, the method I used in cplxpair of sorting and then combining based on approximate equality will probably work fine 99.99% of the time, maybe I'm worrying about nothing.

I don't know how to do the sos filtering with initial conditions, though. do they have to be passed through from one stage to the next or what? I don't even know how to use it with a single-section filter.

probably should discuss this on #2444 though

@coveralls

Coverage Status

Coverage remained the same when pulling 56e5ba6 on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@rgommers
SciPy member

Looks like a good solution for now to me. Is this done except for tests?

@ewmoore ewmoore and 1 other commented on an outdated diff Nov 27, 2013
scipy/signal/filter_design.py
+ # Duplicate poles and zeros and shift from baseband to +wo and -wo
+ p = concatenate((pbw + sqrt(pbw**2 - wo**2), pbw - sqrt(pbw**2 - wo**2)))
+ z = concatenate((zbw + sqrt(zbw**2 - wo**2), zbw - sqrt(zbw**2 - wo**2)))
+
+ # Move half of any zeros that were at infinity to the origin for BPF
+ z = append(z, zeros(degree))
+
+ # Cancel out gain change from frequency scaling
+ k *= bw**degree
+
+ return z, p, k
+
+
+def _zpklp2bs(z, p, k, wo=1.0, bw=1.0):
+ """
+ Transform a lowpass filter prototype to a highpass filter.
@ewmoore
SciPy member
ewmoore added a note Nov 27, 2013

highpass -> band-stop

@endolith
endolith added a note Nov 27, 2013

oops

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@endolith

Yes, I should still add some tests, but otherwise done.

@coveralls

Coverage Status

Coverage remained the same when pulling 4e98af0 on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@endolith endolith and 1 other commented on an outdated diff Dec 2, 2013
scipy/signal/filter_design.py
@@ -1624,6 +1894,10 @@ def cheb1ap(N, rp):
defined as the point at which the gain first drops below -`rp`.
"""
+ if abs(int(N)) != N:
+ raise ValueError("Filter order must be a nonnegative integer")
+ elif N == 0:
+ return numpy.array([]), numpy.array([]), 1
@endolith
endolith added a note Dec 2, 2013

now that I am writing tests for these degenerate cases, I'm not sure about this. for a 0-order cheby1 or ellip, should k=1, or k=rp?

@endolith
endolith added a note Dec 2, 2013

Testing with cheby1 and ellip at higher orders, the DC gain is 1 if N is odd, and -rp dB if N is even, so it would seem that since 0 is even, 0-order gain should also be -rp, but I'm not sure. This is what ellip already did before my changes, while cheb1y had a divide by zero error.

@rgommers
SciPy member
rgommers added a note Dec 2, 2013

Looks like it shouldn't be 1. It would help to give definitions of the filter; I guess from a quick look that rp equals 1/sqrt(1+eps**2) in the definition of https://en.wikipedia.org/wiki/Chebyshev_filter?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@endolith
endolith commented Dec 3, 2013

Ok, I added a lot of tests, including a bunch from web examples I could find because I don't have matlab and that way there's a reference to back them up anyway. Let me know if I did anything wrong.

Examples of matlab behavior were taken from:

@coveralls

Coverage Status

Coverage remained the same when pulling 8288cee on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@coveralls

Coverage Status

Coverage remained the same when pulling 124b496 on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@argriffing

fwiw there is also an exp10 function

Where? from numpy import exp10ImportError: cannot import name exp10

it's in scipy.special

@argriffing

fwiw I think assert_allclose is preferred except when comparing to zero

hmm. well that's not imported in this file, and it uses a relative tolerance instead of decimal points. for most of these tests decimal points make more sense I think, because of the precision of the original values we're comparing to?

OK this was just a reflex suggestion, because it's now the preferred general purpose comparison function. But if you already know about it and there is some reason to not use it then that's OK.

@rgommers rgommers and 1 other commented on an outdated diff Dec 16, 2013
scipy/signal/tests/test_filter_design.py
+ # https://github.com/scipy/scipy/pull/441
+ for N in range(7):
+ z, p, k = cheb2ap(N, 20)
+ assert_(isinstance(z, np.ndarray))
+ assert_(isinstance(p, np.ndarray))
+
+
+class TestEllipap(TestCase):
+
+ def test_output_type(self):
+ # Should consistently output arrays, not lists
+ # https://github.com/scipy/scipy/pull/441
+ for N in range(7):
+ z, p, k = ellipap(N, 1, 20)
+ assert_(isinstance(z, np.ndarray))
+ assert_(isinstance(p, np.ndarray))
@rgommers
SciPy member
rgommers added a note Dec 16, 2013

All the above tests can be combined into one with an extra for-loop:

for func in [cheb2ap, ellipap, ...]:
    ...
@endolith
endolith added a note Dec 27, 2013

What should the test class or test function be called in that case? Each of these classes should have other tests added to them, too.

@endolith
endolith added a note Dec 27, 2013

Actually, no, the function values are implicitly tested by cheby1, butter, etc. so the only thing that needs to be tested on the prototype functions themselves is their output type, not their values.

However, these functions have different parameters, so for func in list won't work. How about this?

    for func in (buttap, 
                 besselap, 
                 lambda N: cheb1ap(N, 1), 
                 lambda N: cheb2ap(N, 20),
                 lambda N: ellipap(N, 1, 20)):
        for N in range(7):
            z, p, k = func(N)
            assert_(isinstance(z, np.ndarray))
            assert_(isinstance(p, np.ndarray))
@rgommers
SciPy member
rgommers added a note Dec 30, 2013

That looks good to me.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers and 2 others commented on an outdated diff Dec 16, 2013
scipy/signal/tests/test_filter_design.py
+ assert_equal(cheb1ord(1, 1.2, 1, 80, analog=True)[0], 17)
+
+ # From http://www.mathworks.com/help/signal/ref/cheb1ord.html
+ n, Wp = cheb1ord(40/500, 150/500, 3, 60)
+ assert_equal(n, 4)
+ assert_array_almost_equal(Wp, 0.0800, decimal=4)
+
+ n, Wp = cheb1ord([60/500, 200/500], [50/500, 250/500], 3, 40)
+ assert_equal(n, 7)
+ assert_array_almost_equal(Wp, [0.1200, 0.4000], decimal=4)
+
+
+class TestCheb2ord(TestCase):
+
+ def test_basic(self):
+ # From http://www.mathworks.com/help/signal/ref/cheb2ord.html
@rgommers
SciPy member
rgommers added a note Dec 16, 2013

I don't think literally copying an example from the Matlab docs is a good idea - copyright and all. Can you change it?

@endolith
endolith added a note Dec 17, 2013

Well the intent is to show that we produce the same output for the same inputs, since it's meant to function identically. I would think that this falls under "facts and data are not subject to copyright", as I didn't copy any of the creative content like the text. But alternatively, someone with access to Matlab could generate some arbitrary examples to use instead, and then we'd have more precision for the values, too?

@rgommers
SciPy member
rgommers added a note Dec 17, 2013

Yeah, it's technically not copywriteable I think (but IANAL). On the other hand, it's something you really don't want to find out first-hand as a bunch of volunteers dealing with a large company:)

I may be able to find someone with Matlab access. Or just post it on the list, if it's not more than running an m-file it's simple for anyone to do.

@Eric89GXL
Eric89GXL added a note Dec 17, 2013

I have MATLAB access. If you ping me with your ideal Python test, I can run the equivalent in MATLAB and give you the results.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff Dec 16, 2013
scipy/signal/tests/test_filter_design.py
+
+class TestButter(TestCase):
+
+ def test_degenerate(self):
+ # 0-order filter is just a passthrough
+ b, a = butter(0, 1, analog=True)
+ assert_array_equal(b, [1])
+ assert_array_equal(a, [1])
+
+ # 1-order filter is same for all types
+ b, a = butter(1, 1, analog=True)
+ assert_array_almost_equal(b, [1])
+ assert_array_almost_equal(a, [1, 1])
+
+ def test_basic(self):
+ # Requires https://github.com/scipy/scipy/pull/3085 to pass
@rgommers
SciPy member
rgommers added a note Dec 16, 2013

This comment can be removed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff Dec 16, 2013
scipy/signal/tests/test_filter_design.py
+ wn = 0.01
+ z, p, k = butter(N, wn, 'low', analog=True, output='zpk')
+ assert_array_almost_equal([], z)
+ assert_(len(p) == N)
+ # All poles should be at distance wn from origin
+ assert_array_almost_equal(wn, abs(p))
+ assert_(all(np.real(p) <= 0)) # No poles in right half of S-plane
+ assert_array_almost_equal(wn**N, k)
+
+ for N in range(25):
+ wn = 0.01
+ z, p, k = butter(N, wn, 'high', analog=False, output='zpk')
+ assert_array_equal(np.ones(N), z) # All zeros exactly at DC
+ assert_(all(np.abs(p) <= 1)) # No poles outside unit circle
+
+ # From http://dsp.etfbl.net/filtri/aproksimacija
@rgommers
SciPy member
rgommers added a note Dec 16, 2013

This too, and similar comments. Links will break over time. Maybe comment on the PR "I used some examples from sites A, B and C to validate ...."

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
SciPy member

Do you have a reference for the zpk transforms used (preferably accessible online)? I could check _zpkbilinear but not the other ones. Would also be good to add a solid reference to one of the docstrings.

@rgommers
SciPy member

Test coverage looks good, although lp2lp and co are not run anymore. Would you mind adding some basic tests that exercise those functions?

@rgommers
SciPy member

I noticed by the way that Matlab has lp2lp and co in tf and ss format, but not zpk. Any idea why?

@rgommers
SciPy member

Would it make sense to take one of the plots in this PR (the bandpass one in the description for example) and add it to the iirfilter docstring?

@rgommers
SciPy member

Overall looks good, but not being an expert I'd like to have that reference to check. @Eric89GXL did you review the private functions?

@rgommers
SciPy member

Last request: could you add a few sentences to the 0.14.0 release notes? This is an important enough improvement to mention there I think.

@Eric89GXL

@rgommers I have not reviewed the private functions yet. This should be the sort of back-end work that good testing would cover, but I can take a look in the next few days to see if it makes sense.

@rgommers
SciPy member

Thanks, that would help. Tests indeed give a reasonable amount of confidence, but it's not a substitute for review.

@Eric89GXL

@endolith let me know if you have a couple MATLAB examples you want me to run and query the output of. I'll do that, let you implement the corresponding test, and then take a deeper look at the code. Sound good?

@endolith endolith MAINT: More descriptive variable names which don't modify input param…
…eters, check for improper transfer function in all transforms, calculate fs*2 only once
81b330c
@endolith

Ok, I fixed the issues with filter_design, and changed the variable names in the process, hopefully it's more clear.

As for references, there are bunch of bits and pieces around:

but moving the poles and zeros around isn't really spelled out explicitly anywhere that I could find, so I tried out stackedit.io by writing up the whole ugly derivation there. :)

I'll make a list of matlab tests next.

@endolith

in tf and ss format, but not zpk. Any idea why?

Maybe ss format so it can handle MIMO systems? Or maybe it's more numerically stable if input is already in ss format?

@coveralls

Coverage Status

Coverage remained the same when pulling 81b330c on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@coveralls

Coverage Status

Coverage remained the same when pulling ebe441c on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@rgommers
SciPy member

These changes look good to me.

@rgommers
SciPy member

@endolith I've rebased this branch and added a fix for the plot in the iirfilter docstring not showing up in the html docs: https://github.com/rgommers/scipy/tree/pr/3085. Could you take over that branch or cherry-pick the last commit?

@rgommers rgommers and 2 others commented on an outdated diff Dec 31, 2013
scipy/signal/filter_design.py
+
+ # Scale poles and zeros to desired bandwidth
+ z_lp = z * bw/2
+ p_lp = p * bw/2
+
+ # Square root needs to produce complex result, not NaN
+ z_lp = z_lp.astype(complex)
+ p_lp = p_lp.astype(complex)
+
+ # Duplicate poles and zeros and shift from baseband to +wo and -wo
+ z_bp = concatenate((z_lp + sqrt(z_lp**2 - wo**2),
+ z_lp - sqrt(z_lp**2 - wo**2)))
+ p_bp = concatenate((p_lp + sqrt(p_lp**2 - wo**2),
+ p_lp - sqrt(p_lp**2 - wo**2)))
+
+ # Move half of any zeros that were at infinity to the origin for BPF
@rgommers
SciPy member
rgommers added a note Dec 31, 2013

half of should be removed from this comment I think, you're adding degree zeros here. It's the same s^{n-m} as in _zpklp2hp.

@endolith
endolith added a note Dec 31, 2013

Hmm, I mean that it's doubling the number of zeros at infinity (degree*2), and then only moving half of them (degree) to the origin.

"Move degree zeros to origin, leaving degree zeros at infinity"?

@Eric89GXL
Eric89GXL added a note Dec 31, 2013

That is clearer to me.

@rgommers
SciPy member
rgommers added a note Dec 31, 2013

Sounds fine to me.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
SciPy member

@endolith your stackedit write-up is quite nice (and useful site too, had never seen it before). I checked a number (not all) of your derivations, looks good and matches the code in this PR. I've also tried some of the examples in other issues, all work now. So this is in good shape. When the update of the tests is done it can be merged.

@rgommers
SciPy member

Using the top part of the stackedit document as the start of a filter design section for http://docs.scipy.org/doc/scipy/reference/tutorial/signal.html might make sense. Right now the tutorial has exactly nothing on filter design.

@endolith

@Eric89GXL Ok here's a bunch of tests for matlab, trying to cover each feature with some arbitrary values:

Use format long e or something like num2str(p, 18) to dump the output arrays with high precision?

[z, p, k] = besself(0, 0.3) % degenerate (octave doesn't allow this)
[z, p, k] = besself(1, 0.3) % first order is same for all types
[z, p, k] = besself(24, 0.43) % high even order
[z, p, k] = besself(23, 0.56) % high odd order
% with transfer function conversion,  without digital conversion
[b, a] = besself(4, 300, 's')

[z, p, k] = butter(0, 0.3) % degenerate (octave doesn't allow this)
[z, p, k] = butter(1, 0.3) % first order is same for all types
[z, p, k] = butter(28, 0.43, 'high') % high even order
[z, p, k] = butter(27, 0.56, 'high') % high odd order
[z, p, k] = butter(8, [0.25 0.33]) % bandpass
[z, p, k] = butter(7, [0.45 0.56], 'stop') % bandstop
% with transfer function conversion,  without digital conversion
[b, a] = butter(4, [100 300], 'bandpass', 's')
[z, p, k] = butter(4, [90.5 110.5], 's') % Paarmann Example 9.11

[z, p, k] = cheby1(0, 1, 0.6) % degenerate (octave doesn't allow this)
[z, p, k] = cheby1(1, 0.1, 0.3) % first order is same for all types
[z, p, k] = cheby1(24, 0.7, 0.2, 'high') % high even order
[z, p, k] = cheby1(23, 0.8, 0.8, 'high') % high odd order
[z, p, k] = cheby1(8, 1, [0.3 0.4]) % bandpass
[z, p, k] = cheby1(7, 1, [0.5 0.6], 'stop') % bandstop
% with transfer function conversion,  without digital conversion
[b, a] = cheby1(5, 0.9, [210 310], 'stop', 's')
[z, p, k] = cheby1(10, 1, 1000, 'high', 's') % Paarmann Example 9.6

[z, p, k] = cheby2(0, 40, 0.3) % degenerate (octave doesn't allow this)
[z, p, k] = cheby2(1, 50, 0.3) % first order is same for all types
[z, p, k] = cheby2(26, 60, 0.3, 'high') % high even order
[z, p, k] = cheby2(25, 80, 0.5, 'high') % high odd order
[z, p, k] = cheby2(9, 40, [0.07 0.2]) % bandpass
[z, p, k] = cheby2(6, 55, [0.1 0.9], 'stop') % bandstop
% with transfer function conversion,  without digital conversion
[b, a] = cheby2(5, 20, [2010 2100], 'stop', 's')

[z, p, k] = ellip(0, 1, 45, 0.1) % degenerate (octave doesn't allow this)
[z, p, k] = ellip(1, 1, 55, 0.3) % first order is same for all types
[z, p, k] = ellip(24, 1, 80, 0.3, 'high') % high even order
[z, p, k] = ellip(23, 1, 70, 0.5, 'high') % high odd order
[z, p, k] = ellip(7, 1, 40, [0.07 0.2]) % bandpass
[z, p, k] = ellip(8, 1, 65, [0.2 0.4], 'stop') % bandstop
% with transfer function conversion,  without digital conversion
[b, a] = ellip(5, 1, 40, [201 240], 'stop', 's')
[z, p, k] = ellip(5, 1, 75, [90.5 110.5], 's') % Paarmann Example 9.12

% minimum order calculations (digital):
[n, Wn] = buttord(0.2, 0.3, 3, 60) % lowpass
[n, Wn] = buttord(0.3, 0.2, 3, 70) % highpass
[n, Wn] = buttord([0.2, 0.5], [0.1, 0.6], 3, 80) % bandpass
[n, Wn] = buttord([0.1, 0.6], [0.2, 0.5], 3, 90) % bandstop

[n, Wn] = cheb1ord(0.2, 0.3, 3, 60) % lowpass
[n, Wn] = cheb1ord(0.3, 0.2, 3, 70) % highpass
[n, Wn] = cheb1ord([0.2, 0.5], [0.1, 0.6], 3, 80) % bandpass
[n, Wn] = cheb1ord([0.1, 0.6], [0.2, 0.5], 3, 90) % bandstop

[n, Wn] = cheb2ord(0.2, 0.3, 3, 60) % lowpass
[n, Wn] = cheb2ord(0.3, 0.2, 3, 70) % highpass
[n, Wn] = cheb2ord([0.2, 0.5], [0.1, 0.6], 3, 80) % bandpass
[n, Wn] = cheb2ord([0.1, 0.6], [0.2, 0.5], 3, 90) % bandstop

[n, Wn] = ellipord(0.2, 0.3, 3, 60) % lowpass
[n, Wn] = ellipord(0.3, 0.2, 3, 70) % highpass
[n, Wn] = ellipord([0.2, 0.5], [0.1, 0.6], 3, 80) % bandpass
[n, Wn] = ellipord([0.1, 0.6], [0.2, 0.5], 3, 90) % bandstop

% analog
[n, Wn] = buttord(200, 600, 3, 60, 's') % lowpass
[n, Wn] = cheb1ord(700, 100, 3, 70, 's') % highpass
[n, Wn] = cheb2ord([20, 50], [10, 60], 3, 80, 's') % bandpass
[n, Wn] = ellipord([1000, 6000], [2000, 5000], 3, 90, 's') % bandstop

% octave does not match matlab's documentation:
[n, Ws] = cheb2ord(40/500, 150/500, 3, 60)
[n, Ws] = cheb2ord([60/500, 200/500], [50/500, 250/500], 3, 40)
@endolith

I've rebased this branch and added a fix for the plot in the iirfilter docstring not showing up in the html docs:

@rgommers

That's weird that it doesn't show up. The other docstrings like bessel do show up, right? I copied one of those and modified it for iirfilter, so I would expect it to work the same way?

@rgommers
SciPy member
rgommers commented Jan 1, 2014

The issue with that example was that it was missing the ... at the beginning of the continuation line. Which is invalid doctest syntax.

edit: in addition I changed the syntax to use the object-oriented MPL interface, which should generally be preferred over the matlab-style one.

@Eric89GXL

@endolith I had to comment out all of the "degenerate" cases in your example, since they cause MATLAB to throw errors complaining that the order had to be positive non-zero. The script I ran (with very slight modifications) and output are here:

https://gist.github.com/Eric89GXL/8239281

@endolith
endolith commented Jan 7, 2014

Thanks! I converted these to unit tests, and there are some discrepancies, so I'll try to figure out what's causing those. I think #3081 will fix some. Also, I didn't realize besself is always analog, so my 's' example is nonsense. I'll make some new examples for that.

@Eric89GXL

No problem, it took all of 2 minutes to tweak and run that script at my end. Let me know if there are other things it would be useful for me to run.

@endolith endolith DOC: Make plot in signal.iirfilter docstring appear,
add note about highpass symmetry, change comment wording
63564c3
@endolith

So the order functions don't produce the same output. For instance, cheb2ord([0.2, 0.5], [0.1, 0.6], 3, 80) produces 9, [ 0.1488, 0.5975], while matlab gives you back the original frequencies you entered: n = 9
Wn = [0.1, 0.6]
I don't know how they work, so I made them pass by setting decimal=1 (approximately the same) and it can be dealt with in another PR?

Also commented out ellipord([1000, 6000], [2000, 5000], 3, 90, True), which produces 8, [ 1666.6666, 5999.9999] while matlab produces 9, [1000, 6000]. I don't know why, but these aren't really related to this PR. Is that ok? Or should I use knownfailureif?

@endolith

@Eric89GXL can you redo these:

[z, p, k] = ellip(5, 1, 40, [201 240], 'stop', 's')
[z, p, k] = cheby2(5, 20, [2010 2100], 'stop', 's')
[z, p, k] = cheby1(23, 0.8, 0.3, 'high') % high odd order
[z, p, k] = besself(1, 1)
[z, p, k] = besself(24, 100)
[z, p, k] = besself(23, 1000)
[b, a] = besself(4, 300)
@rgommers
SciPy member

@endolith I think adding the test and decorating them with knownfailureif makes sense to me. Plus opening a separate issue then that explains that the difference with Matlab needs investigating.

@rgommers
SciPy member

You say gh-3081 will fix some issues, but you closed that PR again?

@endolith

Yeah sorry, I merged those changes into this PR, but haven't pushed them yet.

@rgommers
SciPy member

OK thanks.

@endolith

Even with @dec.knownfailureif(True, running the test file says FAILED (errors=1). I thought it would ignore it? Should I use @dec.skipif( instead? Then it says OK (SKIP=1)

@rgommers
SciPy member

Are you running the tests by python test_filterdesign.py? That bypasses part of the numpy test machinery, and KnownFailure is defined in numpy.testing. So nose reports knownfailures as errors. If you use signal.test() it will work as expected.

@rgommers
SciPy member

I'll send a PR to fix that issue in numpy in a minute.

endolith added some commits Jan 13, 2014
@endolith endolith BUG: buttap() should produce a purely real pole
Previously buttap for odd orders produced a pole like -1+6.1e-17j rather than -1+0j exactly.  This made poly(p) complex, though all the imaginary parts were negligible.  This then made the ba representation complex, though it is a purely real filter.

Now the calculation is symmetrical, so the middle term works out to -exp(0j), which is exactly -1+0j.
d103964
@endolith endolith BUG: Change cheb1ap() to produce a purely real pole
Change theta of calculation to make it symmetrical about 0, so that the lone pole for odd-ordered filters is an exact real number, and the pole pairs are exact conjugates, so that ba representation (from numpy.poly) is purely real.

Also some simplifications:
arcsinh(1 / eps) is equivalent to log((1 + sqrt(1 + eps * eps)) / eps)
-sinh(mu + 1j*theta) is equivalent to -sinh(mu) * cos(theta) + 1j * cosh(mu) * -sin(theta)
which is the rotated version of the original
-sinh(mu) * sin(theta) + 1j * cosh(mu) * cos(theta)
265fb61
@endolith endolith BUG: Change cheb2ap() to produce a purely real pole
Change theta of calculation to make it symmetrical about 0, so that the lone pole for odd-ordered filters is an exact real number, and the pole pairs are exact conjugates, so that ba representation (from numpy.poly) is purely real.
29aea64
@endolith endolith TST: Add test for purely real output of IIR filters 061f0d5
@endolith endolith TST: Add tests for filter design functions 5a1a36e
@endolith endolith BUG: int arguments to filter design functions result in np.int wrapar…
…ound, convert to float first like the original lp2lp
ab9b4fc
@endolith endolith TST: Test for integer frequencies causing wraparound 175a094
@endolith endolith BUG: zpkt2f sometimes produces complex output with negligible imagina…
…ry parts

Fixed in numpy/numpy#4073, but we can't depend on a specific numpy version, so copying it here, too.
e26799c
@endolith endolith TST: Add test for real ba output type of IIR filters e30c56a
@endolith endolith ENH: Raise an error if digital filter frequencies are out of range (O…
…ctave does this, too)
4dc5034
@endolith

From #3081 (comment)

I think this needs a fix that doesn't rely on a change in numpy, given that we have to support numpy >= 1.5.1.

Ok, I copied the code from my change to numpy.poly over to zpk2tf, so it will work regardless of numpy version. In the future when numpy.poly does the right thing, it should probably be removed?

@endolith

Ok added all the tests for filter design functions, found and fixed some bugs, and incorporated #3081. Some of the 'ba' tests were actually failing because Matlab produces inaccurate results. :) For instance, ellip(5, 1, 40, [201, 240], 'stop', True)[0] should be:

array([  1.0000e+00,   0.0000e+00,   2.4266e+05,   0.0000e+00,
         2.3482e+10,   0.0000e+00,   1.1328e+15,   0.0000e+00,
         2.7240e+19,   0.0000e+00,   2.6124e+23])

but matlab produces

array([  1.0000e+00,   1.7435e-13,   2.4266e+05,   3.4594e-08,
         2.3482e+10,   2.5592e-03,   1.1328e+15,   8.3632e+01,
         2.7240e+19,   1.0187e+06,   2.6124e+23])
@coveralls

Coverage Status

Coverage remained the same when pulling 4dc5034 on endolith:zpkiirfilter into 572aaf0 on scipy:master.

endolith added some commits Jan 21, 2014
@endolith endolith ENH: Raise an error if multiple frequencies are fed to lowpass or hig…
…hpass filter design
aa1a573
@endolith endolith TST: Test that invalid Wn values raise an error 1518f17
@endolith endolith TST: Check frequency response for filter design tools
The *ord functions design a filter to match a specification, so we should check whether the resulting frequency response actually meets the specification.

Also change the tests for functions that don't match matlab, considering their current behavior to be correct.  If scipy#3235 changes their behavior, the tests will be changed to match.
4afdc53
@coveralls

Coverage Status

Coverage remained the same when pulling 4afdc53 on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@endolith endolith TST: Add some basic tests for lp2lp, etc.
Using some examples found online
ae94ba0
@coveralls

Coverage Status

Coverage remained the same when pulling ae94ba0 on endolith:zpkiirfilter into 572aaf0 on scipy:master.

@rgommers
SciPy member

Since this is ready and it's about time to merge it, in it goes. Last few added commits look good. Thanks a lot @endolith.

@rgommers rgommers closed this Jan 31, 2014
@rgommers rgommers reopened this Jan 31, 2014
@rgommers rgommers merged commit 233ad82 into scipy:master Jan 31, 2014

1 check passed

Details default The Travis CI build passed
@rgommers
SciPy member

gh-3235 should be rebased and the tests updated there. I'll comment on that PR.

@rgommers
SciPy member

The 0.14.0 release notes update still has to be done, I can take that along when preparing the first beta.

@endolith endolith deleted the endolith:zpkiirfilter branch Feb 2, 2014
@ghego
ghego commented May 26, 2014

I'm stilly having problems with higher order digital filters. Am I doing something wrong or the issue is not fixed yet?

http://stackoverflow.com/questions/21862777/bandpass-butterworth-filter-frequencies-in-scipy

@pv
SciPy member
pv commented May 26, 2014

@ghego: AFAIK, the problem is that the b, a representation itself is numerically unstable, and your code uses it. The iirfilter function was fixed to use zpk internally, and it seems butter and some other functions can take kwargs output="zpk"

@ghego
ghego commented May 27, 2014

Yes, but how do you actually use the zpk representation to filter a signal?

@endolith

You convert it into second-order sections and then use those in b,a format. I posted a link on your stackoverflow question to code cloned from Octave that does this.

@ghego
ghego commented May 27, 2014

oh, ok, thanks. So there's no way to use zpk form directly?

@endolith

No, but I think that would be a good thing to add: #3259

Matlab is the same, as far as I know, like Matlab's freqz takes ba or sos input, but not zpk or ss.

@ghego
ghego commented May 27, 2014
@endolith endolith commented on the diff Jun 14, 2016
scipy/signal/filter_design.py
@@ -311,6 +312,29 @@ def zpk2tf(z, p, k):
else:
b = k * poly(z)
a = atleast_1d(poly(p))
+
+ # Use real output if possible. Copied from numpy.poly, since
@endolith
endolith added a note Jun 14, 2016

This should be in _numpy_compat.py then?

@e-q
e-q added a note Jun 14, 2016

Yes, that would be good.

@endolith
endolith added a note Jun 15, 2016 edited

actually numpy/numpy#4073 never got merged :D

so I can't say if NumpyVersion > than the version it got merged in

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.