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

josef.pktd@gmai... josef.pktd@gmai...
Thu Jan 14 22:53:19 CST 2010


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
>>>>
>>>>
>>>
>>
>
-------------- next part --------------
# -*- coding: utf-8 -*-
""" scipy.signal.hilbert bugfix

Created on Thu Jan 14 22:54:53 2010
changes by josef-pktd
"""
import numpy as np
from numpy import asarray, zeros, newaxis, real
from scipy import fftpack


def hilbert(x, N=None, axis=0):
    """Compute the analytic signal.

    The transformation is done along the first axis.

    Parameters
    ----------
    x : array-like
        Signal data
    N : int, optional
        Number of Fourier components. Default: ``x.shape[0]``
    axis : int, optional

    Returns
    -------
    xa : ndarray, shape (N,) + x.shape[1:]
        Analytic signal of `x`

    Notes
    -----
    The analytic signal `x_a(t)` of `x(t)` is::

        x_a = F^{-1}(F(x) 2U) = x + i y

    where ``F`` is the Fourier transform, ``U`` the unit step function,
    and ``y`` the Hilbert transform of ``x``. [1]

    References
    ----------
    .. [1] Wikipedia, "Analytic signal".
           http://en.wikipedia.org/wiki/Analytic_signal

    """
    x = asarray(x)
    if N is None:
        N = x.shape[axis]
    if N <=0:
        raise ValueError, "N must be positive."
    if np.iscomplexobj(x):
        print "Warning: imaginary part of x ignored."
        x = real(x)
    Xf = fftpack.fft(x,N,axis=axis)
    h = zeros(N)
    if N % 2 == 0:
        h[0] = h[N/2] = 1
        h[1:N/2] = 2
    else:
        h[0] = 1
        h[1:(N+1)/2] = 2

    if len(x.shape) > 1:
        ind = [newaxis]*x.ndim
        ind[axis] = slice(None)
        h = h[ind]
    x = fftpack.ifft(Xf*h, axis=axis)
    return x

aa = np.random.randn(10,3)
print np.abs(hilbert(aa))[:,0]
print np.abs(hilbert(aa[:,:1])).T
print np.abs(hilbert(aa,axis=1))[0]
print np.abs(hilbert(aa[0]))


More information about the SciPy-Dev mailing list