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

josef.pktd@gmai... josef.pktd@gmai...
Sun Jan 17 21:49:18 CST 2010


On Sun, Jan 17, 2010 at 12:31 AM,  <josef.pktd@gmail.com> wrote:
> On Fri, Jan 15, 2010 at 4:20 PM, Ariel Rokem <arokem@berkeley.edu> wrote:
>> Hi - I've never done this before, so it would be great if I could 'look over
>> your shoulder' (in the sense that I know how this ticket came about :D), as
>> you submit a ticket on this.
>
> Done in http://projects.scipy.org/scipy/ticket/1093
> tests pass at 14 decimals


While adding some tests, I got one more question

According to Wikipedia http://en.wikipedia.org/wiki/Analytic_signal#Definition
the imaginary part of the analytical signal is equal to the Hilbert transform

however, fftpack.hilbert has the opposite sign.

>From the examples signal.hilbert looks correct, so does
fftpack.hilbert have a sign mistake or is it based on a different
definition?


>>> r = np.random.randn(20)
>>> fftpack.hilbert(r)
array([-0.27285468, -1.39747965,  1.7991044 , -0.16609304, -1.84459577,
        0.48696479, -0.33190553,  0.59383033,  2.15361055, -0.89341275,
       -0.13730369,  0.84046658,  1.38110384, -1.7595949 , -0.04869402,
        0.59871558, -1.09627219,  0.59375139, -1.6021929 ,  1.10285168])
>>> hilbert(r).imag
array([ 0.27285468,  1.39747965, -1.7991044 ,  0.16609304,  1.84459577,
       -0.48696479,  0.33190553, -0.59383033, -2.15361055,  0.89341275,
        0.13730369, -0.84046658, -1.38110384,  1.7595949 ,  0.04869402,
       -0.59871558,  1.09627219, -0.59375139,  1.6021929 , -1.10285168])


I'm just checking definitions for the tests,

Josef
>
> Josef
>
>>
>> Thanks --
>>
>> Ariel
>>
>> On Fri, Jan 15, 2010 at 11:48 AM, <josef.pktd@gmail.com> wrote:
>>>
>>> On Fri, Jan 15, 2010 at 2:34 PM, Ariel Rokem <arokem@berkeley.edu> wrote:
>>> > 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.
>>>
>>> nice test cases, I like theoretical tests even better than verified
>>> numbers from other packages.
>>>
>>> Besides some cosmetic changes to get them into a test function, the
>>> only part to add is the precision of the tests.
>>> The default precision of assert_almost_equal is only 6 decimals.
>>>
>>> For these kind of cases, I usually go to 12 to 15 depending on the
>>> numerical precision of the algorithm. Usually, I go by trial and error
>>> until the test breaks, or calculate max abs of the error.
>>>
>>> I can add some tests for the axis argument.
>>>
>>> Can you open a ticket for the record or shall I ?
>>>
>>> Josef
>>>
>>>
>>> >
>>> > Cheers,
>>> >
>>> > Ariel
>>> >
>>> > On Thu, Jan 14, 2010 at 11:10 PM, <josef.pktd@gmail.com> wrote:
>>> >>
>>> >> On Fri, Jan 15, 2010 at 1:44 AM, Ariel Rokem <arokem@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@gmail.com> wrote:
>>> >> >>
>>> >> >> 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
>>> >> >> >>>>
>>> >> >> >>>>
>>> >> >> >>>
>>> >> >> >>
>>> >> >> >
>>> >> >>
>>> >> >> _______________________________________________
>>> >> >> SciPy-Dev mailing list
>>> >> >> SciPy-Dev@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@scipy.org
>>> >> > http://mail.scipy.org/mailman/listinfo/scipy-dev
>>> >> >
>>> >> >
>>> >> _______________________________________________
>>> >> SciPy-Dev mailing list
>>> >> SciPy-Dev@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@scipy.org
>>> > http://mail.scipy.org/mailman/listinfo/scipy-dev
>>> >
>>> >
>>> _______________________________________________
>>> SciPy-Dev mailing list
>>> SciPy-Dev@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@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>
>>
>


More information about the SciPy-Dev mailing list