power spectral density - How to average coherences estimated using Welch's method in MATLAB

I'm trying to estimate the average coherence of discontinuous snippets of a signal in MATLAB. For Welch's spectrum estimates, the following works:

w = blackmanharris(50);  % window
n = numel(w);            % fft length
o = 0;                   % overlap

x = rand(200,1); x1 = x(1:100); x2 = x(101:200);

s  = pwelch(x,  w, o, n);
s1 = pwelch(x1, w, o, n);
s2 = pwelch(x2, w, o, n);
sa = (s1+s2)/2;


max(abs(s-sa))  % almost zero

However, when estimating coherences using Welch's method, the arithmetic mean yields an incorrect result:

y = rand(200,1); y1 = y(1:100); y2 = y(101:200);

c  = mscohere(x,  y,  w, o, n);
c1 = mscohere(x1, y1, w, o, n);
c2 = mscohere(x2, y2, w, o, n);
ca = (c1+c2)/2;


max(abs(c-ca))  % very different from zero

How can I average two coherence estimates correctly?

Answer

If you look at the documentation of mscohere which contains the definition of magnitude squared coherence, you will see that it is the absolute squared of the cross spectral density, divided by the product of the two separate spectral densities. Since power densities enter numerator and denominator of the expression, the MSC estimate from a full time series is not simply the average of the estimates from the parts.

To get a good "average" MSC estimate, do the following: estimate in each segment separately cross and single spectral densities, average these three over the segments, and then combine them in the end into an average MSC. Continuing your code:

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;
ca = abs(sxya) .^ 2 ./ sxa ./ sya;

Now max(abs(c-ca)) is almost zero.