Plotting Hilbert Spectrum from Hilbert Huang Transform

Question

Before anyone calls it a repeat post from past posts. I'd like to point out I've read all of those; however, I have questions on making a hard coded version of the hilbert spectrum.

I do indeed have all the ingredients aka the hilbert transform and empirical mode decompositions.

Here is my code below:

def hilbert_huang(self, imf, res):
    # Reconstruct the signal
    sig = np.sum(imf, axis=0) + res
    # Hilbert Transform of IMFs
    hil_imf = np.array([sp.fftpack.hilbert(imf[i]) for i in range(len(imf))])
    # Instantaneous Amplitude
    inst_amp = np.array([np.sqrt(imf[i]**2 + hil_imf[i]**2) for i in range(len(imf))])
    # Instaneous Phase
    inst_phase = np.array([np.arctan(hil_imf[i]/imf[i]) for i in range(len(imf))])
    # Instaneous Temporal Frequency
    inst_freq = np.array([(1/(2*np.pi))* np.gradient(inst_phase[i]) for i in range(len(imf))])

    return inst_amp, inst_phase, inst_freq, sig

However, I'm not sure how to construct the spectrum from all this, I know that the formula is given by $$ H(\omega, t) = \Re\left(\sum_{i=1}^{n} a_{i}(t) e^{j \int \omega_{i} (t) dt}\right) $$

But I'm not sure how to fix $\omega$ , given that $\omega$ is an array with values , in the wikipedia article it says to let $\omega = \omega_{j}$ where $j$ indicates which IMF it's from but again $\omega$ is a vector, so do I just find the max? Let me know if I'm not clear or need to provide any more mathematical expressions or code.

Answer 1

The formula for the Hilbert spectrum in your question does not make much sense. Its right-hand side it is not a function of frequency, instead it is simply the sum of the intrinsic mode functions (IMFs), i.e., it represents the original signal. The Hilbert spectrum is the sum of the individual Hilbert spectra of each IMF:

$$H_j(\omega,t)=\begin{cases}a_j(t),&\omega=\omega_j(t)\\0,&\text{otherwise}\end{cases}$$

where $a_j(t)$ and $\omega_j(t)$ are the instantaneous amplitude and the instantaneous frequency (IF) of the $j^{\,th}$IMF, respectively.

Apart from that, note that your computation of the instantaneous phase and frequency will not work. For the phase you have to use arctan2 instead of arctan, and for computing the IF from the instantaneous phase, you need to unwrap the phase. However, it is generally better to compute the IF directly from the signal, as discussed in this answer.