A model is proposed to simulate multivariate weakly stationary Gaussian stochastic processes based on the spectral representation theorem. In this model, the amplitude, phase angle, and frequency involved in the harmonic function are random so that the generated samples are real stochastic processes. Three algorithms are then adopted to improve the simulation efficiency. A uniform cubic B-spline interpolation method is employed to fit the target factorized power spectral density function curves. A recursive algorithm for the Cholesky factorization is utilized to decompose the cross-power spectral density matrices. Some redundant cosine terms are cut off to decrease the computation quantity of superposition. Finally, an example involving simulation of turbulent wind velocity fluctuations is given to validate the capability and accuracy of the proposed model as well as the efficiency of the optimal algorithms.