[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