We propose a spectral turning-bands approach for the simulation of second-order stationary vector Gaussian random fields. The approach improves existing spectral methods through coupling with importance sampling techniques. A notable insight is that one can simulate any vector random field whose direct and cross-covariance functions are continuous and absolutely integrable, provided that one knows the analytical expression of their spectral densities, without the need for these spectral densities to have a bounded support. The simulation algorithm is computationally faster than circulant-embedding techniques, lends itself to parallel computing and has a low memory storage requirement. Numerical examples with varied spatial correlation structures are presented to demonstrate the accuracy and versatility of the proposal.