How to run Welch’s method on cross-spectral density
Let's define Cross-Spectral Density; you can do that via the usual &quot;Fourier transform of the cross-correlation&quot;, 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 &quot;observation&quot; 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.