相位同步指数（phase synchronization index，PSI）

＾哪来的＆永远～

已于 2023-03-28 21:50:27 修改

阅读量1.8k

点赞数 2

文章标签：人工智能算法

于 2023-03-28 21:48:28 首次发布

本文链接： https://blog.csdn.net/qq_45538220/article/details/129825070

版权

之前详细写过PLV和 PTE的原理及计算过程，发现很受大家欢迎。所以就想着后续有时间就以相同的方式（我觉得好理解）给大家继续介绍其他功能连接的算法。
这次介绍的指标是相位同步指数，介绍这个算法的目的是因为在PLV的基础上，有些人关注的不是某一时刻两个信号相位差的一致性，而是关注某一时刻两个信号的相位一致性。因此在这里给出两个信号相位一致性的原理及计算过程。其实这里我没有找到参考文献，查网上的资料也没有详细计算这个指标的过程，所以有一部分是我根据其他功能连接方法的思路去写的，如果有不对或疑惑的的地方，希望大家指出，可以私信我或在评论区评论。

一、PSI计算步骤

上述是计算PSI的过程，根据推导我觉得这个公式到这里还没有结束，因为结果不是我想要的结果，因此根据我的目的，我就想到了以下改进方法。

二、改进的PSI（暂且称为：dPSI）

代码及改进结果如下：
1、两个信号的相位差为：pi，dPSI = 0.1064

2、两个信号的相位差为：0，dPSI = 1

归一化的最大值为什么选择pi/2，最小值为什么选择0.1578？
因为当计算出来的：PSI = 0，反余切=pi/2；PSI = 2pi，反余切=0.1578。

三、dPSI用于EEG信号并绘制热力图

代码：

clc
clear
data_on = pop_loadset('epoch_65_s5.set'); % 加载自己的数据
% eegData：[62, 650, 5], 62表示62导，5表示有5个trail，650表示每一个trial的时长
eegData = data_on.data;
srate = 500; %采样频率  产生的信号的总时间长度 t = dt*采样数 = 采样数/Fs = 650/500 = 1.3s = 1300ms

%% 绘制原始信号(叠加平均)
figure;
t = (0:size(eegData, 2)-1)/srate;
channel_17 = squeeze(eegData(17, :, :));
channel_17 = squeeze(mean(channel_17,2));
plot(t, channel_17,'b','linewidth',2);
hold on;
channel_20 = squeeze(eegData(20, :, :));
channel_20 = squeeze(mean(channel_20, 2));
plot(t, channel_20,'r','linewidth',2);
xlabel('Time (s)'); ylabel('Amplitude');
title("Original Signal");
grid on

%% 计算信号瞬时相位，这里计算的是全频段的dPSI，其实最好分频，至于原因，请看我PLV那篇博客。
phaseData = Calculate_Instantaneous_Phase(eegData);

%% 计算dPSI
numChannels = size(eegData, 1);
numTrials = size(eegData, 3);      % 刺激数
Max = 1.5708;
Min = 0.1578;

dPSI = zeros(size(phaseData, 2), numChannels, numChannels);
for channelCount = 1:numChannels-1 % 取第一个导联所有trail的相位角
    channelData = squeeze(phaseData(channelCount, :, :));   % [650,5]
   
    for compareChannelCount = channelCount+1:numChannels % 取第二个导联的相位角
        compareChannelData = squeeze(phaseData(compareChannelCount, :, :)); % [650,5]
        
        % 计算相位差
        dif_PSI = abs(sum(channelData - compareChannelData, 2))/numTrials;  % [650,1]
        % 获取反余切值
        acot_PSI = acot(dif_PSI);
        % 归一化，对角矩阵
        dPSI(:,channelCount,compareChannelCount) = round((acot_PSI-Min)/(Max-Min),2);  % 归一化并取4位小数
        dPSI(:,compareChannelCount,channelCount) = round((acot_PSI-Min)/(Max-Min),2);  % 归一化并取4位小数
    end
end

%% 绘制通道17和20之间的dPSI，
figure; 
plot(t, dPSI(:, 17, 20),'b','linewidth',2);
xlabel('Time (s)'); ylabel('dPSI');
title("dPSI of X1 and X2");
grid on

%% 绘制热力图，先对所有时间点求均值
figure
mean_dPSI = squeeze(mean(dPSI,1));  % 这里绘制的是整个时间段的dPSI，可以自行选择时间区间
h = heatmap(mean_dPSI);
colormap('jet');

 %% ---函数功能:计算信号的瞬时相位
function caculate_phase = Calculate_Instantaneous_Phase(data)
    numChannels = size(data,1);
    caculate_phase = zeros(size(data));
    for channelCount = 1:numChannels
        % 取每一个导联的所有trail并进行希尔伯特变换，最后利用angle()函数计算相位角
        % squeeze删除单一维度，矩阵压缩
        caculate_phase(channelCount, :, :) = angle(hilbert(squeeze(data(channelCount, :, :))));
    end
end