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

Ariel Rokem arokem@berkeley....
Fri Jan 15 00:44:54 CST 2010


Yes - looks good. Except I would prefer to eventually set the axis to
default to -1, to be consistent with signal.fft (and also np.fft.fft) which
has axis=-1.

As for whether it's doing what it's supposed to do, for what it's worth - it
seems to do similar things to what Matlab's 'hilbert' function does on a few
simple examples I tried out.

Cheers,

Ariel



On Thu, Jan 14, 2010 at 8:53 PM, <josef.pktd@gmail.com> wrote:

> On Thu, Jan 14, 2010 at 11:27 PM,  <josef.pktd@gmail.com> wrote:
> > On Thu, Jan 14, 2010 at 11:02 PM,  <josef.pktd@gmail.com> wrote:
> >> On Thu, Jan 14, 2010 at 10:54 PM,  <josef.pktd@gmail.com> wrote:
> >>> On Thu, Jan 14, 2010 at 10:24 PM, Ariel Rokem <arokem@berkeley.edu>
> wrote:
> >>>> 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
> >>>
> >>> I think your change would break the currently advertised behavior,
> >>> axis=0 (The transformation is done along the first axis)
> >>>
> >>> but fft and ifft have default axis=-1
> >>>
> >>> fft in hilbert uses axis=0 as in docstring
> >>> but ifft uses default axis=-1
> >>>
> >>> so, I would think the fix should be  x = ifft(Xf*h, axis=0)
> >>>
> >>> But as it currently looks like the axis argument doesn't work anyway,
> >>> there wouldn't be much breakage if the axis would be included as an
> >>> argument and default to -1.
> >>> However, I don't know what the "standard" for scipy.signal is for
> default axis.
> >>>
> >>> Josef
> >>
> >> after adding axis to ifft:
> >>>>> print hilbert(aa).real
> >> [[ 0.82584851  0.15215031  0.14767381]
> >>  [ 0.95021675  0.16803995  0.43562964]
> >>  [ 0.13033881  0.06198952  0.70729614]
> >>  [ 0.69409563  0.06962778  0.72552601]
> >>  [ 0.34297612  0.50579001  0.86463304]
> >>  [ 0.28355261  0.21626889  0.85165102]
> >>  [ 0.49481491  0.21290645  0.71416814]
> >>  [ 0.2645843   0.95783096  0.77514016]
> >>  [ 0.38735994  0.14274852  0.56344808]
> >>  [ 0.88084015  0.39879649  0.64949951]]
> >>>>> print hilbert(aa[:,:1]).real
> >> [[ 0.82584851]
> >>  [ 0.95021675]
> >>  [ 0.13033881]
> >>  [ 0.69409563]
> >>  [ 0.34297612]
> >>  [ 0.28355261]
> >>  [ 0.49481491]
> >>  [ 0.2645843 ]
> >>  [ 0.38735994]
> >>  [ 0.88084015]]
> >>
> >> but it treats a 1d array as row vector and transforms along zero axis
> >> of length 1, and not along the length of the array.
> >> so another fix to handle 1d arrays correctly should be done
> >>
> >>>>> print hilbert(aa[:,1]).real
> >> [ 0.15215031  0.16803995  0.06198952  0.06962778  0.50579001  0.21626889
> >>  0.21290645  0.95783096  0.14274852  0.39879649]
> >>>>> aa[:,1]
> >> array([ 0.15215031,  0.16803995,  0.06198952,  0.06962778,  0.50579001,
> >>        0.21626889,  0.21290645,  0.95783096,  0.14274852,  0.39879649])
> >>>>>
> >
> > there's something wrong with my example, the real part is the same
> > which confused me
> >
> > it works correctly with 1d
> >
> >>>> np.abs(hilbert(aa[:,0]))
> > array([ 0.83251128,  1.04487091,  0.27702083,  0.69901499,  0.49170197,
> >        0.31227114,  0.49505637,  0.26461488,  0.61385196,  0.90716272])
> >
> >>>> np.abs(hilbert(aa[:,:1])).T
> > array([[ 0.83251128,  1.04487091,  0.27702083,  0.69901499,  0.49170197,
> >         0.31227114,  0.49505637,  0.26461488,  0.61385196,  0.90716272]])
> >
> >>>> np.abs(hilbert(aa))[:,0]
> > array([ 0.83251128,  1.04487091,  0.27702083,  0.69901499,  0.49170197,
> >        0.31227114,  0.49505637,  0.26461488,  0.61385196,  0.90716272])
> >
> > besides reading the docstring, I don't know what hilbert is supposed
> > to be good for.
>
> Would something like the function in the attachment do ?
>
>
>
> > Josef
> >
> >
> >> Josef
> >>
> >>
> >>>
> >>>>
> >>>>
> >>>> Cheers,
> >>>>
> >>>> Ariel
> >>>> --
> >>>> Ariel Rokem
> >>>> Helen Wills Neuroscience Institute
> >>>> University of California, Berkeley
> >>>> http://argentum.ucbso.berkeley.edu/ariel
> >>>>
> >>>> _______________________________________________
> >>>> 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
>
>


-- 
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.scipy.org/pipermail/scipy-dev/attachments/20100114/78e58b84/attachment-0001.html 


More information about the SciPy-Dev mailing list