# Phase locking value / phase synchronization

I have pairs of EEG signals filtered to give me just alpha bands. For a machine learning problem I have found that calculating the phase locking value between the alpha signals of these 2 channels would be a useful classification feature.

I'm using this paper as a reference: http://www.ncbi.nlm.nih.gov/pubmed/22661936

I'm wondering what the correct method/implementation is to calculate PLV/synchrony in Python?

I have not been able to find any modules/code that would help me with this problem. That paper, for example, calculates synchrony between independent components after ICA, but I'm not employing that method here so not sure if their formula translates or if any adjustments are necessary for filtered signals

The closest thing I have found is this: https://www.nbtwiki.net/doku.php?id=tutorial:phase_locking_value#.VchjR6ZVhBc

However, while the code is laid out nearer the bottom of the page, it appears to be MATLAB code. It looks like there are some MATLAB methods being employed and I'm not sure how that would translate over to Python/numpy/scipy

I'm not particularly familiar with this domain and don't trust myself to port the code over without making logic errors, so it would be useful if someone could point me toward a Python implementation of this.

There can be several ways to calculate the Phase locking value (PLV). For relatively mono-component and high SNR (well filtered)-Time domain signal can be converted into analytical signal using Hilbert transform to calculate the phase difference. For the right signal it is a very powerful technique as is shown in the tutorial you have referenced. Here is a code in python : y1 and y2 are the input signals. The function gives average phase difference between the two signals even when they are not relatively well behaved. If you run it in a loop, it can give you phase difference in parts or even sample by sample (instantaneous phase).

import numpy as np
import scipy.signal as sig

def hilphase(y1,y2):
sig1_hill=sig.hilbert(y1)
sig2_hill=sig.hilbert(y2)
pdt=(np.inner(sig1_hill,np.conj(sig2_hill))/(np.sqrt(np.inner(sig1_hill,
np.conj(sig1_hill))*np.inner(sig2_hill,np.conj(sig2_hill)))))
phase = np.angle(pdt)

return phase


For well behaved signal as shown in the tutorial even a more simple instantaneous phase program can be easily written :

import numpy as np
import scipy.signal as sig

def hilphase(y1,y2):
sig1_hill=sig.hilbert(y1)
sig2_hill=sig.hilbert(y2)
phase_y1=np.unwrap(np.angle(sig1_hill))
phase_y2=np.unwrap(np.angle(sig2_hill))
Inst_phase_diff=phase_y1-phase_y2
avg_phase=np.average(Inst_phase_diff)
return Inst_phase_diff,avg_phase

• What is the phase return variable in your first function? – Simon Jun 27 '17 at 11:45
• Sorry abt that edited now. – Perscitius Jun 27 '17 at 11:54

If you have two numpy arrays of phase data theta1 and theta2 (in radians), you can calculate phase locking value in numpy without too much effort:

import numpy as np

def phase_locking_value(theta1, theta2):
complex_phase_diff = np.exp(np.complex(0,1)*(theta1 - theta2))
plv = np.abs(np.sum(complex_phase_diff))/len(theta1)
return plv


I would also strongly recommend that you use the circular equivalent of pearson's correlation coefficient instead of phase locking value. From a mathematical perspective circular correlation coefficient is preferable since it has a well understood statistical interpretation. From a feature engineering perspective, circular correlation is preferable because it ranges from -1 to 1 instead of 0 to 1. (Ranging from 0 to 1 can misrepresent variations for weakly phase locked signals).

Have you had a look at mne-python? It has routines for inter-trail PL-factor:

def _induced_power_cwt(data, sfreq, frequencies, use_fft=True, n_cycles=7,
decim=1, n_jobs=1, zero_mean=False):
"""Compute time induced power and inter-trial phase-locking factor


Separately, this spectral connectivity routine provides for PLV:

    'plv' : Phase-Locking Value (PLV) [2]_ given by: PLV = |E[Sxy/|Sxy|]|

class _PLVEst(_EpochMeanConEstBase):
"""PLV Estimator"""


Here is an example from their examples page, computing power and phase lock in source space:

• I have tried to use MNE as a module but its really messy because it requires data to flow in a specific format which is rather limiting. This is why I was looking for a fresh implementation of the formula given in the paper/website I cited. I have also taken a look at the code used by MNE but again, difficult to follow the process because they requires preprocessing (I think CSD is required first) before calculating of PLV – Simon Aug 28 '15 at 5:29
• You did not mention that you'd already investigated & disqualified mne-python. I wasn't able to find any other PLV python implementations (on github, google code or elsewhere). However, I just inspected the MATLAB code that you mentioned, and I think that the code is simple enough that it could be directly translated to python numpy, with minimal fuss. I think you should try it out, you'd need numpy & scipy.signal.hilbert and if you get stuck there's always StackOverflow. – ruoho ruotsi Aug 28 '15 at 6:01