[SciPy-dev] Is this a bugfix for scipy.hilbert?
josef.pktd@gmai...
josef.pktd@gmai...
Sun Jan 17 23:54:05 CST 2010
On Sun, Jan 17, 2010 at 10:49 PM, <josef.pktd@gmail.com> wrote:
> 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,
committed in http://projects.scipy.org/scipy/changeset/6205
There is also a hilbert2 for 2d convolution. It is not in the
documentation (at least in not my old ones), and I didn't try whether
it is correct.
Josef
>
> 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