从公式到代码详细解析PLV（锁相值）

1、作用：估计两个神经电或生物磁信号之间的瞬时相位关系。
2、步骤：①选择感兴趣的频段对原始信号进行带通滤波
②利用希尔伯特变换计算信号的瞬时相位,如图Fig1
③计算PLV，如图Fig2
思考为什么要分频：因为利用希尔伯特变换计算窄带信号的瞬时相位会更准确。

3、公式及解释

4、trial和electrode之间的关系如下：

5、举个例子：


6、计算自己数据的PLV—matlab代码

clc
clear
data_on = pop_loadset('epoch_65_s5.set'); % 加载自己的数据
% eegData：[62, 650, 5],62表示62导，5表示有5个trail/epoch，650表示每一个epoch的时长
eegData = data_on.data;
srate = 500; %采样频率  产生的信号的总时间长度 t = dt*采样数 = 采样数/Fs = 650/500 = 1.3s = 1300ms
filtSpec.order = 50;
filtSpec.range = [4 7]; %theta波段
%% 绘制原始信号(叠加平均)
figure;
Data_17 = squeeze(eegData(17, :, :));
Data_17 = squeeze(mean(Data_17,2));
plot((0:size(eegData, 2)-1)/srate, Data_17,'r');
hold on;
Data_20 = squeeze(eegData(20, :, :));
Data_20 = squeeze(mean(Data_20, 2));
plot((0:size(eegData, 2)-1)/srate, Data_20);
xlabel('Time (s)'); ylabel('Original Signal');
set(gca,'Box','off')
saveas(gcf, 'SIGNAL', 'png')
%% 开始计算PLV

numChannels = size(eegData, 1);
numTrials = size(eegData, 3); % 刺激数
%% 创建逻辑列
if ~exist('dataSelectArr', 'var')
    dataSelectArr = true(numTrials, 1); 
    % true(n, m)：该函数创建n*m的矩阵，该方阵的所有元素为逻辑真，即1，这里是5行1列
else
    if ~islogical(dataSelectArr) % islogical判断对象是否为逻辑类型的数据（true/false）
        error('Data selection array must be a logical');
    end
end
numConditions = size(dataSelectArr, 2); % 是指返回dataSelectArr的列数，故numConditions=1
%% 滤波
filtPts = fir1(filtSpec.order, 2/srate*filtSpec.range);  % 滤波窗函数
filteredData = filter(filtPts, 1, eegData, [], 2);
%% 计算瞬时相位
for channelCount = 1:numChannels
    % 取滤波后每一个导联的所有trail并进行希尔伯特变换，最后利用angle()函数计算相位角
    % squeeze删除单一维度，矩阵压缩
    filteredData(channelCount, :, :) = angle(hilbert(squeeze(filteredData(channelCount, :, :))));
end
%% 计算PLV
disp(size(filteredData, 2));
plv = zeros(size(filteredData, 2), numChannels, numChannels, numConditions);
for channelCount = 1:numChannels-1 % 取第一个导联所有trail的相位角
    channelData = squeeze(filteredData(channelCount, :, :)); 
   
    for compareChannelCount = channelCount+1:numChannels % 取第二个导联的相位角
        compareChannelData = squeeze(filteredData(compareChannelCount, :, :));
        for conditionCount = 1:numConditions
            % sum(A,2)按行求和,sum(dataSelectArr(:,conditionCount))=5,表示有5个trail，即：N=5
            % exp(1i*theta)返回解析函数
            plv(:, channelCount, compareChannelCount, conditionCount) = ...
            abs(sum(exp(1i*(channelData(:, dataSelectArr(:, conditionCount)) - ...
            compareChannelData(:, dataSelectArr(:, conditionCount)))), 2))/...
            sum(dataSelectArr(:, conditionCount));% 计算PLV，sum(x,2);%行求和
        end
    end
    
   %% 下面这一块代码就是先求幅度abs，再求和sum，你会发现plv都等于1
%     for compareChannelCount = channelCount+1:numChannels % 取第二个导联的相位角
%         compareChannelData = squeeze(filteredData(compareChannelCount, :, :));
%         for conditionCount = 1:numConditions
%             % sum(A,2)按行求和,sum(dataSelectArr(:,conditionCount))=5,表示有5个trail，即：N=5
%             % exp(1i*theta)返回解析函数
%             plv(:, channelCount, compareChannelCount, conditionCount) = ...
%             sum(abs(exp(1i*(channelData(:, dataSelectArr(:, conditionCount)) - ...
%             compareChannelData(:, dataSelectArr(:, conditionCount))))), 2)/...
%             sum(dataSelectArr(:, conditionCount));% 计算PLV，sum(x,2);%行求和
%         end
%     end
    
    %% 遮一部分是对上面plv公式的拆解，更好理解每一部分计算出来的结果分别表示什么
    aaaaa = channelData(:, dataSelectArr(:, 1));
    bbbbb = channelData(:, dataSelectArr(:, conditionCount)) - ...
            compareChannelData(:, dataSelectArr(:, conditionCount));
    ccccc = exp(1i*(channelData(:, dataSelectArr(:, conditionCount)) - ...
            compareChannelData(:, dataSelectArr(:, conditionCount))));
    ddddd = sum(exp(1i*(channelData(:, dataSelectArr(:, conditionCount)) - ...
            compareChannelData(:, dataSelectArr(:, conditionCount)))), 2);
    eeeee = abs(sum(exp(1i*(channelData(:, dataSelectArr(:, conditionCount)) - ...
            compareChannelData(:, dataSelectArr(:, conditionCount)))), 2));
    kkkkk = sum(dataSelectArr(:, conditionCount));
    
end
plv = squeeze(plv);
%% 绘制通道17和20之间的PLV，
figure; plot((0:size(eegData, 2)-1)/srate, squeeze(plv(:, 17, 20, :)));
xlabel('Time (s)'); ylabel('Plase Locking Value');
set(gca,'Box','off')
saveas(gcf, 'PLV', 'png')

7、核心代码解释：

8、上述计算过程如下图所示：

9、PLV代码（function）

链接：https://pan.baidu.com/s/1gYOW36tZBQdmzcBWklC6IA
提取码：0pyg