spexmat = function(specobj) {
#
# add the cross spectral density matrix to a spectrum object
#
# returns the original object with $xmat member, which is an array
#   of dimension ns * ns * nf
# where ns = number of series and nf = number of frequencies
#
	if (is.null(specobj[["coh"]]))
		return(specobj);
	spec = specobj[["spec"]];
	coh = specobj[["coh"]];
	phase = specobj[["phase"]];
	p = ncol(spec);
	nf = nrow(spec);
	xmat = array(NA, c(p, p, nf));
	dimnames(xmat) = list(specobj[["snames"]], specobj[["snames"]],
				character(0));
	for (k in 1:nf) {
		for (i in 1:p) {
			xmat[i, i, k] = spec[k, i];
			for (j in (i:p)[-1]) {
				l = i + (j - 1) * (j - 2)/2;
				xmat[i, j, k] =
					sqrt(spec[k, i] * spec[k, j] *
						coh[k, l]) *
					complex(real = cos(phase[k, l]),
						imag = sin(phase[k, l]));
				xmat[j, i, k] = Conj(xmat[i, j, k]);
			}
		}
	}
	specobj[["xmat"]] = xmat;
	specobj;
}