Skip to content

Commit

Permalink
Merge pull request #85 from Wirkungstreffer/BS_dev
Browse files Browse the repository at this point in the history 
Implementation of parallel toolbox
  • Loading branch information
@nguyen-td
nguyen-td authored Mar 14, 2024
2 parents 4a6d1df + 046ce4a commit d15ff51
Show file tree
Hide file tree
Showing 5 changed files with 140 additions and 113 deletions.
88 changes: 43 additions & 45 deletions libs/pellegrini/fp_data2bs_event_uni.m
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [cs,nave]=fp_data2bs_event_uni(data,segleng,segshift,epleng,freqpairs,nshuf)
function [cs,nave]=fp_data2bs_event_uni(data,segleng,segshift,epleng,freqpairs,ishuf)
% Calculates bispectral-tensors from data for event-related measurement
% and their null distributions using a shuffling approach.
%
Expand All @@ -15,7 +15,7 @@
% e.g. segshift=segleng/2 makes overlapping segments
% epleng: leng of each epoch
% freqpairs: pairs of frequency in bins
% nshuf: required number of samples (shuffles) in the null distribution
% ishuf: one shuffle in the null distribution
% para: structure which is eventually used later
%
% output:
Expand All @@ -41,54 +41,52 @@
nseg=floor((epleng-segleng)/segshift)+1; %total number of segments
assert(nseg==1,'only possible with 1 segment')

cs=zeros(nchan,nchan,nchan,2,nshuf);
cs=zeros(nchan,nchan,nchan,2);

%get Fourier coefficients
coeffs = fp_fft_coeffs(data,segleng,segshift,epleng,freqpairs);

for ishuf = 1:nshuf
nave=0;
csloc1=zeros(nchan,nchan,nchan);
csloc2=zeros(nchan,nchan,nchan);
cs1=zeros(nchan,nchan,nchan);
cs2=zeros(nchan,nchan,nchan);
%pre-allocation
nave=0;
csloc1=zeros(nchan,nchan,nchan);
csloc2=zeros(nchan,nchan,nchan);
cs1=zeros(nchan,nchan,nchan);
cs2=zeros(nchan,nchan,nchan);

if ishuf == 1 %the first shuffle contains the true order
inds = 1:nep;
else
inds = randperm(nep,nep); %indices for shuffling of epochs for channel 2
end
if ishuf == 1 %the first shuffle contains the true order
inds = 1:nep;
else
inds = randperm(nep,nep); %indices for shuffling of epochs for channel 2
end

for j=1:nep %loop over epochs

%bispec of f1 (low peak), f2-f1 (left side lobe), f2 (high peak)
csloc1(1,1,1)=transpose(coeffs(1,1,j)) *coeffs(2,1,j) *conj(coeffs(3,1,j));
csloc1(2,1,1)=transpose(coeffs(1,2,inds(j))) *coeffs(2,1,j) *conj(coeffs(3,1,j));
csloc1(1,2,1)=transpose(coeffs(1,1,j)) *coeffs(2,2,inds(j))*conj(coeffs(3,1,j));
csloc1(1,1,2)=transpose(coeffs(1,1,j)) *coeffs(2,1,j) *conj(coeffs(3,2,inds(j)));
csloc1(2,2,1)=transpose(coeffs(1,2,inds(j))) *coeffs(2,2,inds(j))*conj(coeffs(3,1,j));
csloc1(1,2,2)=transpose(coeffs(1,1,j)) *coeffs(2,2,inds(j))*conj(coeffs(3,2,inds(j)));
csloc1(2,2,2)=transpose(coeffs(1,2,inds(j))) *coeffs(2,2,inds(j))*conj(coeffs(3,2,inds(j)));
csloc1(2,1,2)=transpose(coeffs(1,2,inds(j))) *coeffs(2,1,j) *conj(coeffs(3,2,inds(j)));

for j=1:nep %loop over epochs

%bispec of f1 (low peak), f2-f1 (left side lobe), f2 (high peak)
csloc1(1,1,1)=transpose(coeffs(1,1,j)) *coeffs(2,1,j) *conj(coeffs(3,1,j));
csloc1(2,1,1)=transpose(coeffs(1,2,inds(j))) *coeffs(2,1,j) *conj(coeffs(3,1,j));
csloc1(1,2,1)=transpose(coeffs(1,1,j)) *coeffs(2,2,inds(j))*conj(coeffs(3,1,j));
csloc1(1,1,2)=transpose(coeffs(1,1,j)) *coeffs(2,1,j) *conj(coeffs(3,2,inds(j)));
csloc1(2,2,1)=transpose(coeffs(1,2,inds(j))) *coeffs(2,2,inds(j))*conj(coeffs(3,1,j));
csloc1(1,2,2)=transpose(coeffs(1,1,j)) *coeffs(2,2,inds(j))*conj(coeffs(3,2,inds(j)));
csloc1(2,2,2)=transpose(coeffs(1,2,inds(j))) *coeffs(2,2,inds(j))*conj(coeffs(3,2,inds(j)));
csloc1(2,1,2)=transpose(coeffs(1,2,inds(j))) *coeffs(2,1,j) *conj(coeffs(3,2,inds(j)));

%bispec of f1 (low peak), f2 high peak), f1+f2 (right side lobe)
csloc2(1,1,1)=transpose(coeffs(1,1,j)) *coeffs(3,1,j) *conj(coeffs(4,1,j));
csloc2(2,1,1)=transpose(coeffs(1,2,inds(j))) *coeffs(3,1,j) *conj(coeffs(4,1,j));
csloc2(1,2,1)=transpose(coeffs(1,1,j)) *coeffs(3,2,inds(j))*conj(coeffs(4,1,j));
csloc2(1,1,2)=transpose(coeffs(1,1,j)) *coeffs(3,1,j) *conj(coeffs(4,2,inds(j)));
csloc2(2,2,1)=transpose(coeffs(1,2,inds(j))) *coeffs(3,2,inds(j))*conj(coeffs(4,1,j));
csloc2(1,2,2)=transpose(coeffs(1,1,j)) *coeffs(3,2,inds(j))*conj(coeffs(4,2,inds(j)));
csloc2(2,2,2)=transpose(coeffs(1,2,inds(j))) *coeffs(3,2,inds(j))*conj(coeffs(4,2,inds(j)));
csloc2(2,1,2)=transpose(coeffs(1,2,inds(j))) *coeffs(3,1,j) *conj(coeffs(4,2,inds(j)));

cs1=cs1+csloc1;
cs2=cs2+csloc2;

nave=nave+1;
end
%bispec of f1 (low peak), f2 high peak), f1+f2 (right side lobe)
csloc2(1,1,1)=transpose(coeffs(1,1,j)) *coeffs(3,1,j) *conj(coeffs(4,1,j));
csloc2(2,1,1)=transpose(coeffs(1,2,inds(j))) *coeffs(3,1,j) *conj(coeffs(4,1,j));
csloc2(1,2,1)=transpose(coeffs(1,1,j)) *coeffs(3,2,inds(j))*conj(coeffs(4,1,j));
csloc2(1,1,2)=transpose(coeffs(1,1,j)) *coeffs(3,1,j) *conj(coeffs(4,2,inds(j)));
csloc2(2,2,1)=transpose(coeffs(1,2,inds(j))) *coeffs(3,2,inds(j))*conj(coeffs(4,1,j));
csloc2(1,2,2)=transpose(coeffs(1,1,j)) *coeffs(3,2,inds(j))*conj(coeffs(4,2,inds(j)));
csloc2(2,2,2)=transpose(coeffs(1,2,inds(j))) *coeffs(3,2,inds(j))*conj(coeffs(4,2,inds(j)));
csloc2(2,1,2)=transpose(coeffs(1,2,inds(j))) *coeffs(3,1,j) *conj(coeffs(4,2,inds(j)));

%shape cs: chan x chan x chan x peak_combination x shuffles
cs(:,:,:,1,ishuf) = cs1./nave;
cs(:,:,:,2,ishuf) = cs2./nave;
cs1=cs1+csloc1;
cs2=cs2+csloc2;

end
nave=nave+1;
end

%shape cs: chan x chan x chan x peak_combination x shuffle
cs(:,:,:,1) = cs1./nave;
cs(:,:,:,2) = cs2./nave;
8 changes: 5 additions & 3 deletions test_pipes/test_mim_stats.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,13 +22,15 @@
EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel);

%% Create null distribution
EEG1 = pop_roi_connect(EEG, 'methods', {'CS' 'cCOH', 'wPLI' 'MIM'}, 'freqresolution', 200, 'roi_selection', {1, 3, 7}, 'conn_stats', 'on', 'nshuf', 1001); % takes very long!
pop_roi_statsplot(EEG1, 'measure', 'MIM', 'freqrange', [8 13]);
EEG1_MIM = pop_roi_connect(EEG, 'methods', {'CS' 'cCOH', 'wPLI' 'MIM'}, 'freqresolution', 200, 'roi_selection', {1, 3, 7}, 'conn_stats', 'on', 'nshuf', 1000, 'poolsize', 12);
EEG2_MIM = pop_roi_connect(EEG, 'methods', {'CS' 'cCOH', 'wPLI' 'MIM'}, 'conn_stats', 'on', 'nshuf', 1000); % takes very long!
% load('test_pipes/MIM_shuf.mat')
pop_roi_statsplot(EEG1_MIM, 'measure', 'MIM', 'freqrange', [8 13]);

%% Unit test
% Test if the first shuffle of the surrogate connectivity matrix in EEG1 is in fact the true matrix
EEG1 = pop_roi_connect(EEG, 'methods', {'MIM'}, 'conn_stats', 'on', 'nshuf', 3);
EEG2 = pop_roi_connect(EEG, 'methods', {'MIM'}, 'conn_stats', 'off');
if ~isequal(squeeze(EEG1.roi.MIM(:, :, :, 1)), EEG2.roi.MIM)
error 'The first shuffle in the surrogate connectivity array is not the true matrix.'
end
end
31 changes: 22 additions & 9 deletions test_pipes/test_pac.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,11 +32,20 @@
fcomb.high = high;


%EEG1 = pop_roi_connect(EEG, 'methods', {'PAC', 'MIM', 'COH'}, 'fcomb', fcomb); % test all 3 connectivity functions (data2spwctrgc, data2strgcmim, roi_pac)
EEG2 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'bs_outopts', 5, 'conn_stats', 'on', 'nshuf', 4); % compute only b_anti, b_anti_norm
EEG1 = pop_roi_connect(EEG, 'methods', {'PAC', 'MIM', 'COH'}, 'fcomb', fcomb); % test all 3 connectivity functions (data2spwctrgc, data2strgcmim, roi_pac)
tic
EEG3 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'bs_outopts', 5); % compute only b_anti, b_anti_norm
EEG2 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'bs_outopts', 1, 'conn_stats', 'on', 'nshuf', 3, 'poolsize', 12); % compute only b_anti, b_anti_norm
toc
EEG3 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'bs_outopts', 1); % compute only b_anti, b_anti_norm

if ~isequal(squeeze(EEG2.roi.PAC.b_orig(:, :, 1)), EEG3.roi.PAC.b_orig)
error 'The first shuffle in the surrogate connectivity array is not the true matrix.'
end

if ~isequal(squeeze(EEG2.roi.PAC.b_anti(:, :, 1)), EEG3.roi.PAC.b_anti)
error 'The first shuffle in the surrogate connectivity array is not the true matrix.'
end


%% Test bispectrum for frequency band inputs
low = [4 8];
Expand All @@ -49,21 +58,25 @@
EEG4 = pop_roi_connect(EEG, 'methods', {'PAC', 'MIM', 'COH'}, 'fcomb', fcomb); % test all 3 connectivity functions (data2spwctrgc, data2strgcmim, roi_pac)toc
toc
tic
EEG5 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'conn_stats', 'off', 'nshuf', 2, 'bs_outopts', 5);
EEG5 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'conn_stats', 'off', 'nshuf', 3, 'bs_outopts', 1);
toc

if ~isequal(squeeze(EEG5.roi.PAC.b_anti(:, :, 1)), EEG4.roi.PAC.b_anti)
error 'The first shuffle in the surrogate connectivity array is not the true matrix.'
end

%% Test PAC plotting
% Test for single frequency inputs
pop_roi_connectplot(EEG1, 'measure', 'pac', 'plotmatrix', 'on');
pop_roi_connectplot(EEG1, 'measure', 'pac_anti', 'plotmatrix', 'on');
pop_roi_connectplot(EEG3, 'measure', 'pac', 'plotmatrix', 'on');
pop_roi_connectplot(EEG3, 'measure', 'pac_anti', 'plotmatrix', 'on');

% Provoke errors by plotting bispectral tensors that do not exist
pop_roi_connectplot(EEG2, 'measure', 'pac_anti', 'plotmatrix', 'on'); % bs_outopts 4 means only original bispectra are computed, so cannot plot anti
pop_roi_connectplot(EEG3, 'measure', 'pac_anti', 'plotmatrix', 'on'); % bs_outopts 4 means only original bispectra are computed, so cannot plot anti
pop_roi_connectplot(EEG3, 'measure', 'pac', 'plotmatrix', 'on'); % bs_outopts 5 means only antisymm. bispectra are computed, so cannot plot normal bispectrum

% Test for frequency bands
pop_roi_connectplot(EEG4, 'measure', 'pac', 'plotmatrix', 'on');
pop_roi_connectplot(EEG4, 'measure', 'pac_anti', 'plotmatrix', 'on');
pop_roi_connectplot(EEG5, 'measure', 'pac', 'plotmatrix', 'on');
pop_roi_connectplot(EEG5, 'measure', 'pac_anti', 'plotmatrix', 'on');

% plot MIM and COH as a sanity check
pop_roi_connectplot(EEG1, 'measure', 'mim', 'plotmatrix', 'on');
Expand Down
18 changes: 9 additions & 9 deletions utils/calc_pac.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,14 @@
% computation of PAC_km(f1, f2) = 0.5 * |Bkmm(f1, f2-f1)| + 0.5 * |Bkmm(f1, f2)|.
%
% Inputs:
% BS - [array] (nchan x nchan x nchan x nfreqpairs x nshuf) Cross-bispectral tensor
% BS - [array] (nchan x nchan x nchan x nfreqpairs) Cross-bispectral tensor
% RTP - [array] (nchan x nchan x nchan) Threenorm tensor to normalize the cross-bispectrum, used to compute cross-bicoherence
%
% Outputs:
% biv_orig - [array] (nchan x nshuf) Absolute value of the original cross-bispectrum, corresponding to [Bkmm(f1, f2), Bmkk(f1, f2)]
% biv_anti - [array] (nchan x nshuf) Absolute value of the antisymmetrized cross-bispectrum
% biv_orig_norm - [array] (nchan x nshuf) Absolute value of the normalized original cross-bispectrum (cross-bicoherence)
% biv_anti_norm - [array] (nchan x nshuf) Absolute value of the normalized antisymmetrized cross-bispectrum (antisymmetrized cross-bicoherence)
% biv_orig - [array] (1 x nchan) Absolute value of the original cross-bispectrum, corresponding to [Bkmm(f1, f2), Bmkk(f1, f2)]
% biv_anti - [array] (1 x nchan) Absolute value of the antisymmetrized cross-bispectrum
% biv_orig_norm - [array] (1 x nchan) Absolute value of the normalized original cross-bispectrum (cross-bicoherence)
% biv_anti_norm - [array] (1 x nchan) Absolute value of the normalized antisymmetrized cross-bispectrum (antisymmetrized cross-bicoherence)
%
% Authors:
% Tien Dung Nguyen, [email protected],
Expand All @@ -19,14 +19,14 @@
function [biv_orig, biv_anti, biv_orig_norm, biv_anti_norm] = calc_pac(BS, RTP)

% Calculate absolute values of the cross-bispectrum
biv_orig = squeeze(([abs(mean(BS(1, 2, 2, :, :), 4)) abs(mean(BS(2, 1, 1, :, :), 4))])); % [Bkmm, Bmkk], average over frequency bands (4th dimension)
biv_orig = squeeze(([abs(mean(BS(1, 2, 2, :), 4)) abs(mean(BS(2, 1, 1, :), 4))])); % [Bkmm, Bmkk], average over frequency bands (4th dimension)
xx = BS - permute(BS, [2 1 3 4 5]);
biv_anti = squeeze(([abs(mean(xx(1, 2, 2, :,:), 4)) abs(mean(xx(2, 1, 1, :, :), 4))]));
biv_anti = squeeze(([abs(mean(xx(1, 2, 2, :), 4)) abs(mean(xx(2, 1, 1, :), 4))]));

% Calculate absolute values of the cross-bicoherence
bicoh = BS ./ RTP;
biv_orig_norm = squeeze(([abs(mean(bicoh(1, 2, 2, :, :), 4)) abs(mean(bicoh(2, 1, 1, :, :), 4))])); % [Bkmm, Bmkk], average over frequency bands
biv_orig_norm = squeeze(([abs(mean(bicoh(1, 2, 2, :), 4)) abs(mean(bicoh(2, 1, 1, :), 4))])); % [Bkmm, Bmkk], average over frequency bands
xx = bicoh - permute(bicoh, [2 1 3 4 5]);
biv_anti_norm = squeeze(([abs(mean(xx(1, 2, 2, :, :), 4)) abs(mean(xx(2, 1, 1,:, :), 4))]));
biv_anti_norm = squeeze(([abs(mean(xx(1, 2, 2, :), 4)) abs(mean(xx(2, 1, 1,:), 4))]));

end
108 changes: 61 additions & 47 deletions utils/shuffle_BS.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@

[nchan, ndat, nepo] = size(data);

% [inds, PCA_inds] = fp_npcs2inds(npcs);
% [inds, PCA_inds] = fp_npcs2inds(npcs);
if isempty(g.roi_selection)
nROI = nchan/npcs(1);
else
Expand Down Expand Up @@ -86,21 +86,6 @@
fprintf('Generating null distribution using %d shuffles...\n', nshuf)
fprintf('Progress of %d:\n', nshuf);


% Check if Parallel Processing Toolbox is available and licensed
if license('test', 'Distrib_Computing_Toolbox') && ~isempty(ver('parallel'))
% If g.poolsize is defined, create a parallel pool of that size
if isfield(g, 'poolsize') && isnumeric(g.poolsize) && g.poolsize > 0
% Check if there's already an existing parallel pool
currentPool = gcp('nocreate');
if isempty(currentPool)
parpool(g.poolsize);
end
end
else
disp('Parallel Processing Toolbox is not installed or licensed.');
end

% Define bispectrum parameters
fres = fs;
frqs = sfreqs(fres, fs);
Expand Down Expand Up @@ -132,41 +117,70 @@
freqinds_up(i,:) = [find(frqs == low) find(frqs == high)];
end

% Initialize variables to store the PAC results
PAC_orig = zeros(nROI, nROI, nshuf);
PAC_anti = zeros(nROI, nROI, nshuf);
PAC_orig_norm = zeros(nROI, nROI, nshuf);
PAC_anti_norm = zeros(nROI, nROI, nshuf);

% Iterate over ROI pairs
for proi = 1:nROI
for aroi = proi:nROI
% Compute bispectrum % nchan by nchan by nchan by number_of_peaks by number_of_shuffles
X = data([proi aroi],:,:); % number of regions x epoch length x trials
[BS, ~] = fp_data2bs_event_uni(X(:, :)', ndat, floor(ndat/2), ndat, freqinds_up, nshuf); % pass (f1,f2) through freqinds_up
BS_low = BS(:,:,:,1,:);
BS_up = BS(:,:,:,2,:);

[RTP_low,~] = data2bs_threenorm(X(:, :)', ndat, floor(ndat/2), ndat, freqinds_low);
[RTP_up,~] = data2bs_threenorm(X(:, :)', ndat, floor(ndat/2), ndat, freqinds_up);

[biv_orig_low, biv_anti_low, biv_orig_low_norm, biv_anti_low_norm] = calc_pac(BS_low, RTP_low);
[biv_orig_up, biv_anti_up, biv_orig_up_norm, biv_anti_up_norm] = calc_pac(BS_up, RTP_up);

% PAC_km(f1, f2) = 0.5 * |Bkmm(f1, f2-f1)| + 0.5 * |Bkmm(f1, f2)|
b_orig(aroi, proi, :) = mean(vertcat(biv_orig_up(1, :), biv_orig_low(1, :)), 1);
b_orig(proi, aroi, :) = mean(vertcat(biv_orig_up(2, :), biv_orig_low(2, :)), 1);
b_anti(aroi, proi, :) = mean(vertcat(biv_anti_up(1, :), biv_anti_low(1, :)), 1);
b_anti(proi, aroi, :) = mean(vertcat(biv_anti_up(2, :), biv_anti_low(2, :)), 1);
% check if Parallel Processing Toolbox is available and licensed
if license('test', 'Distrib_Computing_Toolbox') && ~isempty(ver('parallel'))
if isfield(g, 'poolsize') && isnumeric(g.poolsize) && g.poolsize > 0
% check if there's already an existing parallel pool
currentPool = gcp('nocreate');
if isempty(currentPool)
parpool(g.poolsize);
end
end
else
disp('Parallel Processing Toolbox is not installed or licensed.');
end

% Pre-allocate variables outside the parfor loop
b_orig = zeros(nROI, nROI, nshuf);
b_anti = zeros(nROI, nROI, nshuf);
b_orig_norm = zeros(nROI, nROI, nshuf);
b_anti_norm = zeros(nROI, nROI, nshuf);

for ishuf = 1:nshuf
b_orig_ishuf = zeros(nROI, nROI);
b_anti_ishuf = zeros(nROI, nROI);
b_orig_norm_ishuf = zeros(nROI, nROI);
b_anti_norm_ishuf = zeros(nROI, nROI);

% Iterate over ROI pairs
for proi = 1:nROI
for aroi = proi:nROI
% Compute bispectrum % nchan by nchan by nchan by number_of_peaks by number_of_shuffles
X = data([proi aroi],:,:); % number of regions x epoch length x ishuf
[BS, ~] = fp_data2bs_event_uni(X(:, :)', ndat, floor(ndat/2), ndat, freqinds_up, ishuf); % pass (f1,f2) through freqinds_up
BS_low = BS(:,:,:,1);
BS_up = BS(:,:,:,2);

[RTP_low,~] = data2bs_threenorm(X(:, :)', ndat, floor(ndat/2), ndat, freqinds_low);
[RTP_up,~] = data2bs_threenorm(X(:, :)', ndat, floor(ndat/2), ndat, freqinds_up);

[biv_orig_low, biv_anti_low, biv_orig_low_norm, biv_anti_low_norm] = calc_pac(BS_low, RTP_low);
[biv_orig_up, biv_anti_up, biv_orig_up_norm, biv_anti_up_norm] = calc_pac(BS_up, RTP_up);

% PAC_km(f1, f2) = 0.5 * |Bkmm(f1, f2-f1)| + 0.5 * |Bkmm(f1, f2)|
b_orig_ishuf(aroi,proi) = mean([biv_orig_up(1) biv_orig_low(1)]);
b_orig_ishuf(proi,aroi) = mean([biv_orig_up(2) biv_orig_low(2)]);
b_anti_ishuf(aroi,proi) = mean([biv_anti_up(1) biv_anti_low(1)]);
b_anti_ishuf(proi,aroi) = mean([biv_anti_up(2) biv_anti_low(2)]);

% normalized versions (for bicoherence)
b_orig_norm(aroi,proi) = mean([biv_orig_up_norm(1) biv_orig_low_norm(1)]);
b_orig_norm(proi,aroi) = mean([biv_orig_up_norm(2) biv_orig_low_norm(2)]);
b_anti_norm(aroi,proi) = mean([biv_anti_up_norm(1) biv_anti_low_norm(1)]);
b_anti_norm(proi,aroi) = mean([biv_anti_up_norm(2) biv_anti_low_norm(2)]);
end

% normalized versions (for bicoherence)
b_orig_norm(aroi, proi, :) = mean(vertcat(biv_orig_up_norm(1, :), biv_orig_low_norm(1, :)), 1);
b_orig_norm(proi, aroi, :) = mean(vertcat(biv_orig_up_norm(2, :), biv_orig_low_norm(2, :)), 1);
b_anti_norm(aroi, proi, :) = mean(vertcat(biv_anti_up_norm(1, :), biv_anti_low_norm(1, :)), 1);
b_anti_norm(proi, aroi, :) = mean(vertcat(biv_anti_up_norm(2, :), biv_anti_low_norm(2, :)), 1);
end

% Store the shuffle information
b_orig(:,:, ishuf) = b_orig_ishuf;
b_anti(:,:, ishuf) = b_anti_ishuf;
b_orig_norm(:,:, ishuf) = b_orig_norm_ishuf;
b_anti_norm(:,:, ishuf) = b_anti_norm_ishuf;

end


clear out

% Save PAC results in the output structure
Expand Down

0 comments on commit d15ff51

Please sign in to comment.