The effects of Q are currently incorporated into numerical time-domain simulations of wave propagation by using a set of discrete relaxation functions. The behavior of these relaxation functions is controlled by the selection of the weight coefficients and relaxation times. In this article we present an approach to determine these parameters efficiently and accurately. In our approach we first determine two sets of weight coefficients and one set of relaxation times to model the smallest and largest target Q. We then derive an empirical formula to interpolate the weight coefficients for modeling an arbitrary Q that can be a function of frequency. A single set of relaxation times is used for any Q. We have considered two cases: (1) the complex viscoelastic material modulus is computed by directly summing of relaxation functions; (2) the modulus is implicitly represented as the harmonic average of the viscoelastic modulus given by individual relaxation functions. For large Q (>50), both approaches can be applied to coarse-grained system presented by Day and Bradley (2001). However, for lower Q, the second case works better.