[SciPy-dev] Is this a bugfix for scipy.hilbert?

[Ariel Rokem]

Hi everyone,

I have been trying to use scipy.signal.hilbert and I got the following
puzzling result:

In [22]: import scipy

In [23]: scipy.__version__ #I have r6182
Out[23]: '0.8.0.dev'

In [24]: import scipy.signal as signal

In [25]: a = np.random.rand(100,100)

In [26]: np.abs(signal.hilbert(a[-1]))
Out[26]:
array([ 0.57567681,  0.25918624,  0.50207097,  0.51834052,  0.24293389,
        0.5779464 ,  0.6515758 ,  0.89973173,  1.00275444,  0.37352935,
        0.62332717,  0.93599749,  0.40651376,  0.65088756,  0.8332281 ,
        0.5770101 ,  0.9288512 ,  0.46671906,  0.41536055,  0.71418068,
        0.81250913,  0.07652627,  0.72939072,  0.26755626,  0.36396146,
        0.59725999,  1.02264694,  0.41227986,  0.98122853,  0.71906675,
        0.58582611,  0.77288117,  0.3217015 ,  0.65261394,  0.11947618,
        0.75632703,  0.43432935,  0.52182485,  1.0277177 ,  1.01104986,
        0.3023265 ,  0.6024772 ,  0.69257548,  0.55418735,  0.46259052,
        0.25832231,  0.38278355,  0.45508532,  0.26215872,  0.34207947,
        0.80704729,  0.80755477,  0.95317178,  0.97458885,  0.58762294,
        0.82540618,  0.62005585,  0.82494646,  1.04221293,  0.14983027,
        1.01571579,  0.99381328,  0.24158714,  0.84256569,  0.53418924,
        0.24067628,  0.90489883,  1.02217747,  0.34988034,  0.5310065 ,
        0.48135002,  1.03020269,  0.6013679 ,  0.46062485,  0.3918485 ,
        0.21554545,  0.31704519,  0.04868385,  0.1787766 ,  0.37361852,
        0.21977912,  0.7649772 ,  0.77867281,  0.37684278,  0.64432638,
        0.77494951,  0.87106309,  0.77611484,  0.52666801,  0.88683667,
        0.69164967,  0.98618191,  0.84811375,  0.35934198,  0.32650478,
        0.1752677 ,  0.60574454,  0.5109132 ,  0.52332287,  0.99777805])

In [27]: np.abs(signal.hilbert(a))[-1]
Out[27]:
array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])


----------------------------------------------------------------------

I was expecting both of these to have the same values - am I missing
something?

I think that the following solves this issue, but now I am not that sure
whether it does what it is supposed to do and I couldn't find a test for
this in test_signaltools.py. Does anyone know of a good test-case for the
analytic signal, that I could create for this?

Index: scipy/signal/signaltools.py
===================================================================
--- scipy/signal/signaltools.py    (revision 6182)
+++ scipy/signal/signaltools.py    (working copy)
@@ -1062,13 +1062,13 @@
     """
     x = asarray(x)
     if N is None:
-        N = len(x)
+        N = x.shape[-1]
     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=-1)
     h = zeros(N)
     if N % 2 == 0:
         h[0] = h[N/2] = 1
@@ -1078,7 +1078,7 @@
         h[1:(N+1)/2] = 2

     if len(x.shape) > 1:
-        h = h[:, newaxis]
+        h = h[newaxis,:]
     x = ifft(Xf*h)
     return x



Cheers,

Ariel
-- 
Ariel Rokem
Helen Wills Neuroscience Institute
University of California, Berkeley
http://argentum.ucbso.berkeley.edu/ariel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20100114/f260a190/attachment.html>