2
$\begingroup$

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:

enter image description here

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?

$\endgroup$

1 Answer 1

3
$\begingroup$

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{1 - (-1)^n}{\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 non-zero 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 FIR taps (the FIR filter order is $L$) and, in my centered and symmetrical case, $L$ must be even. $\beta$ is a "shape parameter", maybe around 5 or 6.

Because of the half-band symmetry, every other tap has a zero coefficient. So $\frac{L}2+1$ is the number of non-zero taps and $\frac{L}2$ is an odd integer. That means $L=2+4m$ where $m \in \mathbb{Z} \ge 0$.

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 having non-zero coefficients is not really $L+1$ but is $\frac{L}2+1$. There are $\frac{L}2$ equally-spaced taps having zero-valued coefficients that can be skipped over in the FIR.

$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Duck
    Commented May 1, 2019 at 10:00
  • $\begingroup$ 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. $\endgroup$ Commented May 1, 2019 at 19:13
  • 1
    $\begingroup$ 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] $$ $\endgroup$ Commented May 1, 2019 at 22:25
  • $\begingroup$ 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. $\endgroup$ Commented May 2, 2019 at 3:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.