  We present a fast and simple algorithm that allows the extraction of multiple
exponential signals from the temporal dependence of correlation functions
evaluated on the lattice including the statistical fluctuations of each signal
and treating properly backward signals. The basic steps of the method are the
inversion of appropriate matrices and the determination of the roots of an
appropriate polynomial, constructed using discretized derivatives of the
correlation function. The method is tested strictly using fake data generated
assuming a fixed number of exponential signals included in the correlation
function with a controlled numerical precision and within given statistical
fluctuations. All the exponential signals together with their statistical
uncertainties are determined exactly by the algorithm. The only limiting factor
is the numerical rounding off. In the case of correlation functions evaluated
by large-scale QCD simulations on the lattice various sources of noise, other
than the numerical rounding, can affect the correlation function and they
represent the crucial factor limiting the number of exponential signals,
related to the hadronic spectral decomposition of the correlation function,
that can be properly extracted. The algorithm can be applied to a large variety
of correlation functions typically encountered in QCD or QCD+QED simulations on
the lattice, including the case of exponential signals corresponding to poles
with arbitrary multiplicity and/or the case of oscillating signals. The method
is able to to detect the specific structure of the multiple exponential signals
without any a priori assumption and it determines accurately the ground-state
signal without the need that the lattice temporal extension is large enough to
allow the ground-state signal to be isolated.
