计算脑电信号相位锁定值 Phase Locking Value

PLV：相位同步信息在时域的计算值

分别计算两个通道数据的希尔伯特变换，然后计算数据的瞬时相位差值，然后根据所有点的相位差计算自然常数为底的指数，将指数的绝对值求和平均，得到PLV的值

值的区间是0-1，0说明相位完全没同步

Matlab函数：文末有附
https://www.mathworks.com/matlabcentral/fileexchange/31600-phase-locking-value

用法：
脑电数据格式是16个通道x840个数据点x80个trial
采样率1200Hz
滤波区间8-12Hz
绘制3通道和4通道的PLV值

  eegData=samples(:,:,find(labels==0));srate = 1200; %HzfiltSpec.order = 50;filtSpec.range = [8 12]; %Hz[plv] = pn_eegPLV(eegData, srate, filtSpec);figure;plot(squeeze(plv(:,3,4)));xlabel('Time (s)'); ylabel('Plase Locking Value');

结果：

函数：

function [plv] = pn_eegPLV(eegData, srate, filtSpec, dataSelectArr)
% Computes the Phase Locking Value (PLV) for an EEG dataset.
%
% Input parameters:
%   eegData is a 3D matrix numChannels x numTimePoints x numTrials
%   srate is the sampling rate of the EEG data
%   filtSpec is the filter specification to filter the EEG signal in the
%     desired frequency band of interest. It is a structure with two
%     fields, order and range. 
%       Range specifies the limits of the frequency
%     band, for example, put filtSpec.range = [35 45] for gamma band.
%       Specify the order of the FIR filter in filtSpec.order. A useful
%     rule of thumb can be to include about 4 to 5 cycles of the desired
%     signal. For example, filtSpec.order = 50 for eeg data sampled at
%     500 Hz corresponds to 100 ms and contains ~4 cycles of gamma band
%     (40 Hz).
%   dataSelectArr (OPTIONAL) is a logical 2D matrix of size - numTrials x
%     numConditions. For example, if you have a 250 trials in your EEG
%     dataset and the first 125 correspond to the 'attend' condition and
%     the last 125 correspond to the 'ignore' condition, then use
%     dataSelectArr = [[true(125, 1); false(125, 1)],...
%       [false(125, 1); true(125, 1)]];
%
% Output parameters:
%   plv is a 4D matrix - 
%     numTimePoints x numChannels x numChannels x numConditions
%   If 'dataSelectArr' is not specified, then it is assumed that there is
%   only one condition and all trials belong to that condition.
%
%--------------------------------------------------------------------------
% Example: Consider a 28 channel EEG data sampled @ 500 Hz with 231 trials,
% where each trial lasts for 2 seconds. You are required to plot the phase
% locking value in the gamma band between channels Fz (17) and Oz (20) for
% two conditions (say, attend and ignore). Below is an example of how to
% use this function.
%
%   eegData = rand(28, 1000, 231); 
%   srate = 500; %Hz
%   filtSpec.order = 50;
%   filtSpec.range = [35 45]; %Hz
%   dataSelectArr = rand(231, 1) >= 0.5; % attend trials
%   dataSelectArr(:, 2) = ~dataSelectArr(:, 1); % ignore trials
%   [plv] = pn_eegPLV(eegData, srate, filtSpec, dataSelectArr);
%   figure; plot((0:size(eegData, 2)-1)/srate, squeeze(plv(:, 17, 20, :)));
%   xlabel('Time (s)'); ylabel('Plase Locking Value');
%
% NOTE:
% As you have probably noticed in the plot from the above example, the PLV 
% between two random signals is spuriously large in the first 100 ms. While 
% using FIR filtering and/or hilbert transform, it is good practice to 
% discard both ends of the signal (same number of samples as the order of 
% the FIR filter, or more).
% 
% Also note that in order to extract the PLV between channels 17 and 20, 
% use plv(:, 17, 20, :) and NOT plv(:, 20, 17, :). The smaller channel 
% number is to be used first.
%--------------------------------------------------------------------------
% 
% Reference:
%   Lachaux, J P, E Rodriguez, J Martinerie, and F J Varela. 
%   �Measuring phase synchrony in brain signals.� 
%   Human brain mapping 8, no. 4 (January 1999): 194-208. 
%   http://www.ncbi.nlm.nih.gov/pubmed/10619414.
% 
%--------------------------------------------------------------------------
% Written by: 
% Praneeth Namburi
% Cognitive Neuroscience Lab, DUKE-NUS
% 01 Dec 2009
% 
% Present address: Neuroscience Graduate Program, MIT
% email:           praneeth@mit.edunumChannels = size(eegData, 1);
numTrials = size(eegData, 3);
if ~exist('dataSelectArr', 'var')dataSelectArr = true(numTrials, 1);
elseif ~islogical(dataSelectArr)error('Data selection array must be a logical');end
end
numConditions = size(dataSelectArr, 2);% disp('Filtering data...');
filtPts = fir1(filtSpec.order, 2/srate*filtSpec.range);
filteredData = filter(filtPts, 1, eegData, [], 2);% disp(['Calculating PLV for ' mat2str(sum(dataSelectArr, 1)) ' trials...']);
for channelCount = 1:numChannelsfilteredData(channelCount, :, :) = angle(hilbert(squeeze(filteredData(channelCount, :, :))));
end
plv = zeros(size(filteredData, 2), numChannels, numChannels, numConditions);
for channelCount = 1:numChannels-1channelData = squeeze(filteredData(channelCount, :, :));for compareChannelCount = channelCount+1:numChannelscompareChannelData = squeeze(filteredData(compareChannelCount, :, :));for conditionCount = 1:numConditionsplv(:, channelCount, compareChannelCount, conditionCount) = abs(sum(exp(1i*(channelData(:, dataSelectArr(:, conditionCount)) - compareChannelData(:, dataSelectArr(:, conditionCount)))), 2))/sum(dataSelectArr(:, conditionCount));endend
end
plv = squeeze(plv);
return;