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

Question

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!

Answer 1

Let's define Cross-Spectral Density; you can do that via the usual "Fourier transform of the cross-correlation", but assuming WSS, Wiener-Khinchin allows us to state it as

$$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$).

Now, you could of course first estimate the cross-correlation and then estimate the spectrum of that, and end up with an estimate of the CSD; however, estimation of a cross-correlation is usually done by first transforming the both signals to frequency domain, (conjugate) multiplying there, and then transforming back.

Transforming that back into frequency domain seems like you've basically reverted the last step.

Small problem with using Welch here: it usually ends with squaring the transform of your averaged segments – that erases the phase information necessary in the formula above.