nfft = length(data1);
fs = 500;
% [2] psd
overlap = nfft/2;
window = hann(nfft);
[Pxx1, f1] = pwelch(data1, window, overlap, nfft, fs);
[Pxx2, f2] = pwelch(data2, window, overlap, nfft, fs);
data_psd = [Pxx1(2:nfft/500*45+1), Pxx2(2:nfft/500*45+1)];
[icc_psd, LB, UB, F, df1, df2, p] = ICC(data_psd, 'A-1', 0.05, 0);
ICC_psd_all = [ICC_psd_all; icc_psd];
ICC_wave_matrix(event_index,channel_index) = icc_wave;
ICC_PSD_matrix(event_index,channel_index) = icc_psd;
% [3] 节律
% 从psd中导出
% delta 0.5-4Hz Theta 4-8Hz Alpha 8-14Hz Beta 14-30Hz Gamma 30-80Hz
[icc_delta, icc_theta, icc_alpha, icc_beta, icc_Gamma] = compute_rhythm2 (Pxx1, f1, Pxx2, f2);
ICC_delta_matrix(event_index, channel_index) = icc_delta;
ICC_theta_matrix(event_index, channel_index) = icc_theta;
ICC_alpha_matrix(event_index, channel_index) = icc_alpha;
ICC_beta_matrix(event_index, channel_index) = icc_beta;
ICC_Gamma_matrix(event_index, channel_index) = icc_Gamma;
function [icc_delta, icc_theta, icc_alpha, icc_beta, icc_Gamma] = compute_rhythm2 (Pxx1, f1, Pxx2, f2)
% delta 0.5-4Hz
index1_delta = find(f1==0);index2_delta = find(f1==4);
data1_delta = Pxx1(index1_delta:index2_delta);
index3_delta = find(f2==0);index4_delta = find(f2==4);
data2_delta = Pxx2(index3_delta:index4_delta);
data_delta = [data1_delta, data2_delta];
% Theta 4-8Hz
index1_Theta = find(f1==4);index2_Theta = find(f1==8);
data1_theta = Pxx1(index1_Theta:index2_Theta);
index3_Theta = find(f2==4);index4_Theta = find(f2==8);
data2_theta = Pxx2(index3_Theta:index4_Theta);
data_theta = [data1_theta, data2_theta];
% Alpha 8-14Hz
index1_alpha = find(f1==8);index2_alpha = find(f1==14);
data1_alpha = Pxx1(index1_alpha:index2_alpha);
index3_alpha = find(f2==8);index4_alpha = find(f2==14);
data2_alpha = Pxx2(index3_alpha:index4_alpha);
data_alpha = [data1_alpha, data2_alpha];
% Beta 14-30Hz
index1_beta = find(f1==14);index2_beta = find(f1==30);
data1_beta = Pxx1(index1_beta:index2_beta);
index3_beta = find(f2==14);index4_beta = find(f2==30);
data2_beta = Pxx2(index3_beta:index4_beta);
data_beta = [data1_beta, data2_beta];
% Gamma 30-80Hz
index1_Gamma = find(f1==30);index2_Gamma = find(f1==80);
data1_Gamma = Pxx1(index1_Gamma:index2_Gamma);
index3_Gamma = find(f2==30);index4_Gamma = find(f2==80);
data2_Gamma = Pxx2(index3_Gamma:index4_Gamma);
data_Gamma = [data1_Gamma, data2_Gamma];
[icc_delta, LB, UB, F, df1, df2, p] = ICC(data_delta, 'A-1', 0.05, 0);
[icc_theta, LB, UB, F, df1, df2, p] = ICC(data_theta, 'A-1', 0.05, 0);
[icc_alpha, LB, UB, F, df1, df2, p] = ICC(data_alpha, 'A-1', 0.05, 0);
[icc_beta, LB, UB, F, df1, df2, p] = ICC(data_beta, 'A-1', 0.05, 0);
[icc_Gamma, LB, UB, F, df1, df2, p] = ICC(data_Gamma, 'A-1', 0.05, 0);
end