This article proposes a new method for simulating a Gaussian process, whose spectrum diverges at one or multiple frequencies in [0, 1/2] (not necessarily at zero). The method uses a generalization of the discrete wavelet transform, the discrete wavelet packet transform (DWPT), and requires only explicit knowledge of the spectral density function of the process—not its autocovariance sequence. An orthonormal basis is selected such that the spectrum of the wavelet coefficients is as at as possible across specific frequency intervals, thus producing approximately uncorrelated wavelet coefficients. This method is compared to a popular time-domain technique based on the Levinson–Durbin recursions. Simulations show that the DWPT-based method performs comparably to the time-domain technique for a variety of sample sizes and processes—at significantly reduced computational time. The degree of approximation and reduction in computer time may be adjusted through selection of the orthonormal basis and wavelet filter.