[SciPy-User] how to treat an invalid value, in signal/filter_design.py

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view

[SciPy-User] how to treat an invalid value, in signal/filter_design.py

R Schumacher
In an attempt to computationally invert the effect of an analog RC filter on a data set and reconstruct the signal prior to the analog front end, a co-worker suggested: "Mathematically, you just reverse the a and b parameters. Then the zeros become the poles, but if the new poles are not inside the unit circle, the filter is not stable."

So then to "stabilize" the poles' issue seen, I test for the DIV/0 error and set it to 2./N+0.j in scipy/signal/filter_design.py ~ line 244
    d = polyval(a[::-1], zm1)
    if d[0]==0.0+0.j:
        d[0] = 2./N+0.j
    h = polyval(b[::-1], zm1) / d

- Question is, is this a mathematically valid treatment?
- Is there a better way to invert a Butterworth filter, or work with the DIV/0 that occurs without modifying the signal library?

I noted d[0] > 2./N+0.j makes the zero bin result spike low; 2/N gives a reasonable "extension" of the response curve.
The process in general causes a near-zero offset however, which I remove with a high pass now; In an full FFT of a ~megasample one can see that the first 5 bins have run away.

An example attached...

Ray Schumacher

Ray Schumacher
PO Box 182, Pine Valley, CA 91962

SciPy-User mailing list
[hidden email]

NBP_test.py (4K) Download Attachment