In this paper, we propose algorithms which preserve energy in empirical mode decomposition (EMD), generating finite $n$ number of band limited Intrinsic Mode Functions (IMFs). In the first energy preserving EMD (EPEMD) algorithm, a signal is decomposed into linearly independent (LI), non orthogonal yet energy preserving (LINOEP) IMFs and residue (EPIMFs). It is shown that a vector in an inner product space can be represented as a sum of LI and non orthogonal vectors in such a way that Parseval's type property is satisfied. From the set of $n$ IMFs, through Gram-Schmidt orthogonalization method (GSOM), $n!$ set of orthogonal functions can be obtained. In the second algorithm, we show that if the orthogonalization process proceeds from lowest frequency IMF to highest frequency IMF, then the GSOM yields functions which preserve the properties of IMFs and the energy of a signal. With the Hilbert transform, these IMFs yield instantaneous frequencies and amplitudes as functions of time that reveal the imbedded structures of a signal. The instantaneous frequencies and square of amplitudes as functions of time produce a time-frequency-energy distribution, referred as the Hilbert spectrum, of a signal. Simulations have been carried out for the analysis of various time series and real life signals to show comparison among IMFs produced by EMD, EPEMD, ensemble EMD and multivariate EMD algorithms. Simulation results demonstrate the power of this proposed method.