1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ If the complex addition is destroying the phase information, why are you trying to combine them? $\endgroup$ – user28715 Mar 10 '18 at 10:14
  • $\begingroup$ The order of the averaging makes a difference - imagine two CPSDs which are very similar but are approx. polarity inverted. When I average them in the way proposed by A. Donda, much of the CPSD will cancel. However, if I first take the magnitude squared BEFORE averaging, this won't happen. As such, I am unable to find a way to correctly average MSCs from multiple pairs of time series. $\endgroup$ – Matt Vowels Mar 10 '18 at 15:50
  • $\begingroup$ If there is a polarity difference, you need to rotate one signal to coherently align to the other. If the CSPD phases are not coherent, you should explain why you want to average them. A Donda is correct for the question asked. If there is some random phase perturbation present, they really aren’t coherent, so why is that not correct? $\endgroup$ – user28715 Mar 10 '18 at 16:22
  • $\begingroup$ An MSC of a single pair of time series is agnostic to the phase - it will tell you whether the spectral content is similar, even if they are in reverse polarity. If we wish to know the phase as well, then the angle of the CPSD can be used, but this is separate. Hence, I wish to estimate the average spectral coherence, regardless of phase. Using A Donda's method means that any CPSDs that are 'out of phase' will cancel. So, I really want to take the averages of the magnitudes of the CPSDs before the average MSC is computed. I hope that makes some sense? $\endgroup$ – Matt Vowels Mar 10 '18 at 16:33
  • $\begingroup$ You seem to want an incoherent coherence. Good Luck with that. Cliff Carter wrote some very good papers on coherence and time delay estimation. $\endgroup$ – user28715 Mar 10 '18 at 16:53
-1
$\begingroup$

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.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.