BUG: Use zpk in IIR filter design to improve accuracy #3085
Conversation
(since 'bp' is there and functions like 'lp2hp' imply it)
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.
Previously it said "IndexError: 0-d arrays can't be indexed" Combine bandpass and bandstop calculations for DRY Also add quotes to error messages
|
Ahh @endolith I didn't see this PR, ignore my questions in a (much) older issue. |
| 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, |
There was a problem hiding this comment.
With parens you might as well indent to the open paren, PEP8.
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
also I thought it was discouraged to use r_ except as a shortcut while working on the command line
|
This looks nice (aside from my nitpicks), I'm excited about having ZPK implemented natively! |
|
still needs second-order sections to make high-order digital filters usable, though. |
|
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. |
|
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? |
…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"
|
@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. |
|
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 |
|
Looks like a good solution for now to me. Is this done except for tests? |
|
|
||
| def _zpklp2bs(z, p, k, wo=1.0, bw=1.0): | ||
| """ | ||
| Transform a lowpass filter prototype to a highpass filter. |
|
Yes, I should still add some tests, but otherwise done. |
| if abs(int(N)) != N: | ||
| raise ValueError("Filter order must be a nonnegative integer") | ||
| elif N == 0: | ||
| return numpy.array([]), numpy.array([]), 1 |
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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?
…hpass filter design
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.
Using some examples found online
|
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. |
BUG: Use zpk in IIR filter design to improve accuracy
|
gh-3235 should be rebased and the tests updated there. I'll comment on that PR. |
|
The 0.14.0 release notes update still has to be done, I can take that along when preparing the first beta. |
|
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 |
|
@ghego: AFAIK, the problem is that the |
|
Yes, but how do you actually use the zpk representation to filter a signal? |
|
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. |
|
oh, ok, thanks. So there's no way to use zpk form directly? |
|
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. |
|
Makes sense. Thanks! On Tue, May 27, 2014 at 9:12 AM, endolith notifications@github.com wrote:
Francesco Mosconi, PhD |
| @@ -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 | |||
There was a problem hiding this comment.
This should be in _numpy_compat.py then?
There was a problem hiding this comment.
actually numpy/numpy#4073 never got merged :D
so I can't say if NumpyVersion > than the version it got merged in
Removed code defending against old Numpy: [1] was included in Numpy 1.12.0, and minimum Numpy now is 1.23.5. See [2] for context. [1] numpy/numpy#4073 [2] scipy#3085 (comment)
Removed code defending against old Numpy: [1] was included in Numpy 1.12.0, and minimum Numpy now is 1.23.5. See [2] for context. [1] numpy/numpy#4073 [2] scipy#3085 (comment)
…test for complex k, real p, z (#22213) Removed code defending against old Numpy: [1] was included in Numpy 1.12.0, and minimum Numpy now is 1.23.5. See [2] for context. [1] numpy/numpy#4073 [2] #3085 (comment)
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:
After:
21-pole analog Butterworth lowpass at Wn=0.01
Before:
After:
cheby2 bandpass
Before:
After:
