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

[Josh Lawrence]

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?

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