We developed a general method to generate populations of artificial spike trains (ASTs) that match the statistics of recorded neurons. The method is based on computing a Gaussian local rate function of the recorded spike trains, which results in rate templates from which ASTs are drawn as gamma distributed processes with a refractory period. Multiple instances of spike trains can be sampled from the same rate templates. Importantly, we can manipulate rate-covariances between spike trains by performing simple algorithmic transformations on the rate templates, such as filtering or amplifying specific frequency bands, and adding behavior related rate modulations. The method was examined for accuracy and limitations using surrogate data such as sine wave rate templates, and was then verified for recorded spike trains from cerebellum and cerebral cortex. We found that ASTs generated with this method can closely follow the firing rate and local as well as global spike time variance and power spectrum. The method is primarily intended to generate well-controlled spike train populations as inputs for dynamic clamp studies or biophysically realistic multicompartmental models. Such inputs are essential to study detailed properties of synaptic integration with well-controlled input patterns that mimic the in vivo situation while allowing manipulation of input rate covariances at different time scales.