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