A fast and exact procedure for the numerical synthesis of stationary multivariate Gaussian time series with a priori prescribed and well controlled autoand cross-covariance functions is proposed. It is based on extending the Circulant Embedding technique to the multivariate case and can be viewed as a modification and variation around the Chan and Wood algorithm proposed earlier to solve the same problem. The procedure is shown to yield time series possessing exactly the desired covariance structure, when sufficient conditions are satisfied. Such conditions are discussed theoretically and examined on several examples of multivariate time series models. Issues related to prescribing a priori the spectral structure rather than the covariance one are also discussed. Matlab routines implementing this procedure are publicly available at http://perso.ens-lyon.fr/patrice.abry/software.html.