Spectral Coherence

EEG
#1

Hi all,

I am having a surprisingly difficult time understanding and implementing spectral coherence in MATLAB (actually I need imaginary coherence but one thing at a time). I need to calculate spectral coherence in multiple windows on a single run.

Does this toy code look right to you? Why is coherence all over the place? Shouldn’t this just be a single number? Thanks in advance.

clear all
close all
clc

%%

srate = 1000;
time = 0:1/srate:1;

sine1 = sin(2pi10.*time);
sigx = sine1 + rand(1,1001);
sigy = sine1 + rand(1,1001);

%%

x1 = fft(sigx(1:500))/500;
y1 = fft(sigy(1:500))/500;

xx1 = x1.*conj(x1);
yy1 = y1.*conj(y1);
xy1 = x1.*conj(y1);

%%

x2 = fft(sigx(501:1000))/500;
y2 = fft(sigy(501:1000))/500;

xx2 = x2.*conj(x2);
yy2 = y2.*conj(y2);
xy2 = x2.*conj(y2);

%%

avgXX = (xx1+xx2)/2;
avgYY = (yy1+yy2)/2;
avgXY = (xy1+xy2)/2;

coh = abs(avgXY).^2 ./(avgXX.*avgYY);

figure
hold on
plot(coh)

#2

Hi Thalderson. Spectral coherence is best done with repeated trials that can be averaged over. In this case, the problem is actually with the denominator. Because the signal is so simple, the denominator is basically the same as the numerator, which basically cuts out the peak spectral component and makes the rest of the spectrum highly unstable (dividing 0/0 plus noise).

Observe, for example:
hz = linspace(0,1,length(coh));
plot(hz,avgXX, hz,avgYY, hz,abs(avgXY).^2), set(gca,‘xlim’,[0 .1])

Just plotting the numerator gives the expected coherence profile. Notice also the DC offset (non-zero 0 Hz component, which comes from the uniform noise (average of .5) instead of normal noise (average of 0).

#3

Thank you for your reply Mike - your books and videos which are proving to be an indispensable resource.