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

## 4 Answers

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

• @ 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! – user55777 Jul 1 at 10:49
• Yes, it seems the behavior of Convolve[] has changed in succeeding versions. I don't know how to fix that. – J. M. is away Jul 28 at 2:32

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 := HilbertTransform[#, v, w] & /@
{Sin[v], Cos[v], 1/(1 + v^2), Sinc[v], DiracDelta[v]}
Out = {-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 := {};
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]}] • 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]. – J. M. is away Jan 22 '12 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. – Thies Heidecke Jan 22 '12 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! – Maxim Sakharov Apr 19 '17 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}];
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,
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}]]]
]
)
]

• 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}]]] – J. M. is away Jan 21 '12 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! – user55777 Jul 1 at 9:31

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.

## protected by Kuba♦Apr 19 '17 at 12:26

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?