Plotting Hilbert Spectrum from Hilbert Huang Transform
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),&amp;\omega=\omega_j(t)\\0,&amp;\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.