123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100 |
- % do within-stage PLI surrogates after geting bump locations and define tempral windows for connectivity
- clear all
- close all
- clc
- nSurro = 200;
- buOff = 2; % FC stges everything between bumps except for bumps
- withBL = true;
- dimBL = [10, 6]; % form 40 - 10 to onset - 6 samples -> BL with 25 samples, 250 millisecons
- dataHSMM = '/hmmBu4RroiAveMap3StagBL.mat';
- savename = '/surroCirchmmBu4RroiAveMap3StagBL.mat';
- blOff = 40; % 400 milliseconds baseline
- nSj = 18;
- fs = 100;
- nRoi = 68;
- nFC = (nRoi^2)/2 - nRoi/2;
- savePath = './dataSave';
- dataPath = './data';
- dataMEG = './data/thetaROIs_R.mat';
- load([dataPath dataHSMM], 'proBump', 'x','y')
- load(dataMEG, 'roiaveR','xR', 'yR', 'subjtR')
- nStage = size(proBump, 3)+1;
- nStageBL= nStage;
- buOff_i = buOff * ones(1, nStage-1);
- buOff_i = [0, buOff_i];
- buOff_e = buOff * ones(1, nStage-1);
- buOff_e = [buOff_e, 0];
- if withBL,
- nStageBL = nStage + 1;
- buOff_i = [dimBL(1), buOff_i];
- buOff_e = [dimBL(2), buOff_e];
- end
- roiR = single(roiaveR);
- xd = xR;
- yd = yR;
- subj = subjtR;
- clear roiaveR xR yR subjtR
- pairs = combnk([1:nRoi],2);
- tooShort = [];
- pliSurro = zeros(nStageBL, nFC, nSurro, 'single');
- for sj = 1:nSj
- xj = x(subj==sj); % bumps
- yj = y(subj==sj);
- xdj = xd(subj==sj); % data
- ydj = yd(subj==sj);
- nTr = length(xj);
- lens = yj - xj + 1;
- pliTr = zeros(nStageBL, nFC, nSurro, nTr,'single');
- maxN = max(lens);
- buSampl = ones(nTr, nStage,'single');
- % location of bumps given by HsMM
- for bu = 1:nStage-1
- buSampl(:,bu+1) = round( [1:size(proBump,1)] * squeeze(proBump(:, subj==sj, bu)));
- end
- buSampl = horzcat(buSampl, lens);
- buSampl = buSampl + blOff;
- if withBL,
- buSampl = horzcat(ones(size(buSampl,1), 1), buSampl);
- end
- for tr = 1:nTr
- % pli: corss spectral power
- data = hilbert(roiR(xdj(tr):ydj(tr),:));
- data = repmat(data,1,1,nSurro);
- dataSu = zeros(size(data));
- raIx = uint32(randi(lens(tr)-1,1,nSurro));
- for si = 1:nSurro
- dataSu(:,:,si) = circshift(data(:,:,1), raIx(si));
- end
- % corss spectral density
- csd = data(:, pairs(:,1),:) .* conj(dataSu(:, pairs(:,2),:));
- % sample phase differences sign
- siCsd = sign( imag( csd )); % [samples, FC, surros]
- % stage mean PLI
- for s = 1:size(buSampl,2)-1
- if (buSampl(tr,s+1)-buOff_e(s)) - (buSampl(tr,s)+buOff_i(s)) <= 0,
- pliTr(s,:,:,tr) = 0; % [stages, FC, surro, trials]
- tooShort = vertcat(tooShort, [s,tr]);
- end
- pliTr(s,:,:,tr) = mean( siCsd( buSampl(tr,s)+buOff_i(s):buSampl(tr,s+1)-buOff_e(s) ,:,:), 1);
- end
- end
- if any(any(any(any(isnan(pliTr)))))
- aNaN = 1;
- end
- pliSurro(:,:,:) = squeeze(pliSurro(:,:,:)) + squeeze( mean(pliTr, 4) ) ./ nSj; % mean over subjects
- end
- p95a = prctile(abs(pliSurro),95,3); % minimum FC respects to Surrogates
- p98a = prctile(abs(pliSurro),98,3); % minimum FC respects to Surrogates
- %save([savePath savename], 'pliSurro','p95a', 'p98a', 'tooShort')
-
|