Implementing discrete and continuous Hilbert transforms

Question

What is an efficient and accurate Mathematica implementation of the Hilbert transform, for both continuous and especially discretely sampled functions?

This transform relates phase and amplitude in minimum phase systems.

Answer 1

Here's a direct implementation of the formula

$$\mathcal H(u)(t) = \frac1{\pi} -\hspace{-1.1em}\int_{-\infty}^\infty \frac{u(\tau)}{t-\tau}\, \mathrm d\tau$$

hilbertTransform[f_, u_, t_] :=
       FullSimplify[Convolve[f, 1/u, u, t, PrincipalValue -> True]/π]

Try it out:

hilbertTransform[#, v, w] & /@ {Sin[v], Cos[v], 1/(1 + v^2), Sinc[v], DiracDelta[v]}
   {-Cos[w], Sin[w], w/(1 + w^2), (1 - Cos[w])/w, 1/(π w)}

---

For the discrete Hilbert transform, here is a Mathematica routine:

hilbert[data_?VectorQ] := Module[{fopts = FourierParameters -> {1, -1}, e, n},
   e = Boole[EvenQ[n = Length[data]]]; 
   Im[InverseFourier[Fourier[data, fopts] * 
                     PadRight[ArrayPad[ConstantArray[2, Quotient[n, 2] - e], {1, e}, 1], n],
                     fopts]]] /; And @@ Thread[Im[data] == 0]

(making everything completely analogous to FourierTransform[] and Fourier[]). The algorithm is based on the routine in Marple's paper, and is essentially the same algorithm used by the function hilbert() in MATLAB's Signal Processing Toolbox.

Examples:

hilbert[{1, -2, 1}]
   {1.73205, 0., -1.73205}

hilbert[{1, -2, 1, 2}]
   {2., 0., -2., 0.}

Answer 2

You can realize a discrete Hilbert transform by convolving your discrete signal with a Hilbert kernel. The convolution is implemented with least effort in the frequency domain, where the spectrum of the Hilbert kernel is $$\sigma_H(\omega)=-i\cdot\mathrm{sgn}(\omega)$$ where $\omega$ is the angular frequency.

Continuous case

We define a function to perform the convolution in the frequency domain for us

HilbertTransform[f_,u_,t_] := Module[{fp = FourierParameters -> {1, -1}, x},
    InverseFourierTransform[
      -I (2 HeavisideTheta[x] - 1) FourierTransform[f, u, x, fp],
      x, t, fp
    ]
]

by exploiting the convolution theorem and getting the result back to our domain of interest via an inverse Fourier transform. In this case the HeavisideTheta function turns out to be better suited for modeling $\mathrm{sgn}(\omega)$ than the more intuitive Sign function, which can give wrong results in extreme cases like DiracDelta.

Some results for the test cases from J.M.'s answer:

In[1] := HilbertTransform[#, v, w] & /@
             {Sin[v], Cos[v], 1/(1 + v^2), Sinc[v], DiracDelta[v]}
Out[1] = {-Cos[w], Sin[w], w/(w^2+1), 1/w-Cos[w]/w, 1/(π w)}

Discrete case

In the discrete case we prepare the phase shifts for the occurring frequencies for a given data size n

HilbertSpectrum[0] := {};
HilbertSpectrum[n_Integer?Positive] := With[{nhalf = Quotient[n - 1, 2]},
  Join[{0}, ConstantArray[-I, nhalf], If[EvenQ[n], {0}, {}]
          , ConstantArray[ I, nhalf]
  ]
];

which we can then use in

Hilbert[data_?VectorQ, padding_Integer?NonNegative] := Module[
  {fp = FourierParameters -> {1, -1}, n = Length[data], m, paddeddata},
  m = n + padding;
  paddeddata = PadRight[data, m];
  Re @ InverseFourier[ HilbertSpectrum[m] Fourier[paddeddata, fp], fp][[;;n]]
]

which, like its continuous counterpart, performs the convolution in frequency space. To reduce artifacts due to circular convolution in case of nonperiodic signals we can choose to pad the data with zeros before convolution and cut to its original length after the transform. If we apply it to a simple sine wave

With[{n = 128}, sinedata = N@Table[Sin[2 π 3 k/n], {k, n}];]
     ListLinePlot[{sinedata, Hilbert[sinedata, 0]}]

this results in a good approximation of an inverted cosine:

But we can also try the transform on more interesting waveforms, for example a sine with continuously changing frequency

With[{n = 128}, sinedata = N@Table[Sin[2 π k/15 k/n], {k, n}];]
     ListLinePlot[{sinedata, Hilbert[sinedata, 256]}]

or a transition between different waveforms

With[{n = 128},
  data = N@Join[
    Table[TriangleWave[3 k/n], {k, n}], 
    Table[SquareWave[3 k/n], {k, n}]
  ];
]
ListLinePlot[{data, Hilbert[data, 256]}]

Answer 3

For continuous signal, I think is easy, took this in a course. For discrete, hard for me, we did not study it at school (yet).

But if the signal is continuous, this gives the Hilbert transform of the signal:

f[t_] := Sin[t];
g = FourierTransform[f[t], t, omega];
InverseFourierTransform[I* Sign[omega]*g, omega, t]

==>

-Cos[t]

P.S. I just saw a reference for the original paper for the implementation of discrete Hilbert transform. (link reference) by Kak, 1970.

note(1)

Just did a search on the net, and found what seems like a good reference with lots of Mathematica code for Hilbert transform, and a code for the discrete one. Reference is Handbook of geophysical exploration, volume 1 by Klaus Helbig and Sven Treitel. And on page 162, they give the implementation:

hilbert[RealTrace_] := (lh = Length[RealTrace];
   If[OddQ[lh], lh2 = Ceiling[lh/2], lh2 = lh/2];
   fd = Fourier[RealTrace, FourierParameters -> {1, -1}];
   Do[fd[[i]] = fd[[i]]*E^(+I*Pi/2), {i, lh2 + 1, lh}];
   Do[fd[[i]] = fd[[i]]*E^(-I*Pi/2), {i, 1, lh2}];
   QuadratureTrace = 
    Re[Chop[InverseFourier[fd, FourierParameters -> {1, -1}]]];);

Here is a link to the on-line version I saw the above (link here) if you like to read more about it.

note(2)

I wrote the function from the book above to make it easy for others to use, here it is below, and put a quick Manipulate around it, used the triangle and square functions shown in the other answer here to make it to compare.

Manipulate[Module[{data1, data2, k, opts},
  opts = {
    ImageSize -> {300}, ImagePadding -> 20, AxesLabel -> {"n", "ht(n)"}, 
    PlotLabel -> Text@Row[{Style["data", Blue], Spacer[10], 
    Style["Hilbert", Red]}], PlotStyle -> {Blue, Red}, 
    PlotMarkers -> {Automatic, Tiny}, Joined -> True
  };
  data1 = N@Table[TriangleWave[3 k/n], {k, n}];
  data2 = N@Table[SquareWave[3 k/n], {k, n}];

  Grid[{
    {ListPlot[{data1, hilbert[data1]}, opts]},
    {ListPlot[{data2, hilbert[data2]}, opts]}
    }, Frame -> All]
  ],    
 {{n, 32, "number of points"}, 32, 256, 1, Appearance -> "Labeled"},
 ContinuousAction -> False,
 Initialization :>
  (       
   (*function below is from Handbook of geophysical exploration,       
     volume 1 by Klaus Helbig and Sven Treitel, page 162*)

   hilbert[RealTrace_] := Module[{lh, fd, lh2, i},
     lh = Length[RealTrace];
     If[OddQ[lh], lh2 = Ceiling[lh/2], lh2 = lh/2];
     fd = Fourier[RealTrace, FourierParameters -> {1, -1}];
     Do[fd[[i]] = fd[[i]]*E^(I*Pi/2), {i, lh2 + 1, lh}];
     Do[fd[[i]] = fd[[i]]*E^(-I*Pi/2), {i, 1, lh2}];
     QuadratureTrace = 
      Re[Chop[InverseFourier[fd, FourierParameters -> {1, -1}]]]
     ]
   )
 ]

Answer 4

All the functions presented here make use of Fourier transforms to calculate the discrete Hilbert transform in frequency space. One can of course also compute the Hilbert transform directly in position space:

HilbertDirect = Compile[{{data, _Real, 1}},
   Table[
         2/\[Pi] 
         Sum[ 
              data[[i]]/(j - i)
         , {i, 1 + Mod[j, 2],Length[data] - Mod[j, 2], 2}]
   , {j, 1, Length[data]}]
]

Or, with slightly improved sensitivity for higher frequency components:

HilbertDirect2 = Compile[{{data, _Real, 1},{nn, _Integer}},
   Table[
        2/nn 
        Sum[
            data[[i]] Cot[((j - i) \[Pi])/nn]
        , {i,1 + Mod[j, 2], Length[data] - Mod[j, 2], 2}]
   , {j, 1, Length[data]}]
]

(where one should set nn>=Length[data]).

In most cases, this is of course slower than the Fourier transform approach, since the above scales directly with Length[data]. However, the quality of this finite impulse response Hilbert transform seems pretty good. Especially the edge effects for non-periodic cases are better than low padding frequency space results.