[SciPy-Dev] scipy.signal.normalize problems?
Josh Lawrence
josh.k.lawrence at gmail.com
Wed May 30 13:45:33 EDT 2012
Okay,
I might take a few days to get around to writing a test for it. But I'll update the pull request when I do.
--Josh
On May 30, 2012, at 1:38 PM, Ralf Gommers wrote:
>
>
> On Wed, May 30, 2012 at 7:20 PM, Josh Lawrence <josh.k.lawrence at gmail.com> wrote:
> Okay,
>
> I made a pull request:
>
> https://github.com/scipy/scipy/pull/233
>
> Do I need to file a ticket on Trac or is the PR good enough?
>
> The PR should be good enough. Two suggestions though for the PR:
> - in a comment, link to this discussion (gmane archive)
> - please add a unit test. You already have an example that didn't work before, and now gives correct results (or at least matches Matlab). You can convert that into a test, with the output of normalize() checked against the hardcoded correct numerical values.
>
> Thanks,
> Ralf
>
>
> --Josh
>
> On May 29, 2012, at 4:08 PM, Ralf Gommers wrote:
>
>>
>>
>> On Mon, May 21, 2012 at 9:11 PM, Josh Lawrence <josh.k.lawrence at gmail.com> wrote:
>> If I change the allclose lines to have
>>
>> allclose(0, outb[:,0], atol=1e-14)
>>
>> it works. I think that is what the original goal was, anyways. From the documentation of allclose, what I have written above result in ensuring
>>
>> np.abs(0 - outb[:,0]) > atol + rtol * np.abs(outb[:,0])
>>
>> with rtol defaulting to 1e-5. I'm still not sure about how things have been written for the 'b' argument of normalize being rank-2, so I can't guarantee that my fix makes things work for that. Should I file a bug report/submit a patch/send a git pull request, etc?
>>
>> Sounds like you understood the issue and have a fix that does what you expect (matches Matlab), so a pull request would be nice.
>>
>> For the rank-2 thing, you can just note it on the PR if you can't figure it out.
>>
>> Ralf
>>
>>
>> Cheers,
>>
>> Josh Lawrence
>>
>> On May 21, 2012, at 2:45 PM, Josh Lawrence wrote:
>>
>> > Hey all,
>> >
>> > I've been having some problems designing a Chebyshev filter and I think I have narrowed down the hang-up to scipy.signal.normalize. I think what's going on in my case is that line 286 of filter_design.py (the first allclose call in the normalize function) is producing a false positive. Here's the function definition:
>> >
>> > def normalize(b, a):
>> > """Normalize polynomial representation of a transfer function.
>> >
>> > If values of b are too close to 0, they are removed. In that case, a
>> > BadCoefficients warning is emitted.
>> > """
>> > b, a = map(atleast_1d, (b, a))
>> > if len(a.shape) != 1:
>> > raise ValueError("Denominator polynomial must be rank-1 array.")
>> > if len(b.shape) > 2:
>> > raise ValueError("Numerator polynomial must be rank-1 or"
>> > " rank-2 array.")
>> > if len(b.shape) == 1:
>> > b = asarray([b], b.dtype.char)
>> > while a[0] == 0.0 and len(a) > 1:
>> > a = a[1:]
>> > outb = b * (1.0) / a[0]
>> > outa = a * (1.0) / a[0]
>> > if allclose(outb[:, 0], 0, rtol=1e-14): <------------------ Line 286
>> > warnings.warn("Badly conditioned filter coefficients (numerator): the "
>> > "results may be meaningless", BadCoefficients)
>> > while allclose(outb[:, 0], 0, rtol=1e-14) and (outb.shape[-1] > 1):
>> > outb = outb[:, 1:]
>> > if outb.shape[0] == 1:
>> > outb = outb[0]
>> > return outb, outa
>> >
>> > I marked line 286. If I reproduce all the steps carried out by scipy.signal.iirdesign, I end up with a (b, a) pair which results of scipy.signal.lp2lp and looks like this:
>> >
>> > In [106]: b_lp2
>> > Out[106]: array([ 1.55431359e-06+0.j])
>> >
>> > In [107]: a_lp2
>> > Out[107]:
>> > array([ 1.00000000e+00 +0.00000000e+00j,
>> > 3.46306104e-01 -2.01282794e-16j,
>> > 2.42572185e-01 -6.08207573e-17j,
>> > 5.92946943e-02 +0.00000000e+00j,
>> > 1.82069156e-02 +5.55318531e-18j,
>> > 2.89328123e-03 +0.00000000e+00j,
>> > 4.36566281e-04 -2.95766719e-19j,
>> > 3.50842810e-05 -3.19180568e-20j, 1.64641246e-06 -1.00966301e-21j])
>> >
>> > scipy.signal.iirdesign takes b_lp2, a_lp2 (my local variable names to keep track of what's going on) and runs them through scipy.signal.bilinear (in filter_design.py bilinear is called on line 624 within iirfilter. iirdesign calls iirfilter which calls bilinear). Inside bilinear, normalize is called on line 445. I've made my own class with bilinear copied and pasted from filter_design.py to test things. In bilinear, the input to normalize is given by
>> >
>> > b = [ 1.55431359e-06 1.24345087e-05 4.35207804e-05 8.70415608e-05
>> > 1.08801951e-04 8.70415608e-05 4.35207804e-05 1.24345087e-05
>> > 1.55431359e-06]
>> > a = [ 72269.02590913 -562426.61430468 1918276.19173089 -3745112.83646825
>> > 4577612.13937628 -3586970.61385926 1759651.18184723 -494097.93515708
>> > 60799.46134722]
>> >
>> > In normalize, right before the allclose() call, outb is defined by
>> >
>> > outb = [[ 2.15073272e-11 1.72058618e-10 6.02205162e-10 1.20441032e-09
>> > 1.50551290e-09 1.20441032e-09 6.02205162e-10 1.72058618e-10
>> > 2.15073272e-11]]
>> >
>> > From what I read of the function normalize, it should only evaluate true if all of the coefficients in outb are smaller than 1e-14. However, that's not what is going on. I have access to MATLAB and if I go through the same design criteria to design a chebyshev filter, I get the following:
>> >
>> > b =
>> >
>> > 1.0e-08 *
>> >
>> > Columns 1 through 5
>> >
>> > 0.002150733144728 0.017205865157826 0.060220528052392 0.120441056104784 0.150551320130980
>> >
>> > Columns 6 through 9
>> >
>> > 0.120441056104784 0.060220528052392 0.017205865157826 0.002150733144728
>> >
>> > which matches up rather well for several significant figures.
>> >
>> > I apologize if this is not clearly explained, but I'm not sure what to do. I tried messing around with the arguments to allclose (switching it to be allclose(0, outb[:,0], ...), or changing the keyword from rtol to atol). I am also not sure why normalize is setup to run on rank-2 arrays. I looked through filter_design.py and none of the functions contained in it send a rank-2 array to normalize from what I can tell. Any thoughts?
>> >
>> > Cheers,
>> >
>> > Josh Lawrence
>>
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20120530/109e5861/attachment.html>
More information about the SciPy-Dev
mailing list