We present a method to calculate classical lattice-dynamical quantities, such as the temperature-dependent vibrational spectrum, from simulations that do not require an explicit solution of the time evolution. We start from the moment expansion of the relevant time-correlation function for a many-body system and show that it can be conveniently rewritten by using a basis in which the low-order moments are diagonal. This allows us to approximate the main spectral features (i.e., position and width of the phonon peaks) from thermal averages available from any statistical simulation. We illustrate our method with an application to a model system that presents a structural transition and strongly temperature-dependent phonons. Our theory also clarifies the status of previous heuristic schemes to estimate phonon frequencies.