Wave propagation analysis of heterogeneous unbounded media employing domain discretization methods necessitates truncation of the finite computational domain at its truncation boundaries. The Perfectly Matched Layer (PML) is a robust absorbing boundary condition, which is extensively employed to favorably absorb outgoing scattered waves from the interior computational domain, regardless of frequencies and directions of outgoing waves. Nonetheless, most available PMLs are developed for elastic domains, while their developments in domain discretization methods, which consider damping using the Rayleigh damping formulation, are insufficient. In this study, a viscoelastic PML formulation is developed that accommodates the Rayleigh-type damping. The PMLs are basically developed from a coordinate transformation concept using a complex stretching function. Suitable stretching functions can efficiently control the absorptive properties and the spectral character of the PML. Besides, the stretching functions play a crucial role in the long-time stability of transient wave simulations. The standard stretching function usually results in spurious reflections, especially when propagating waves incident to the truncation boundary at near-grazing incidence. On the other hand, the complex-frequency-shifted (CFS) stretching function has been shown to mitigate long-time instability. Therefore, the proposed viscoelastic PML in the time-domain and frequency-domain are derived using the CFS stretching function for the first time. Comprehensive formulations of the proposed viscoelastic CFS-PML are given, where the displacement field is represented as the sole unknown variable of the problem. Hence, the proposed CFS-PML can be incorporated into any domain discretization method. However, the proposed CFS-PML is further developed in the context of the voxel element method, where element-by-element iterative strategies are implemented to solve the governing equations without storing global matrices in the computer memory. Comprehensive numerical experiments are conducted to demonstrate the validity and efficiency of the proposed viscoelastic CFS-PML formulation in analyzing soil profiles with various damping ratios and properties. Moreover, long-time stability for problems involving waveguides and near-grazing wave incidence is studied, where no instabilities are detected.