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