Re: Hilbert transform

[Ozzy Lash]

On Fri, Jul 6, 2012 at 6:03 PM, Sergei Steshenko <address@hidden> 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.

If you make the length of the vector odd, then you see a very small
error. For example:

octave:13> foo = center(linspace(0,1,101));
octave:14> max(abs(foo + imag(hilbert(imag(hilbert(foo))))))
ans =  5.5511e-16
octave:15> min(abs(foo + imag(hilbert(imag(hilbert(foo))))))
ans = 0

I'm not sure if this is expected, or is caused by the implementation.

Bill