Re: Hilbert transform

[Ben Abbott]

On Jul 6, 2012, at 7:03 PM, Sergei Steshenko wrote:

> 
> ----- Original Message -----
>> From: Juan Pablo Carbajal <address@hidden>
>> To: Sergei Steshenko <address@hidden>
>> Cc: Ben Abbott <address@hidden>; "address@hidden" <address@hidden>
>> Sent: Saturday, July 7, 2012 1:49 AM
>> Subject: Re: Hilbert transform
>> 
>> On Sat, Jul 7, 2012 at 12:34 AM, Sergei Steshenko <address@hidden> 
>> wrote:
>>> 
>>> ----- Original Message -----
>>>> From: Ben Abbott <address@hidden>
>>>> To: Sergei Steshenko <address@hidden>
>>>> Cc: "address@hidden" <address@hidden>
>>>> Sent: Friday, July 6, 2012 8:53 PM
>>>> Subject: Re: Hilbert transform
>>>> 
>>>> 
>>>> On Jul 6, 2012, at 12:44 PM, Sergei Steshenko wrote:
>>>> 
>>>>>   Hello,
>>>>> 
>>>>>   i am talking about 'hilbert' function from from
>>>> 'signal-1.1.3/hilbert.m' file, so Octave help list purists are 
>> welcome
>>>> to send me with my uncomfortable questions to octave-dev list.
>>>>> 
>>>>>   But I'll ask my questions here - from my reading (and 
>> recollections of
>>>> what I learned a long long time ago) the issue is 
>> mathematical/computational.
>>>>> 
>>>>>   First a couple of references:
>>>>> 
>>>>>   1) http://w3.msi.vxu.se/exarb/mj_ex.pdf - I think put together 
>> really
>>>> nicely;
>>>>>   2) http://en.wikipedia.org/wiki/Hilbert_transform
>>>>>   .
>>>>> 
>>>>>   Wherever we look, we find that the definition of Hilbert transform 
>> is
>>>> through integral of a _real_ function, i.e.
>>>>> 
>>>>>   hilbert(u(t)) == integral_from_minus_to_plus_inf("u(tau) / (t 
>> -
>>>> tau)", "dtau")
>>>>> 
>>>>>   and as such it should be a _real_ function of 't' provided 
>> u(t) is
>>>> a real function of 't'.
>>>>> 
>>>>> 
>>>>>   Also, it is proven that
>>>>> 
>>>>>   hilbert(hilbert(u(t))) == -u(t)
>>>>>   .
>>>>> 
>>>>>   Now, here is Octave and its package reality:
>>>>> 
>>>>> 
>>>>>   "
>>>>>   octave:1> hilbert([1 2 3 4])
>>>>>   ans =
>>>>> 
>>>>>      1 + 1i   2 - 1i   3 - 1i   4 + 1i
>>>>> 
>>>>>   octave:2> hilbert(hilbert([1 2 3 4]))
>>>>>   warning: HILBERT: ignoring imaginary part of signal
>>>>>   ans =
>>>>> 
>>>>>      1 + 1i   2 - 1i   3 - 1i   4 + 1i
>>>>> 
>>>>>   octave:3>
>>>>>   ".
>>>>> 
>>>>>   Three violations already:
>>>>> 
>>>>>   1) output is complex rather than real;
>>>>>   2) the transform is not invertible;
>>>>>   3) since Hilbert transform is linear, complex input should be 
>> accepted
>>>> according to
>>>>> 
>>>>>   hilbert(foo + i * bar) == hilbert(foo) + i * hilbert(bar)
>>>>>   .
>>>>> 
>>>>>   To put things politically correctly, Hilbert transform is a canine 
>> female
>>>> to calculate - because of the above "/ (t -tau)", and 
>> it's
>>>> problematic to calculate in discrete domain.
>>>>> 
>>>>>   So, my first practical question is: "What does Matlab do 
>> ?".
>>>>> 
>>>>>   Thanks,
>>>>>     Sergei.
>>>> 
>>>> Matlab's online doc for hilbert() is at the link below.
>>>> 
>>>> http://www.mathworks.com/help/toolbox/signal/ref/hilbert.html
>>>> 
>>>> The example below is included on that page.
>>>> 
>>>>      hilbert ([1 2 3 4])
>>>>      ans =   1 + 1i   2 - 1i   3 - 1i   4 + 1i
>>>> 
>>>> It appears that the hilbert() function is *not* a direct implementation 
>> of the
>>>> Hilbert Transform, but the version in the signal package does appear to 
>> be
>>>> consistent with the one which is part of the Matlab Signals toolbox.
>>>> 
>>>> After a quick look, my impression is that the hilbert() function's 
>> output is
>>>> hilbert(x) = 1i * H(x) + x.  Where, H(x) is the Hilbert Transform.  I 
>> don' t
>>>> know why the author (mathworks?) decided to restrict the input to real 
>> values.
>>>> 
>>>> Ben
>>> 
>>> Thanks to all for clarifications and explanations.
>>> 
>>> I still see a problem though. First, as others have pointed out, what the 
>> documentation of signal-1.1.3/hilbert.m says among other things:
>>> 
>>> "
>>> `real(H)' contains the original signal F.  `imag(H)' contains the
>>>       Hilbert transform of F.
>>> ".
>>> 
>>> So, if I want Hilbert transform proper, I need 'imag' - so far so 
>> good.
>>> 
>>> Here is a quick example:
>>> 
>>> "
>>> octave:8> imag(hilbert(imag(hilbert([-3 -1 1 3]))))
>>> ans =
>>> 
>>>     2   2  -2  -2
>>> ".
>>> 
>>> The input ('[-3 -1 1 3]') is a 0 DC vector:
>>> 
>>> "
>>> center([-3 -1 1 3])
>>> ans =
>>> 
>>>    -3  -1   1   3
>>> "
>>> 
>>> - which makes life easier.
>>> 
>>> So, above I expected '3 1 -1 -3' answer according to
>>> 
>>> Hilbert(Hilbert(u(t))) == -u(t)
>>> 
>>> property ('Hilbert' is meant to be true Hilbert transform).
>>> 
>>> 
>>> Obviously the answer I got: "2   2  -2  -2" is not the expected 
>> one. And the difference is not just scaling coefficient.
>>> 
>>> 
>>> Any ideas ?
>>> 
>>> Thanks,
>>>    Sergei.
>>> 
>>> P.S. It looks like signal-1.1.3/hilbert.m follows what Matlab documentation 
>> describes, but the question is whether what is described in the 
>> documentation 
>> (the part returned with 'imag') is true Hilbert transform.
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>>> 
>>> _______________________________________________
>>> Help-octave mailing list
>>> address@hidden
>>> https://mailman.cae.wisc.edu/listinfo/help-octave
>> 
>> You are working with a discrete transformation very sensitive to the
>> length of the signal. If you check with a longer signal everything is
>> fine.
>> 
>> t  = linspace (0,1,100);
>> y  = sin (2*pi*8*t) - sin (2*pi*10*t);
>> yg= imag ( hilbert ( imag ( hilbert (y))));
>> plot (t,y,t,-yg,'o')
>> 
>> Cheers
>> 
>> 
>> -- 
>> M. Sc. Juan Pablo Carbajal
>> -----
>> PhD Student
>> University of Zürich
>> http://ailab.ifi.uzh.ch/carbajal/
> 
> Thanks for the example.
> 
> I'd say with long sequences it's much better, but still:
> 
> "
> octave:37> foo = center(linspace(0,1,100));
> octave:38> max(abs(foo + imag(hilbert(imag(hilbert(foo))))))
> ans =  0.0050505
> "
> 
> - the expected answer is 0.
> 
> And, even more interesting:
> 
> "
> octave:39> min(abs(foo + imag(hilbert(imag(hilbert(foo))))))
> ans =  0.0050505
> ",
> 
> i.e. the error is constant.
> 
> 
> Thanks,
>   Sergei.

The "error" is mostly due to adding "foo"

        foo = center(linspace(0,1,100));
        min (abs (foo))
        ans =  0.0050505

Ben