[Scipy-svn] r4662 - in trunk/scipy/signal: . tests

[scipy-svn at scipy.org]

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()