# How to run Welch’s method on cross-spectral density

I have a multivariate time series and I’m trying to figure out the cross spectral density of all pairs of variables efficiently. I’m unsure of how to implement Welch’s method in this context. Do I calculate the PSD of each series individually using Welch and then use those to compute the cross-spectral density or do I average the cross-spectral density of each sub-band?

If the answer is to average the cross-spectral density of each sub-bad, then how does one write an efficient algorithm to do this for a high dimension?

Thanks!

$$R_{XY}(f) = \lim_{\Delta t \to \infty} {\large\mathbb E} \left\{\frac{X_{[t_0;t_o+\Delta t]}^*(f)Y_{[t_0;t_o+\Delta t]}(f)}{\Delta t}\right\}$$
(abusing the subscript here to denote this is a finite-duration "observation" from point $$t_0$$ to $$t_0+\Delta t$$).