scipy.signal.normalize problems?
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
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
On Mon, May 21, 2012 at 9:11 PM, Josh Lawrence <josh.k.lawrence@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@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
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? --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@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@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Wed, May 30, 2012 at 7:20 PM, Josh Lawrence <josh.k.lawrence@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@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@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
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@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@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@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
I updated the PR to link to this discussion and I included a unit-test. I applied the changes to the current scipy-master and running scipy.signal.test() (where the new test is located) produces no errors or failures. Will this be backported to 0.11? --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@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@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@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Thu, Jun 28, 2012 at 7:01 PM, Josh Lawrence <josh.k.lawrence@gmail.com>wrote:
I updated the PR to link to this discussion and I included a unit-test. I applied the changes to the current scipy-master and running scipy.signal.test() (where the new test is located) produces no errors or failures.
Will this be backported to 0.11?
Could one of the scipy.signal maintainers please have a look at this fix? https://github.com/scipy/scipy/pull/233 Backporting to 0.11.x is no problem. Ralf
Sorry, I had my git repository all messed up. Here is the new pull request. https://github.com/scipy/scipy/pull/259 --Josh On Jun 28, 2012, at 1:06 PM, Ralf Gommers wrote:
On Thu, Jun 28, 2012 at 7:01 PM, Josh Lawrence <josh.k.lawrence@gmail.com> wrote: I updated the PR to link to this discussion and I included a unit-test. I applied the changes to the current scipy-master and running scipy.signal.test() (where the new test is located) produces no errors or failures.
Will this be backported to 0.11?
Could one of the scipy.signal maintainers please have a look at this fix? https://github.com/scipy/scipy/pull/233
Backporting to 0.11.x is no problem.
Ralf _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
participants (2)
-
Josh Lawrence -
Ralf Gommers