Re: Hilbert transform

[Sergei Steshenko]

----- 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.

>