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.