[SciPy-dev] Is this a bugfix for scipy.hilbert?
Ariel Rokem
arokem at berkeley.edu
Fri Jan 15 14:34:34 EST 2010
Hi -
attached is a file with a couple of tests. I am not sure this tests the
issues we were dealing with previously (the axis issues, etc.), but it has
some sensible test-cases, which compare to what Matlab would give you (not
quite 10by3 or 10by6, but as you can see, they make sense). Also - all the
assertions are assert_almost_equal. Do you think that's OK? I think there
are float-precision issues here, which would make assert_equal fail, but I
am not sure - I would be happy to get any general comments on these tests,
in case I am doing this all wrong.
Cheers,
Ariel
On Thu, Jan 14, 2010 at 11:10 PM, <josef.pktd at gmail.com> wrote:
> On Fri, Jan 15, 2010 at 1:44 AM, Ariel Rokem <arokem at berkeley.edu> wrote:
> > 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.
>
> I'm indifferent to the default axis, from a quick look and my
> experience there are not many functions with axis arguments in signal.
> So I'm fine with switching to axis=-1. We should do it with this
> bugfix, since until now the function wasn't correct anyway for 2d.
>
> >
> > 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.
>
> I was reading briefly on wikipedia, and checked with fftpack.hilbert,
> which returns the same array as signal.hilbert(a).imag, but I didn't
> manage to figure out why fftpack.hilbert only allows 1d (i got lost
> starting at convolve.pyf)
>
> Could you write a simple test case compared to matlab, e.g. 10by3 as
> in my example, for both axis, or 10by6 if 10by3 doesn't make sense?
>
> If nobody objects, I can commit the change with axis=-1.
>
> Josef
>
> >
> > Cheers,
> >
> > Ariel
> >
> >
> >
> > On Thu, Jan 14, 2010 at 8:53 PM, <josef.pktd at gmail.com> wrote:
> >>
> >> On Thu, Jan 14, 2010 at 11:27 PM, <josef.pktd at gmail.com> wrote:
> >> > On Thu, Jan 14, 2010 at 11:02 PM, <josef.pktd at gmail.com> wrote:
> >> >> On Thu, Jan 14, 2010 at 10:54 PM, <josef.pktd at gmail.com> wrote:
> >> >>> 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
> >> >>
> >> >> 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 at scipy.org
> >> >>>> http://mail.scipy.org/mailman/listinfo/scipy-dev
> >> >>>>
> >> >>>>
> >> >>>
> >> >>
> >> >
> >>
> >> _______________________________________________
> >> SciPy-Dev mailing list
> >> SciPy-Dev at 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
> >
> > _______________________________________________
> > SciPy-Dev mailing list
> > SciPy-Dev at scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-dev
> >
> >
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at 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.python.org/pipermail/scipy-dev/attachments/20100115/b55b63c1/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test_hilbert.py
Type: application/octet-stream
Size: 1028 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20100115/b55b63c1/attachment.obj>
More information about the SciPy-Dev
mailing list