# 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?

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.

• Sorry to butt in here, and so late! But the reply above contains: sxya = (sxy1 + sxy2) / 2; This means that an arithmetic average is performed BEFORE the absolute value of the CPSD is derived - this means you will lose the phase information. Is there not a better way? – Matt Vowels Mar 10 '18 at 4:11