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

josef.pktd@gmai... josef.pktd@gmai...
Thu Jan 14 22:02:38 CST 2010


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])
>>>

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
>>
>>
>


More information about the SciPy-Dev mailing list