matlab实现试验间相位一致性ITPC

试验间相位一致性（inter-trial phase coherence，ITPC） 是衡量在多个试验中，特定频率成分的相位一致性大小的指标。在EEG或其他神经生理信号分析中，ITPC通常用来评估在一系列试验（如认知任务或刺激反应）中，特定脑区或脑区间的相位同步程度。

ITPC的含义：
1. 相位信息：在信号处理中，一个正弦波信号可以表示为振幅和相位的组合。相位表示波形在时间 上的位置，反映了波的周期性变化；
2. 试验间比较：在多个试验中，如果相同频率成分的相位在时间上保持一致，那么这些试验间的相位被认为是一致的。ITPC通过计算这些相位之间的一致性来量化这种同步；
3. 一致性度量：ITPC的值介于0和1之间。值越接近1表示试验间的相位高度一致，即在多次试验中，特定频率的相位在时间点上表现出相似的波形位置，值接近0表示相位一致性差，即相位在不同试验中随机分布；
4. 神经科学应用：在神经科学研究中，ITPC可以用来评估大脑在执行特定任务或特定刺激反应时的神经震荡同步性。例如，研究者可能会使用ITPC来分析注意力任务中，大脑不同区域的beta波（约13-30Hz）相位是否同步，从而了解这些区域是如何协同工作的；
5. 临床应用：在临床科学中，ITPC也可用来评估脑损伤、意识障碍或其他神经疾病患者的脑功能状态，通过分析特定频带的脑波相位一致性来辅助诊断和治疗评估。

计算ITPC：
1. 傅里叶变换：对每个试验的信号进行傅里叶变换，提取特定频率成分的相位信息；
2. 相位提取：从傅里叶变换的结果中提取目标频率的相位角；
3. 平均相位计算：计算所有试验在目标频率下的相位的平均值；
4. 相位一致性计算：计算每个试验相位与平均相位之间的差异，并使用这些差异来计算ITPC。

a）下面先写有一段数据，然后进行的计算；
b）最后是，写成函数，可直接调用的那种；

a) 数据，计算：

// %ITPC试验间相位一致性
clc;clear;close all;
a = load('EEG111802.mat'); % 204x12000x4,204通道x12000时长x4个trial
data = a.EEG.data;   %取出数据
samples = 200;  % 数据的采样率
channels = size(data,1);  %通道数
lengths = size(data,2);   %采样点数
trials = size(data,3);   %trial个数
f = 5;   % 给出要计算的目标频率

% 计算目标频率对应的索引
f_index = round(f*lengths/samples) + 1;
% 初始化ITPC值
itpc = zeros(1,channels);

% 计算每个通道的相位
for ch = 1:channels
	channelData = squeeze(data(ch,:,:)); % squeeze降维
	
	%初始化
	phase = zeros(1,trials);		
	% 对每个试验计算傅里叶变换
	for t = 1:trials
		% 提取单个试验的数据
		trialData = channelData(:,t); 
		ft = fft(trialData); % 快速傅里叶变换
		phase(t) = angle(ft(f_index));  % 相位
	end
	meanP = mean(phase);  % 计算trial均值
	dev = phase - meanP;  % 计算相位差
	itpc(ch) = abs(mean(exp(1i*dev))); %复数1i，数字1；计算复数平均值并取绝对值
end
%% 计算出来的ITPC是1x204的数值：

% 根据需要可以计算整体的一个均值
itpc_mean = mean(itpc);
% 或者也可以使用itpc画脑地形图
q = [0 1];  % 显示的范围颜色
topoplot(itpc, a.EEG.chanlocs, 'style','both','electrodes','on','numcontour',6,'maplimits',q ,'shading','interp');%
title('ITPC','fontsize',13);

b）函数，调用

% data:要计算的数据，1行是通道，2列是时间点，3分段trial
% f：目标频率（Hz）
% samples:采样率（Hz）
function itpc = calculateITPC(data, f,samples)
	channels = size(data,1);  %通道数
	lengths = size(data,2);   %采样点数
	trials = size(data,3);   %trial个数
	% 计算目标频率对应的索引
    f_index = round(f*lengths / samples) + 1; 
    
    % 初始化
    itpc = zeros(1,channels );
	phase = zeros(1,trials);
    % 计算每个通道的相位
    for ch = 1:channels 
        channelData = squeeze(data(ch,:,:));
        % 对每个试验计算傅里叶变换
        for t = 1:trials 
            % 提取单个试验的数据
            trialData = channelData(:,t);
            % trialData = trialData - mean(trialData); %去直流分量
            ft = fft(trialData); % 快速傅里叶变换
            phase(t) = angle(ft(f_index));  % 相位
        end
        meanP = mean(phase);
        dev = phase-meanP;
        itpc(ch) = abs(mean(exp(1i*dev))); %复数1i，数字1；计算复数平均值并取绝对值              
    end
end

%% 调用时，如下：
只是简单画了个图，根据需要画不同的图，以及后续的统计分析

如要计算30Hz：
f = 30;
samples = 200;
a = load('EEG111802.mat');
data = a.EEG.data(:,:,:);
itpc_values = calculateITPC(data,f,samples);
figure;
plot(itpc_values);
title('Inter_Trial Phase Coherence at 13Hz');
xlabel('channels Point');
ylabel('ITPC Values');