[SciPy-Dev] scipy.signal.normalize problems?
Ralf Gommers
ralf.gommers at googlemail.com
Wed May 30 13:38:45 EDT 2012
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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20120530/f885d99c/attachment.html>
More information about the SciPy-Dev
mailing list