# What is the algorithm to do a Discrete Hilbert Transform?

I have a time series data. I have applied Empirical Mode Decomposition and obtained several IMFs.

Now I need to apply a Hilbert Transform to the IMFs to get the instantaneous frequencies.

I have found two papers explaining how to do that:

this one, that says this:

What Dr. Huang proposed for more accurately calculate the instantaneous frequency is to do normalization of IMF. The process is as followed:

1. Take absolute value of IMF.
2. Find extrema.
3. Based on these extrema, construct envelope.
4. Normalize IMF using the envelope. The FM part of signal becomes almost equal amplitude.
5. Repeat process 2-4 after the amplitude of normalized IMF retains a straight line with identical value.
6. Find the instantaneous frequency on the normalized IMF.

and this one that says this:

Notice that the first one says to take the absolute value of IMF and the other one doesn't.

I have tried both but repeating the loop to make all values be contained between 0 and 1 is doing nothing. Every loop the values are the same.

How do I really do that and how do I find the frequency of each IMF?

The Discrete Hilbert Transform is an "ideal" (implying that this is not the practical implementation) linear time-invariant filter with input $$x[n]$$ having output

$$\mathscr{H}\Big\{ x[n] \Big\} \triangleq \hat{x}[n] = \sum\limits_{i=-\infty}^{\infty} h[i] \, x[n-i]$$

Where the impulse response of a discrete-time Hilbert transformer is:

$$h[n] = \begin{cases} \frac{\big(1 - (-1)^n\big)}{\pi n} \quad & n \ne 0 \\ \\ 0 & n = 0 \end{cases}$$

Because $$1 - (-1)^n = 0$$ for even $$n$$, this can be restated as

$$h[n] = \begin{cases} \frac{2}{\pi n} \quad & n \text{ odd} \\ \\ 0 & n \text{ even} \end{cases}$$

This is not a causal impulse response, nor is it finite in length. To make it finite in length, you would need to window it with a decent window:

$$h[n] = \begin{cases} \frac{2}{\pi n} w[n] \quad & n \text{ odd} \\ \\ 0 & n \text{ even} \end{cases}$$

where $$w[n]$$ is some window function of width $$L+1$$ samples (and $$L$$ is an even positive integer). If it were a Hamming Window, it would be:

$$w[n] = \begin{cases} 0.54 \ + \ 0.46 \cdot \cos\left(\pi \frac{n}{L/2} \right) \qquad & |n| \le L/2 \\ 0 & |n| > L/2 \\ \end{cases}$$

If it were a Kaiser window it would be

$$w[n] = \begin{cases} \frac{1}{I_0(\beta)} \, I_0\left(\beta \sqrt{1 - \left(\frac{n}{L/2}\right)^2 } \right) \qquad & |n| \le L/2 \\ 0 & |n| > L/2 \\ \end{cases}$$

where

$$I_0(u) \triangleq \sum\limits_{k=0}^{\infty} \frac{(-1)^k \big( \tfrac{u}{2} \big)^{2k}}{(k!)^2}$$

$$I_0(x)$$ is the 0th-order modified Bessel function of the 1st kind. $$L+1$$ is the number of non-zero samples or FIR taps (the FIR filter order is $$L$$) and, in my centered and symmetrical case, must be even. $$\beta$$ is a "shape parameter", maybe around 5 or 6, i dunno.

Now, to make this causal, your impulse response has to be delayed by $$\frac{L}2$$ samples to be $$h[n-\frac{L}2]$$, but then should also all other signals that this Hilbert output is compared to, they should also be delayed by $$\frac{L}2$$ samples to keep the phase relationship correct.

With a windowed finite-length impulse response (which is what we call an "FIR"), the Hilbert transform output is

$$\hat{x}[n] = \sum\limits_{i=-L/2}^{L/2} h[i] \, x[n-i]$$

and delaying the output so that the filter is causal you get

$$\hat{x}[n-\tfrac{L}2] = \sum\limits_{i=0}^{L} h[i-\tfrac{L}2] \, x[n-i]$$

but with this delayed output $$\hat{x}[n-\tfrac{L}2]$$, you must compare that only to the like delayed input $$x[n-\tfrac{L}2]$$ in order for the two signals to have their 90° phase relationship (which is fundamentally what the Hilbert Transform is about). Note that every even-indexed sample of $$h[n]$$ is zero, so the number of taps is not really $$L+1$$ but is $$\frac{L}2 - 1$$ taps with non-zero coefficients.

• thanks. I guess I must be dumb because I was unable to understand a single line of your answer. Sorry. Suppose for example I already have decomposed the signal using Empirical Mode Decomposition and obtained several IMFs. How can I do a HT on those IMF and get the frequency of each one, for example. Again, sorry for my ignorance. – Duck May 1 '19 at 10:00
• well, Dog, i haven't really figgered out what the IMFs (i guess they're "Intrinsic Mode Functions") are either. there is a semantic usage issue, but the Discrete Hilbert Transform is well defined in the textbooks and the lit. And the definition is what I said at the top. The windowing and delaying is necessary to make your filter FIR and causal. If you're trying to use a Hilbert Transformer to estimate an instantaneous amplitude envelope or instantaneous frequency, that's good for a whole 'nother question. – robert bristow-johnson May 1 '19 at 19:13
• So Dog, if you have a narrow-band discrete-time signal $x[n]$ and its Discrete Hilbert transform is $\hat{x}[n]$, defined as above, the instantaneous frequency in radians/sample (which is dimensionless) is $$\omega[n] = \operatorname{arctan}\left(\frac{x[n-1]\hat{x}[n] - x[n]\hat{x}[n-1]}{x[n]x[n-1]+\hat{x}[n]\hat{x}[n-1]} \right)$$ The instantaneous phase will advance by $\omega[n]$ for each sample $$\phi[n] = \phi[n-1] + \omega[n]$$ – robert bristow-johnson May 1 '19 at 22:25
• Dog, i had to make a small correction. suppose you were to use the Hamming window (which is easier to understand than the Kaiser), can you translate this math into code and get what you're looking for? remember your result is necessarily delayed by $\frac{L}2$ samples and $L$ needs to be reasonably large, like $L$=50 or something. – robert bristow-johnson May 2 '19 at 3:38