First, download and source mvspec.R, which is spec.pgram with a few changes in the defaults:
function(x, spans = NULL, kernel = NULL, taper = 0, pad = 0, fast = TRUE, demean = TRUE, detrend = FALSE, plot = FALSE, na.action = na.fail,...)and written so you can extract the estimate of the multivariate spectral matrix as fxx. For example, if x contains a p-variate time series (i.e., the p columns of x are time series), and you issue the command spec = mvspec(x, spans=3), then spec$fxx is an array with dimensions dim=c(p,p,nfreq), where nfreq is the number of frequencies used. If you print spec$fxx you'll see nfreq p × p spectral matrix estimates.