[SciPy-Dev] scipy.signal.normalize problems?

Tue May 29 16:08:27 EDT 2012

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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20120529/31767607/attachment.html>