[Scipy-svn] r4662 - in trunk/scipy/signal: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Aug 21 15:06:55 EDT 2008
Author: stefan
Date: 2008-08-21 14:05:50 -0500 (Thu, 21 Aug 2008)
New Revision: 4662
Modified:
trunk/scipy/signal/signaltools.py
trunk/scipy/signal/tests/test_signaltools.py
Log:
Fix calculation of Chebyshev polynomial in chebwin. Closes #581.
Modified: trunk/scipy/signal/signaltools.py
===================================================================
--- trunk/scipy/signal/signaltools.py 2008-08-21 00:35:44 UTC (rev 4661)
+++ trunk/scipy/signal/signaltools.py 2008-08-21 19:05:50 UTC (rev 4662)
@@ -854,30 +854,31 @@
M = M+1
# compute the parameter beta
- beta = cosh(1.0/(M-1.0)*arccosh(10**(at/20.)))
+ order = M - 1.0
+ beta = cosh(1.0/order * arccosh(10**(abs(at)/20.)))
k = r_[0:M]*1.0
x = beta*cos(pi*k/M)
#find the window's DFT coefficients
- p = zeros(x.shape) * 1.0
- for i in range(len(x)):
- if x[i] < 1:
- p[i] = cos((M - 1) * arccos(x[i]))
- else:
- p[i] = cosh((M - 1) * arccosh(x[i]))
+ # Use analytic definition of Chebyshev polynomial instead of expansion
+ # from scipy.special. Using the expansion in scipy.special leads to errors.
+ p = zeros(x.shape)
+ p[x > 1] = cosh(order * arccosh(x[x > 1]))
+ p[x < -1] = (1 - 2*(order%2)) * cosh(order * arccosh(-x[x < -1]))
+ p[numpy.abs(x) <=1 ] = cos(order * arccos(x[numpy.abs(x) <= 1]))
# Appropriate IDFT and filling up
# depending on even/odd M
if M % 2:
- w = real(fft(p));
- n = (M + 1) / 2;
- w = w[:n] / w[0];
+ w = real(fft(p))
+ n = (M + 1) / 2
+ w = w[:n] / w[0]
w = concatenate((w[n - 1:0:-1], w))
else:
p = p * exp(1.j*pi / M * r_[0:M])
- w = real(fft(p));
- n = M / 2 + 1;
- w = w / w[1];
- w = concatenate((w[n - 1:0:-1], w[1:n]));
+ w = real(fft(p))
+ n = M / 2 + 1
+ w = w / w[1]
+ w = concatenate((w[n - 1:0:-1], w[1:n]))
if not sym and not odd:
w = w[:-1]
return w
Modified: trunk/scipy/signal/tests/test_signaltools.py
===================================================================
--- trunk/scipy/signal/tests/test_signaltools.py 2008-08-21 00:35:44 UTC (rev 4661)
+++ trunk/scipy/signal/tests/test_signaltools.py 2008-08-21 19:05:50 UTC (rev 4662)
@@ -55,5 +55,48 @@
assert_array_equal(signal.order_filter([1,2,3],[1,0,1],1),
[2,3,2])
+class TestChebWin:
+ def test_cheb_odd(self):
+ cheb_odd_true = array([0.200938, 0.107729, 0.134941, 0.165348,
+ 0.198891, 0.235450, 0.274846, 0.316836,
+ 0.361119, 0.407338, 0.455079, 0.503883,
+ 0.553248, 0.602637, 0.651489, 0.699227,
+ 0.745266, 0.789028, 0.829947, 0.867485,
+ 0.901138, 0.930448, 0.955010, 0.974482,
+ 0.988591, 0.997138, 1.000000, 0.997138,
+ 0.988591, 0.974482, 0.955010, 0.930448,
+ 0.901138, 0.867485, 0.829947, 0.789028,
+ 0.745266, 0.699227, 0.651489, 0.602637,
+ 0.553248, 0.503883, 0.455079, 0.407338,
+ 0.361119, 0.316836, 0.274846, 0.235450,
+ 0.198891, 0.165348, 0.134941, 0.107729,
+ 0.200938])
+
+ cheb_odd = signal.chebwin(53, at=-40)
+ assert_array_almost_equal(cheb_odd, cheb_odd_true, decimal=4)
+
+ def test_cheb_even(self):
+ cheb_even_true = array([0.203894, 0.107279, 0.133904,
+ 0.163608, 0.196338, 0.231986,
+ 0.270385, 0.311313, 0.354493,
+ 0.399594, 0.446233, 0.493983,
+ 0.542378, 0.590916, 0.639071,
+ 0.686302, 0.732055, 0.775783,
+ 0.816944, 0.855021, 0.889525,
+ 0.920006, 0.946060, 0.967339,
+ 0.983557, 0.994494, 1.000000,
+ 1.000000, 0.994494, 0.983557,
+ 0.967339, 0.946060, 0.920006,
+ 0.889525, 0.855021, 0.816944,
+ 0.775783, 0.732055, 0.686302,
+ 0.639071, 0.590916, 0.542378,
+ 0.493983, 0.446233, 0.399594,
+ 0.354493, 0.311313, 0.270385,
+ 0.231986, 0.196338, 0.163608,
+ 0.133904, 0.107279, 0.203894])
+
+ cheb_even = signal.chebwin(54, at=-40)
+ assert_array_almost_equal(cheb_even, cheb_even_true, decimal=4)
+
if __name__ == "__main__":
run_module_suite()
More information about the Scipy-svn
mailing list