Implementing discrete and continuous Hilbert transforms

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.

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] *
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.}

• @ J.M.is away when I try your hilbertTransform with the first 3 tests, all give me the right answers as yours, but for the last one, hilbertTransform[DiracDelta[v], v, w] gives Convolve[DiracDelta[v], 1/v, v, w, PrincipalValue -> True]/\[Pi], which is unevaluated and different your output. How to fix it? I am based on v11.3.0.0 and win7. Thank you! Commented Jul 1, 2019 at 10:49
• Yes, it seems the behavior of Convolve[] has changed in succeeding versions. I don't know how to fix that. Commented Jul 28, 2019 at 2:32
• @J.M.isinlimbo how to apply the discrete hilbert function on an InterpolatingFunction, can you help with this question?
– lxy
Commented Feb 14, 2020 at 7:45

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},
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]}]


• Your implementation of the positive integer case of HilbertSpectrum[] can be compacted a bit: HilbertSpectrum[n_Integer?Positive] := PadRight[ArrayPad[ConstantArray[-I, Quotient[n - 1, 2]], {1, Boole[EvenQ[n]]}], n, I]. Commented Jan 22, 2012 at 16:16
• Thanks for the suggestion! I like the use of ArrayPad with Boole and I haven't thought of using the more canonical Quotient. I'll adapt the use of Quotient. I think I'll still mainly stick with my current version because it's faster and I think it's a bit easier to see what the code is doing. Commented Jan 22, 2012 at 16:47
• Recently I've been trying to modify Hilbert[data_?VectorQ, padding_Integer?NonNegative] function to work with 2D data. While it's clear how to apply Fourier to matrices I'm lost with HilbertSpectrum. How does it look like for 2D data? Thanks in advance! Commented Apr 19, 2017 at 12:17

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}];
Re[Chop[InverseFourier[fd, FourierParameters -> {1, -1}]]];);


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}];
Re[Chop[InverseFourier[fd, FourierParameters -> {1, -1}]]]
]
)
]

• A little improvement: hilbert[RealTrace_] := Block[{fd, i}, fd = I Fourier[RealTrace, FourierParameters -> {1, -1}]; Do[fd[[i]] *= -1, {i, Ceiling[Length[RealTrace]/2]}]; Re[InverseFourier[fd, FourierParameters -> {1, -1}]]] Commented Jan 21, 2012 at 5:36
• @Nasser Given the property of Hilbert transform: $\mathcal{F}[\mathcal{H}[f]](\omega)=-i sign(\omega)\mathcal{F}[f](\omega)$. So we can take inverse Fourier on both sides to get the Hilbert transformation. But why there is not minus sign before I in InverseFourierTransform[I* Sign[omega]*g, omega, t]? Thank you in advance! Commented Jul 1, 2019 at 9:31
• @Nasser I checked this and tried it and the command f[t_] := Sin[t]; g = FourierTransform[f[t], t, omega]; InverseFourierTransform[I* Sign[omega]*g, omega, t] is exactly the same as FourierTransform[f[t], t, omega]. There is no difference between them. How can Hilbert transform be the same as Fourier transform? Commented Jan 16, 2023 at 10:28
• @Vangsnes Could you show screen shot? in V 13.2 FourierTransform[f[t], t, omega] is not the same as InverseFourierTransform[I*Sign[omega]*g, omega, t] ? Or may be I did not understand you. Here is screen shot !Mathematica graphics Commented Jan 17, 2023 at 0:12
• @Vangsnes Yes that is correct, InverseFourierTransform[I* Sign[omega]*g, omega, t] is the Hilbert transform. Thanks. Commented Jan 17, 2023 at 9:20

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.