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

scipy-svn@scip... scipy-svn@scip...
Sun Jan 17 22:56:47 CST 2010


Author: josef
Date: 2010-01-17 22:56:46 -0600 (Sun, 17 Jan 2010)
New Revision: 6205

Modified:
   trunk/scipy/signal/signaltools.py
   trunk/scipy/signal/tests/test_signaltools.py
Log:
signal hilbert: correction to axis handling, new axis argument, with tests, ticket:1093

Modified: trunk/scipy/signal/signaltools.py
===================================================================
--- trunk/scipy/signal/signaltools.py	2010-01-16 10:29:43 UTC (rev 6204)
+++ trunk/scipy/signal/signaltools.py	2010-01-18 04:56:46 UTC (rev 6205)
@@ -1028,22 +1028,24 @@
     return w
 
 
-def hilbert(x, N=None):
+def hilbert(x, N=None, axis=-1):
     """Compute the analytic signal.
 
-    The transformation is done along the first axis.
+    The transformation is done along the last axis by default.
 
     Parameters
     ----------
     x : array-like
         Signal data
     N : int, optional
-        Number of Fourier components. Default: ``x.shape[0]``
+        Number of Fourier components. Default: ``x.shape[axis]``
+    axis : int, optional
+        
 
     Returns
     -------
-    xa : ndarray, shape (N,) + x.shape[1:]
-        Analytic signal of `x`
+    xa : ndarray
+        Analytic signal of `x`, of each 1d array along axis
 
     Notes
     -----
@@ -1054,6 +1056,8 @@
     where ``F`` is the Fourier transform, ``U`` the unit step function,
     and ``y`` the Hilbert transform of ``x``. [1]
 
+    changes in scipy 0.8.0: new axis argument, new default axis=-1
+
     References
     ----------
     .. [1] Wikipedia, "Analytic signal".
@@ -1062,13 +1066,13 @@
     """
     x = asarray(x)
     if N is None:
-        N = len(x)
+        N = x.shape[axis]
     if N <=0:
         raise ValueError, "N must be positive."
     if iscomplexobj(x):
         print "Warning: imaginary part of x ignored."
         x = real(x)
-    Xf = fft(x,N,axis=0)
+    Xf = fft(x, N, axis=axis)
     h = zeros(N)
     if N % 2 == 0:
         h[0] = h[N/2] = 1
@@ -1078,8 +1082,10 @@
         h[1:(N+1)/2] = 2
 
     if len(x.shape) > 1:
-        h = h[:, newaxis]
-    x = ifft(Xf*h)
+        ind = [newaxis]*x.ndim
+        ind[axis] = slice(None)
+        h = h[ind]
+    x = ifft(Xf*h, axis=axis)
     return x
 
 def hilbert2(x,N=None):

Modified: trunk/scipy/signal/tests/test_signaltools.py
===================================================================
--- trunk/scipy/signal/tests/test_signaltools.py	2010-01-16 10:29:43 UTC (rev 6204)
+++ trunk/scipy/signal/tests/test_signaltools.py	2010-01-18 04:56:46 UTC (rev 6205)
@@ -5,7 +5,7 @@
 from numpy.testing import *
 
 import scipy.signal as signal
-from scipy.signal import lfilter, correlate, convolve, convolve2d
+from scipy.signal import lfilter, correlate, convolve, convolve2d, hilbert
 
 
 from numpy import array, arange
@@ -704,5 +704,85 @@
         x = np.arange(6)
         assert_array_equal(signal.decimate(x, 2, n=1).round(), x[::2])
 
+
+class TestHilbert:
+    def test_hilbert_theoretical(self):
+        #test cases by Ariel Rokem
+        decimal = 14
+
+        pi = np.pi
+        t = np.arange(0, 2*pi, pi/256)
+        a0 = np.sin(t)
+        a1 = np.cos(t)
+        a2 = np.sin(2*t)
+        a3 = np.cos(2*t)
+        a = np.vstack([a0,a1,a2,a3])
+
+        h = hilbert(a)
+        h_abs = np.abs(h)
+        h_angle = np.angle(h)
+        h_real = np.real(h)
+
+        #The real part should be equal to the original signals:
+        assert_almost_equal(h_real, a, decimal)
+        #The absolute value should be one everywhere, for this input:
+        assert_almost_equal(h_abs, np.ones(a.shape), decimal)
+        #For the 'slow' sine - the phase should go from -pi/2 to pi/2 in
+        #the first 256 bins: 
+        assert_almost_equal(h_angle[0,:256], np.arange(-pi/2,pi/2,pi/256),
+                            decimal)
+        #For the 'slow' cosine - the phase should go from 0 to pi in the
+        #same interval:
+        assert_almost_equal(h_angle[1,:256], np.arange(0,pi,pi/256), decimal)
+        #The 'fast' sine should make this phase transition in half the time:
+        assert_almost_equal(h_angle[2,:128], np.arange(-pi/2,pi/2,pi/128),
+                            decimal)
+        #Ditto for the 'fast' cosine:
+        assert_almost_equal(h_angle[3,:128], np.arange(0,pi,pi/128), decimal)
+
+        #The imaginary part of hilbert(cos(t)) = sin(t) Wikipedia
+        assert_almost_equal(h[1].imag, a0, decimal)
+
+    def test_hilbert_axisN(self):
+        # tests for axis and N arguments
+        a = np.arange(18).reshape(3,6)
+        # test axis        
+        aa = hilbert(a, axis=-1)
+        yield assert_equal, hilbert(a.T, axis=0), aa.T
+        # test 1d
+        yield assert_equal, hilbert(a[0]), aa[0]
+        
+        # test N
+        aan = hilbert(a, N=20, axis=-1)
+        yield assert_equal, aan.shape, [3,20]
+        yield assert_equal, hilbert(a.T, N=20, axis=0).shape, [20,3]
+        #the next test is just a regression test,
+        #no idea whether numbers make sense
+        a0hilb = np.array(
+                          [  0.000000000000000e+00-1.72015830311905j ,
+                             1.000000000000000e+00-2.047794505137069j,
+                             1.999999999999999e+00-2.244055555687583j,
+                             3.000000000000000e+00-1.262750302935009j,
+                             4.000000000000000e+00-1.066489252384493j,
+                             5.000000000000000e+00+2.918022706971047j,
+                             8.881784197001253e-17+3.845658908989067j,
+                            -9.444121133484362e-17+0.985044202202061j,
+                            -1.776356839400251e-16+1.332257797702019j,
+                            -3.996802888650564e-16+0.501905089898885j,
+                             1.332267629550188e-16+0.668696078880782j,
+                            -1.192678053963799e-16+0.235487067862679j,
+                            -1.776356839400251e-16+0.286439612812121j,
+                             3.108624468950438e-16+0.031676888064907j,
+                             1.332267629550188e-16-0.019275656884536j,
+                            -2.360035624836702e-16-0.1652588660287j  ,
+                             0.000000000000000e+00-0.332049855010597j,
+                             3.552713678800501e-16-0.403810179797771j,
+                             8.881784197001253e-17-0.751023775297729j,
+                             9.444121133484362e-17-0.79252210110103j ])
+        yield assert_almost_equal, aan[0], a0hilb, 14, 'N regression'
+
+        
+
+
 if __name__ == "__main__":
     run_module_suite()



More information about the Scipy-svn mailing list