function [plvmat, plvts] = PLV(x, fs, freqrange)
%% Phase-Locking Value
% PLV is a measure of functional connectivity, suitable for EEG data. 
% Denote two time series as x and y (for example, two channels.) Step by step: 
%   1) Bandpass the data to a narrow filter
%   2) Use the hilbert transform to get instantaneous phase of the data (where
%      you are at the oscillation around a unit-circle).
%   3) For every time point in every trial, calculate the difference
%      between x and y. Now we need to quantify this difference. Euler's
%      formula states that exp(ix) = cosx + i sinx where i is the imaginary
%      number. Replace x with the phase difference (x := phase(x) -
%      phase(y)). This moves the phase difference to a complex plane where
%      x axis denotes real and y axis denotes complex parts. Let's say that
%      v = a + ib and v is the exp(i(phase(x)-phase(y)). The magnitude of 
%      this vector r = sqrt(a^2 + b^2) would give how big the phase
%      difference is. 
%   4) Average every time point across trials. This gives you whether the
%      phase differences are consistent i.e. two channels are phase-locked
%      or not. 
% 
% Input:
%   x: time series (3D: channels x time x trials)
%   fs: Sampling rate (scalar)
%   freqrange: frequency of interest (for example, [4 7] for theta band
% 
% Output:
%   plvmat: PLV averaged across time (n x n matrix where ns are channels.)
%   plvts: PLV at every time point (n x n x time where ns are channels).
% This function needs fieldtrip toolbox to function.
% Authored by Yasir Çatal
% catalyasir@gmail.com
%% The code
nchan = size(x, 1);
ntime = size(x, 2);
ntrials = size(x, 3);

x = reshape(x, nchan, []); % reshape data to make it 2D
filtered = ft_preproc_bandpassfilter(x, fs, freqrange, [], 'fir');
filtered = reshape(filtered, nchan, [], ntrials); % Reshape again to get trials

plvts = zeros(nchan, nchan, ntime);
plvmat = zeros(nchan);

for a = 1:nchan
    x = angle(hilbert(squeeze(filtered(a, :, :))));
    for b = (a+1):nchan % Note the i + 1. Don't do calculations twice. Only the upper triangle will be relevant.
        y = angle(hilbert(squeeze(filtered(b, :, :))));
        plvts(a, b, :) = abs(mean(exp(1i*(x - y)')));
        plvmat(a, b) = mean(plvts(a, b, :));
    end
end
end