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

[josef.pktd at gmail.com]

On Thu, Jan 14, 2010 at 10:24 PM, Ariel Rokem <arokem at 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


>
>
> Cheers,
>
> Ariel
> --
> Ariel Rokem
> Helen Wills Neuroscience Institute
> University of California, Berkeley
> http://argentum.ucbso.berkeley.edu/ariel
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>