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.
Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this communityWhat 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] *
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.}
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!
$\endgroup$
– user55777
Jul 1 '19 at 10:49
Convolve[]
has changed in succeeding versions. I don't know how to fix that.
$\endgroup$
– J. M.'s ennui♦
Jul 28 '19 at 2:32
hilbert
function on an InterpolatingFunction
, can you help with this question?
$\endgroup$
– jsxs
Feb 14 '20 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},
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]}]
HilbertSpectrum[]
can be compacted a bit: HilbertSpectrum[n_Integer?Positive] := PadRight[ArrayPad[ConstantArray[-I, Quotient[n - 1, 2]], {1, Boole[EvenQ[n]]}], n, I]
.
$\endgroup$
– J. M.'s ennui♦
Jan 22 '12 at 16:16
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.
$\endgroup$
– Thies Heidecke
Jan 22 '12 at 16:47
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}]]]
]
)
]
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}]]]
$\endgroup$
– J. M.'s ennui♦
Jan 21 '12 at 5:36
I
in InverseFourierTransform[I* Sign[omega]*g, omega, t]
? Thank you in advance!
$\endgroup$
– user55777
Jul 1 '19 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.