How to calculate phase synchrony?

hello!

I am trying to see the relationship between two brainwave signals measured at two different locations on the head. The two signals are theta medians going up and down usually between float values 1.0 and 5.0. They are updated 25 times per second.

The animated gif below demonstrates the spectral data of 8 EEG signals (wiggly lines) and 8 theta medians (short horizontal lines) :

eeg-theta-median

So far, I only managed to calculate the amplitude difference between the theta signals with:

(thetaMed0-thetaMed1).abs

However, this measurement does not seem to be useful for seeing phase synchrony as two signals can be synchronous (in-phase) even if their crests and troughs are not the same:

40

Could you help with some guidance on calculating phase synchrony between two signals?

Cheers! k

http://paulbourke.net/miscellaneous/correlate/

The particular bit of magic in the second link is: “The Fourier transform of the cross correlation function is the product of the Fourier transform of the first series and the complex conjugate of the Fourier transform of the second series.”

So let’s manufacture a signal where y is shifted to the right by 25 samples (conversely, x is shifted by -25 samples).

x = Signal.sineFill(512, [1]);
y = x.rotate(25);

The magic formula requires FFTs.

~table = Signal.fftCosTable(x.size);
~xfft = x.fft(Signal.newClear(x.size), ~table);
~yfft = y.fft(Signal.newClear(y.size), ~table);

Then, if we multiply the FFT of the test signal (y) against the complex conjugate of the reference signal (x)'s FFT, we get the FFT of the crosscorrelation.

~crosscorrfft = ~yfft * ~xfft.conjugate;

// so the cross-correlation function is the ifft of that

~crosscorr = ~crosscorrfft.real.ifft(~crosscorrfft.imag, ~table);
~crosscorr.real.plot;
~crosscorr.real.maxIndex;  // 25

And the ‘x’ coordinate of the maximum value of that function is the amount of shift for the test signal. If you shift the y signal to the left by 25 samples, you should get the maximum amount of correlation against the reference.

When x and y are measured from different sources, there may be a number of local maxima, all reflecting amounts of shift where correlation is relatively high.

(I happened to have looked into autocorrelation a couple of years ago :stuck_out_tongue_winking_eye: which is a special case of this where x == y.)

hjh

2 Likes

Or, if you’re only interested in the degree of correlation right now without shifting either signal, sum(x * y) will be high if x and y are highly correlated, and lower if they are not correlated. (The cross-correlation function checks every possible amount of shift.)

hjh

1 Like

Many thanks James,

this is pretty cool!

To understand what your code does I compared 4 examples (maybe useful for others):

// example 1:
// same signal
(
x = Signal.sineFill(512, [1]); //red
y = Signal.sineFill(512, [1]); //green
z=(x*y); //blue
a=sum(x*y).postln; //260
)

//degree of correlation: 260

( //red cannot be seen as green is on top;
var data = [x,y,z];
var plotter = Plotter("plotter",Rect(5,5,630,235));
plotter.plotMode = \points;
plotter.value_(data);
plotter.superpose_(true);
plotter.setProperties(
	\plotColor, [Color.red,Color.green,Color.blue]
);
);

(
~table = Signal.fftCosTable(x.size);
~xfft = x.fft(Signal.newClear(x.size), ~table);
~yfft = y.fft(Signal.newClear(y.size), ~table);
~crosscorrfft = ~yfft * ~xfft.conjugate;
~crosscorr = ~crosscorrfft.real.ifft(~crosscorrfft.imag, ~table);
~crosscorr.real.plot;
~crosscorr.real.maxIndex;  // 0
)

// 0 means no shift of phase between the two signals



//example 2:
//slightly out of phase signals with same amp
(
x = Signal.sineFill(512, [1]); //red
y = x.rotate(25); //green
z=(x*y); //blue
a=sum(x*y).postln; //244
)

//degree of correlation: 244

(
var data = [x,y,z];
var plotter = Plotter("plotter",Rect(5,5,630,235));
plotter.plotMode = \points;
plotter.value_(data);
plotter.superpose_(true);
plotter.setProperties(
	\plotColor, [Color.red,Color.green,Color.blue]
);
);

(
~table = Signal.fftCosTable(x.size);
~xfft = x.fft(Signal.newClear(x.size), ~table);
~yfft = y.fft(Signal.newClear(y.size), ~table);

~crosscorrfft = ~yfft * ~xfft.conjugate;

~crosscorr = ~crosscorrfft.real.ifft(~crosscorrfft.imag, ~table);
~crosscorr.real.plot;
~crosscorr.real.maxIndex;  // 25
)

// 25: means signals are 25 'frames' out of phase;



// example 3:
// same phase, different amp;
(
x = Signal.sineFill(512, [1]); //red
y = Signal.sineFill(512, [1])*0.5; //green
z=(x*y); //blue
a=sum(x*y).postln;//128
)

//degree of correlation: 128

(
var data = [x,y,z];
var plotter = Plotter("plotter",Rect(5,5,630,235));
plotter.plotMode = \points;
plotter.value_(data);
plotter.superpose_(true);
plotter.setProperties(
	\plotColor, [Color.red,Color.green,Color.blue]
);
);

(
~table = Signal.fftCosTable(x.size);
~xfft = x.fft(Signal.newClear(x.size), ~table);
~yfft = y.fft(Signal.newClear(y.size), ~table);

~crosscorrfft = ~yfft * ~xfft.conjugate;

~crosscorr = ~crosscorrfft.real.ifft(~crosscorrfft.imag, ~table);
~crosscorr.real.plot;
~crosscorr.real.maxIndex;  //0
)

// 0: means signals are 0 'frames' out of phase;



// example 4:
// more out of phase and different amp
(
x = Signal.sineFill(512, [1]); //red
y = x.rotate(100)*0.5; //green
z=(x*y); //blue
a=sum(x*y).postln; //43
)

//degree of correlation: 43

(
var data = [x,y,z];
var plotter = Plotter("plotter",Rect(5,5,630,235));
plotter.plotMode = \points;
plotter.value_(data);
plotter.superpose_(true);
plotter.setProperties(
	\plotColor, [Color.red,Color.green,Color.blue]
);
);

(
~table = Signal.fftCosTable(x.size);
~xfft = x.fft(Signal.newClear(x.size), ~table);
~yfft = y.fft(Signal.newClear(y.size), ~table);

~crosscorrfft = ~yfft * ~xfft.conjugate;

~crosscorr = ~crosscorrfft.real.ifft(~crosscorrfft.imag, ~table);
~crosscorr.real.plot;
~crosscorr.real.maxIndex;  // 100
)

// 100: means signals are 100 'frames' out of phase;

Maximum degree of correlation is around sample-rate/2 i.e. in the examples above 512/2.

Maybe turning some of the values absolute before plotting could help see correlations, I will experiment. I will try to implement in the EEG code now…

Cheers! Have a great day, k

I posted here a few minutes ago, removing the reply, sorry, I might have found a clue…

If x == y, then sum(x * y) is the same as sum(x ** 2) – which depends on the magnitudes. There was never any guarantee of a single fixed value representing perfect correlation.

But… if you know that sum(x ** 2) is perfect correlation, then sum(x * y) / sum(x.squared) is a fraction where 1.0 always means perfect correlation (though, if y is on average greater than x, then this will be > 0 so still not necessarily perfect).

Usually auto-correlation and cross-correlation are looking for the amount of shift that will line up the signals as closely as possible (the FFT approach).

The validity of FFT doesn’t depend on the sampling rate at all. 25 samples/second is as valid as 44,100 – the only difference is the maximum frequency that can be represented (half sampling rate).

hjh

1 Like

Actually – I just realized I misrepresented the simpler approach – I should’ve known that sum(x * y) isn’t an absolute measurement. sum(x * y) would be valid if you have one signal ‘x’ and many different ‘y’ signals, and you’re trying to figure out which ‘y’ is most closely correlated against ‘x’. If both x and y are constantly changing, then there’s no stable point of reference and that value doesn’t mean anything.

I guess you would have to normalize both x and y’s center and range. Not well tested:

f = { |array|
	array = array - array.mean;  // remove DC offset
	array = array / (array.abs.maxItem);
};

DC offset would have a cumulative effect in the sum, so, remove it. Then, scale everything to -1.0 … 1.0 (without introducing new DC offset). But, this would be sensitive to spikes in the input.

And…

Your question about FFT made me think of something else – is the EEG equipment filtering the readings before you take the 25 samples per second? If not, then you might be getting aliasing in your readings.

// here we have a sine wave
// that is below the Nyquist frequency
// of a high sampling rate
~high_rate = Signal.fill(88200, { |i| sin((2pi * 44000 / 88200) * i) });

// the plot shows very fast oscillation
~high_rate[0..500].plot;

// and we naively decimate by a factor of 2
~low_rate = ~high_rate[0, 2..];

// but now the plot is very slow!
// because 44000 Hz has folded over
// around the new Nyquist frequency
~low_rate[0..500].plot;

Before downsampling anything, it’s necessary to filter out frequencies above the new Nyquist. The analog to digital converter in your soundcard does this for you, by design. If the EEG is giving you 25 samples per second, then it’s probably filtering (but I don’t know the device specs). If the EEG is giving you a higher data rate and you’re polling 25 times per second, then your results may not be valid.

hjh

1 Like