Averaging magnitude squared coherence across multiple time series In a previous post, A. Donda had suggested that, in order to calculate the average magnitude squared coherence of more than one pair of time series (e.g. y1 and x1, and y2 and x2), one ought to follow the following steps:

sx1 = pwelch(x1, w, o, n);
sy1 = pwelch(y1, w, o, n);

sxy1 = cpsd(x1, y1, w, o, n);

sx2 = pwelch(x2, w, o, n);
sy2 = pwelch(y2, w, o, n);

sxy2 = cpsd(x2, y2, w, o, n); 

sxa = (sx1 + sx2) / 2;
sya = (sy1 + sy2) / 2;

sxya = (sxy1 + sxy2) / 2; % NOTE THIS LINE %
ca = abs(sxya) .^ 2 ./ sxa ./ sya;


...where w,o,n are the window parameters, and the other variables should (hopefully) be self-explanatory.

Let me illustrate why this method can be problematic with an example, in the hope that someone might be able to help me with a solution. 

Imagine two pairs of signals y1, x1, and y2,x2. All four signals contain the same frequency component (e.g. a 440Hz sine wave). In one pair, the sines waves are in phase w.r.t. each other. In the other pair, the sine waves are in inverse polarity.

Separate MSCs of both of these pairs of signals will yield approx. '1' at 440Hz in both results, identically, telling us that 440Hz is common to both signals in both pairs (i.e. the relative phase of each of the signals within a pair is ignored). If we now wish to calculate the average MSC, we must first compute the average PSDs and then the average CPSDs. Ideally, this average MSC for the two pairs would also yield a '1' at 440Hz. BUT the CPSD output is complex, so averaging the two CPSDs of the two pairs of signals according to the method above will result in cancellation, because the averaging is taken BEFORE the magnitude of the CPSDs is calculated. 

I have offered one solution below - is this correct?
I believe I might have a solution... Editing A Donda's work above, I think a simple rearrangement of the last two lines is necessary:

sxya = (abs(sxy1) + abs(sxy2)) / 2; 
ca = sxya.^ 2 ./ sxa ./ sya;

So here, the average of the magnitudes is taken, rather than the average of the complex CPSDs. Then, this average can be squared.