The present paper introduces a new approach to simulate any stationary multivariate Gaussian random field whose cross-covariances are predefined continuous and integrable functions. Such a field is given by convolution of a vector of univariate random fields and a functional matrix which is derived by Cholesky decomposition of the Fourier transform of the predefined cross-covariance matrix. In contrast to common methods, no restrictive model for the cross-covariance is needed. It is stationary and can also be reduced to the isotropic case. The computational effort is very low since fast Fourier transform can be used for simulation. As will be shown the algorithm is computationally faster than a recently published spectral turning bands model. The applicability is demonstrated using a common numerical example with varied spatial correlation structure. The model was developed to support simulation algorithms for mineral microstructures in geoscience.