function [plimat, plits] = PLI(x, fs, freqrange)
%% Phase-Lag Index
% PLI 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. If x is bigger (difference > 0) encode it as +1.
%      Otherwise, -1.
%   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 (consistently lagging or consistently dragging; not quite the 
%      same tempo). 
% Of note, PLI was found to be more resilient against volume conduction
% compared to PLV (see this website for nice figures and references: 
% https://sapienlabs.org/lab-talk/eeg-connectivity-using-phase-lag-index/).
% 
% 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:
%   plimat: PLI averaged across time (n x n matrix where ns are channels.)
%   plits: pli 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

plits = zeros(nchan, nchan, ntime);
plimat = zeros(nchan);

for i = 1:nchan
    x = angle(hilbert(squeeze(filtered(i, :, :))));
    for j = (i+1):nchan % Note the i + 1. Don't do calculations twice. Only the upper triangle will be relevant.
        y = angle(hilbert(squeeze(filtered(j, :, :))));
        plits(i, j, :) = abs(mean(sign(x - y)'));
        plimat(i, j) = mean(plits(i, j, :));
    end
end
end
