A diffuse sound field is conventionally defined as a zero-mean circularly symmetric complex Gaussian random field. A more recent, generalized definition is that of a sound field having mode shapes that are diffuse in the conventional sense, and eigenfrequencies that conform to the Gaussian orthogonal ensemble. Such a generalized diffuse sound field can represent a random ensemble of sound fields that share gross features, such as modal density and total absorption, but otherwise have any possible arrangement of local wave scattering features. The problem of generating realizations or Monte Carlo samples of a conventional diffuse sound field or, equivalently, of the mode shapes of a generalized diffuse sound field, is addressed here. Such realizations can be obtained from an eigenvalue decomposition of the spatial correlation function. A discrete decomposition is numerically expensive when the sound pressures at many locations are of interest, so a fast analytical decomposition based on prolate spheroidal wave functions is developed. The approach is numerically validated by comparison with a detailed room model, where random wave scatterers are explicitly modeled as acoustic point masses with random positions, and good correspondence is observed. Furthermore, applications involving correlated sound sources and sound-structure interaction are presented.