Re: Hilbert transform

[John B. Thoo]

On Jul 6, 2012, at 3:34 PM, Sergei Steshenko 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.

Hi, everyone.  I'm sorry for joining late, and with no answers, but with a 
question.

As Sergei pointed out, with  x0 = [-3,-1,1,3],  H(H(x0)) \neq -x0.  However, 
with  y0 = [-3,-1,3,1],  indeed  H(H(y0)) = -y0.


octave-3.2.3:13> y0 = [-3,-1,3,1]
y0 =

  -3  -1   3   1

octave-3.2.3:14> y1 = imag (hilbert (y0))
y1 =

   1  -3  -1   3

octave-3.2.3:15> imag (hilbert (y1)) + y0
ans =

   0   0   0   0

octave-3.2.3:16> 


Why the difference?

Btw, here is another routine for the Hilbert transform that avoids using "imag".

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% begin hilb.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = hilb(in)
% HILBERT computes real hilbert transform of real input
% Has symbol -i sgn(k)
lx = length (in);
in = fft (in);
in (1:(lx/2+1)) = -i*in (1:(lx/2+1)); 
in (lx/2+2:lx) = i*in (lx/2+2:lx); 
out = real (ifft(in));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% begin hilb.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

---John.

-----------------------------------------------------------------------
"Ten thousand difficulties do not make one doubt....  A man may be annoyed that 
he cannot work out a mathematical problem ... without doubting that it admits 
an answer."

---John Henry Newman [_Apologia_, p. 239 in Project Gutenberg's 
   <http://www.gutenberg.org/ebooks/22088>]