1
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ How would one vary $\omega$ in the spectrum if $\omega$ is already time varying. In my case the signal im putting through is already sampled so $t$ is basically the discrete length of the signal which in my case is 20000 pts long $\endgroup$ – Sam Aug 1 '18 at 14:27
  • $\begingroup$ @Sam: You can evaluate $H_j(\omega,t)$ for any $\omega$ you like, but for a given time $t_k$ it's only non-zero for $\omega_j(t_k)$, i.e., for the current value of the instantaneous frequency. Imagine $H_j(\omega,t)$ as a trace in the time-frequency plane following the function $\omega_j(t)$. For each time instance $t_k$, its value is given by the instantaneous amplitude $a_j(t_k)$. $\endgroup$ – Matt L. Aug 1 '18 at 18:34
  • $\begingroup$ Gotcha, so say I have 3 imfs so there'll be 3 instantaneous frequency vectors with $t = 0, ..., 6000 $, and $ \omega_{j}(t) = ( \omega_{1}(t), \omega_{2}(t), \omega_{3}(t) ) $ say I want the spectrum for just the first imf $j=1$, then we get $H_{1}(\omega, t) = a_{1}(t), for \omega = \omega_{1} (t) $ so then getting numbers say I want the value at $t = 3000$ then it will become $H_{1}(\omega_{1} ( 3000), 3000) = a_{1}(3000) $ , correct? sorry if i was also in unclear in this one. $\endgroup$ – Sam Aug 1 '18 at 20:19
  • $\begingroup$ @Sam: Yes, that's it. $\endgroup$ – Matt L. Aug 2 '18 at 14:53
  • $\begingroup$ @Sam: If the answer was helpful you can accept it by clicking on the green check mark, thanks! $\endgroup$ – Matt L. Aug 3 '18 at 8:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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