diff --git a/libs/nolte/data2bs_event.m b/libs/nolte/data2bs_event.m index e5a476c..fbe2563 100755 --- a/libs/nolte/data2bs_event.m +++ b/libs/nolte/data2bs_event.m @@ -17,7 +17,7 @@ % cs: nchan by nchan by nchan by number_of_frequency_pairs % (by number_of_segments) tensor such that % cs(i,j,k,f)= -% where f1=freqpairs(f,1) and f2=freqpairs(f,1), +% where f1=freqpairs(f,1) and f2=freqpairs(f,2), % x=fft(data) and the average is over epeochs and segments % % if para.fave=0 then cs contains a fifth argument denoting diff --git a/libs/pellegrini/fp_cpsd_mt.m b/libs/pellegrini/fp_cpsd_mt.m deleted file mode 100644 index 3973147..0000000 --- a/libs/pellegrini/fp_cpsd_mt.m +++ /dev/null @@ -1,61 +0,0 @@ -function S = fp_cpsd_mt(X1,X2,ind_1,ind_2,h,window,noverlap,nchunks,taparray) - -% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe - -ind_pow = intersect(ind_1, ind_2); -nfft = 2*(h-1); - -n1 = numel(ind_1); -n2 = numel(ind_2); -n3 = numel(ind_pow); -S = complex(zeros(h,n1,n2)); -pow = zeros(h,n3); - -winstep = window-noverlap; -ntapers = size(taparray,2); - -% compute tapered periodogram with FFT - -for k = 1:nchunks - - XSEG1 = X1((1:window) + (k-1)*winstep,:); - XSEG2 = X2((1:window) + (k-1)*winstep,:); - - % compute periodogram - P1 = fft(taparray.*permute(XSEG1(:,:,ones(1,ntapers)),[1 3 2]),nfft); - P1 = P1(1:h,:,:); - P2 = fft(taparray.*permute(XSEG2(:,:,ones(1,ntapers)),[1 3 2]),nfft); - P2 = P2(1:h,:,:); - - % now make cross-products of them to fill cross-spectrum matrix - - for ii = 1:n1 - o = ind_1(ii); - for jj = 1:n2 - oo = ind_2(jj); - S(:,ii,jj) = S(:,ii,jj) + mean(P1(:,:,o) .* conj(P2(:,:,oo)),2); - end - end - - if ~isempty(ind_pow) - for pp = 1:n3 - u = ind_pow(pp); - pow(:,pp) = pow(:,pp) + mean(P1(:,:,u) .* conj(P1(:,:,u)),2); - end - end - -end - -S = S/nchunks; - -if ~isempty(ind_pow) - pow = pow/nchunks; - - for pp = 1:n3 - clear a b - a = find(ind_1 == ind_pow(pp)); - b = find(ind_2 == ind_pow(pp)); - S(:,a,b) = pow(:,pp); - end -end - diff --git a/libs/pellegrini/fp_cpsd_mt_matlab2015.m b/libs/pellegrini/fp_cpsd_mt_matlab2015.m deleted file mode 100644 index 3973147..0000000 --- a/libs/pellegrini/fp_cpsd_mt_matlab2015.m +++ /dev/null @@ -1,61 +0,0 @@ -function S = fp_cpsd_mt(X1,X2,ind_1,ind_2,h,window,noverlap,nchunks,taparray) - -% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe - -ind_pow = intersect(ind_1, ind_2); -nfft = 2*(h-1); - -n1 = numel(ind_1); -n2 = numel(ind_2); -n3 = numel(ind_pow); -S = complex(zeros(h,n1,n2)); -pow = zeros(h,n3); - -winstep = window-noverlap; -ntapers = size(taparray,2); - -% compute tapered periodogram with FFT - -for k = 1:nchunks - - XSEG1 = X1((1:window) + (k-1)*winstep,:); - XSEG2 = X2((1:window) + (k-1)*winstep,:); - - % compute periodogram - P1 = fft(taparray.*permute(XSEG1(:,:,ones(1,ntapers)),[1 3 2]),nfft); - P1 = P1(1:h,:,:); - P2 = fft(taparray.*permute(XSEG2(:,:,ones(1,ntapers)),[1 3 2]),nfft); - P2 = P2(1:h,:,:); - - % now make cross-products of them to fill cross-spectrum matrix - - for ii = 1:n1 - o = ind_1(ii); - for jj = 1:n2 - oo = ind_2(jj); - S(:,ii,jj) = S(:,ii,jj) + mean(P1(:,:,o) .* conj(P2(:,:,oo)),2); - end - end - - if ~isempty(ind_pow) - for pp = 1:n3 - u = ind_pow(pp); - pow(:,pp) = pow(:,pp) + mean(P1(:,:,u) .* conj(P1(:,:,u)),2); - end - end - -end - -S = S/nchunks; - -if ~isempty(ind_pow) - pow = pow/nchunks; - - for pp = 1:n3 - clear a b - a = find(ind_1 == ind_pow(pp)); - b = find(ind_2 == ind_pow(pp)); - S(:,a,b) = pow(:,pp); - end -end - diff --git a/libs/pellegrini/fp_cpsd_mt_matlab2019.m b/libs/pellegrini/fp_cpsd_mt_matlab2019.m deleted file mode 100644 index 9552a73..0000000 --- a/libs/pellegrini/fp_cpsd_mt_matlab2019.m +++ /dev/null @@ -1,58 +0,0 @@ -function S = fp_cpsd_mt_matlab2019(X1,X2,ind_1,ind_2,h,window,noverlap,nchunks,taparray) - -% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe - -ind_pow = intersect(ind_1, ind_2); -nfft = 2*(h-1); - -n1 = numel(ind_1); -n2 = numel(ind_2); -n3 = numel(ind_pow); -S = complex(zeros(h,n1,n2)); -pow = zeros(h,n3); - -winstep = window-noverlap; -ntapers = size(taparray,2); - -% compute tapered periodogram with FFT - -for k = 1:nchunks - - XSEG1 = X1((1:window) + (k-1)*winstep,:); - XSEG2 = X2((1:window) + (k-1)*winstep,:); - - % compute periodogram - P1 = fft(taparray.*permute(XSEG1(:,:,ones(1,ntapers)),[1 3 2]),nfft); - P1 = P1(1:h,:,:); - P2 = fft(taparray.*permute(XSEG2(:,:,ones(1,ntapers)),[1 3 2]),nfft); - P2 = P2(1:h,:,:); - - % now make cross-products of them to fill cross-spectrum matrix - - for ii = 1:n1 - o = ind_1(ii); - S(:,ii,:) = S(:,ii,:) + mean(repmat(P1(:,:,o), 1, 1, length(ind_2)) .* conj(P2(:,:,ind_2)),2); - end - - if ~isempty(ind_pow) - for pp = 1:n3 - u = ind_pow(pp); - pow(:,pp) = pow(:,pp) + mean(P1(:,:,u) .* conj(P1(:,:,u)),2); - end - end - -end - -S = S/nchunks; - -if ~isempty(ind_pow) - pow = pow/nchunks; - - for pp = 1:n3 - clear a b - a = find(ind_1 == ind_pow(pp)); - b = find(ind_2 == ind_pow(pp)); - S(:,a,b) = pow(:,pp); - end -end - diff --git a/libs/pellegrini/fp_cpsd_welch.m b/libs/pellegrini/fp_cpsd_welch.m deleted file mode 100644 index a3b8096..0000000 --- a/libs/pellegrini/fp_cpsd_welch.m +++ /dev/null @@ -1,26 +0,0 @@ -function S = fp_cpsd_welch(X_1, X_2,ind_1,ind_2,h,window,noverlap) -%X_1 and X_2 can be the same - then the usual CS is -%calculated. Otherwise X_2 should contain data from another trial than -%X_1. ind_1 and ind_2 contain channels of interest. The output S will be of -%size length(ind_1) x length(ind_2) x nfreq. - -% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe - -ind_pow = intersect(ind_1, ind_2); -nfft = 2*(h-1); - -n1 = numel(ind_1); -n2 = numel(ind_2); -S = complex(zeros(n1,n2,h)); - -for ii = 1:n1 - o = ind_1(ii); - S(ii,:,:) = transpose(cpsd(X_1(:,o),X_2(:,ind_2),window,noverlap,nfft)); % cross-spectra - if ismember(o,ind_pow) - clear b - b = find(ind_2 == o); - S(ii,b,:) = transpose(cpsd(X_1(:,o),X_1(:,o),window,noverlap,nfft)); - end - -end - \ No newline at end of file diff --git a/libs/pellegrini/fp_data2bs_event_uni.m b/libs/pellegrini/fp_data2bs_event_uni.m new file mode 100644 index 0000000..fc86d06 --- /dev/null +++ b/libs/pellegrini/fp_data2bs_event_uni.m @@ -0,0 +1,94 @@ +function [cs,nave]=fp_data2bs_event_uni(data,segleng,segshift,epleng,freqpairs,nshuf) +% Calculates bispectral-tensors from data for event-related measurement +% and their null distributions using a shuffling approach. +% +% Based on the data2bs_event function of the METH toolbox by Guido Nolte. +% Modified by (c) 2023 Franziska Pellegrini and Stefan Haufe +% +% usage: [cs,nave]=data2bs_event(data,segleng,segshift,epleng,freqpairs,para); +% +% input: +% data: ndat times nchan matrix each colum is the time-series in one +% channel; +% segleng: length of each segment in bins, e.g. segleng=1000; +% segshift: numer of bins by which neighboring segments are shifted; +% 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 +% para: structure which is eventually used later +% +% output: +% cs: nchan by nchan by nchan by peak_combination (=2) by number of shuffles +% tensor such that +% cs(i,j,k,ishuf)= +% where f1=freqpairs(f,1) and f2=freqpairs(f,2), +% x=fft(data) and the average is over epeochs and segments. +% +% NOTE that this function is restricted to one frequency combination only! +% The addition for frequency BANDS will be added in the future. +% +% if para.fave=0 then cs contains a fifth argument denoting +% the segment. +% if para.fave=1 or ommited, then cs was averaged over segments. + +% nave: number of averages + +[ndat,nchan]=size(data); + +nep=floor(ndat/epleng); + +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); + +%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); + + 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))); + + %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 + + %shape cs: chan x chan x chan x peak_combination x shuffles + cs(:,:,:,1,ishuf) = cs1./nave; + cs(:,:,:,2,ishuf) = cs2./nave; + +end \ No newline at end of file diff --git a/libs/pellegrini/fp_fc_shuffletest.m b/libs/pellegrini/fp_fc_shuffletest.m deleted file mode 100644 index 072efcb..0000000 --- a/libs/pellegrini/fp_fc_shuffletest.m +++ /dev/null @@ -1,23 +0,0 @@ -function fp_fc_shuffletest(isb) -%calculates true MIM and generates MIM null distribution -%isb = subject index -% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe - -DIRIN = './Data_MI/'; - -%subjects with high performance classification -subs = [3 4 5 8 9 11 12 14 15 16 17 18 19 21 22 23 25 27 28 29 30 31 33 34 35 37]; -nshuf = 1001; %first shuffle = true MIM, then 1000 shuffles - -%% load preprocessed EEG -sub = ['vp' num2str(subs(isb))]; -load([DIRIN sub '.mat']) - -%% shuffle -%one shuffle ~ 1 min on local machine -%first shuffle is true MIM -npcs = repmat(3,1,nroi); -MIM_s = fp_shuffle_MIM(data,npcs,fs,filt1, nshuf); - -%% save -save([DIRIN 'RDE_shuf_' num2str(isb) '.mat'],'-v7.3') diff --git a/libs/pellegrini/fp_fft_coeffs.m b/libs/pellegrini/fp_fft_coeffs.m new file mode 100644 index 0000000..5e24efb --- /dev/null +++ b/libs/pellegrini/fp_fft_coeffs.m @@ -0,0 +1,26 @@ +function coeffs = fp_fft_coeffs(data,segleng,segshift,epleng,freqpairs) +% Get Fourier coefficients for epoched data. +% The output coeffs is peaks x nchan x ntrials. +% The 4 peaks correspond to 1) the low-frequeny peak, 2) the left side +% lobe, 3) the high-frequency peak, 4) the right side lobe. +% +% Copyright (c) 2023 Franziska Pellegrini and Stefan Haufe + +[ndat,nchan]=size(data); + +mywindow=repmat(hanning(segleng),1,nchan); +nep=floor(ndat/epleng); + +nseg=floor((epleng-segleng)/segshift)+1; %total number of segments +assert(nseg==1,'only possible with 1 segment') + +for j=1:nep + %disp(j) + dataloc=data((j-1)*epleng+1:j*epleng,:); + datalocfft(:,:,j)=fft(detrend(dataloc).*mywindow); +end + +coeffs(1,:,:) = datalocfft(freqpairs(1),:,:); %f1 +coeffs(2,:,:) = datalocfft(freqpairs(2)-freqpairs(1)+1,:,:); %f2-f1 +coeffs(3,:,:) = datalocfft(freqpairs(2),:,:); %f2 +coeffs(4,:,:) = datalocfft(freqpairs(2)+freqpairs(1)-1,:,:); %f1+f2 \ No newline at end of file diff --git a/libs/pellegrini/fp_npcs2inds.m b/libs/pellegrini/fp_npcs2inds.m deleted file mode 100644 index 65c96ec..0000000 --- a/libs/pellegrini/fp_npcs2inds.m +++ /dev/null @@ -1,19 +0,0 @@ -function [inds, PCA_inds] = fp_npcs2inds(npcs) - -% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe - -%% -beg_inds = cumsum([1 npcs(1:end-1)]); -end_inds = cumsum([npcs]); - -for iroi = 1:numel(npcs) - PCA_inds{iroi} = beg_inds(iroi):end_inds(iroi); -end - -inds = {}; ninds = 0; -for iroi = 1:numel(npcs) - for jroi = (iroi+1):numel(npcs) - inds{ninds+1} = {PCA_inds{iroi}, PCA_inds{jroi}}; - ninds = ninds + 1; - end -end \ No newline at end of file diff --git a/libs/pellegrini/fp_plot_fc_shuffletest.m b/libs/pellegrini/fp_plot_fc_shuffletest.m deleted file mode 100644 index 6e1fbbb..0000000 --- a/libs/pellegrini/fp_plot_fc_shuffletest.m +++ /dev/null @@ -1,36 +0,0 @@ -function fp_plot_fc_shuffletest -%Function that generates p-values from the MIM null distribution and plots -%them -% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe - -DIRIN = '~/Dropbox/Franziska/Data_MEG_Project/RDE_shuffletest/right_MI/'; -DIRFIG = '~/Desktop/'; - -nsub = 26; - -%% generate p-values by comparing true MIM to null distribution -for isb = 1:nsub - in = load([DIRIN 'RDE_shuf_' num2str(isb) '.mat']); - - %average over one region dimension to obtain netMIM - MIM_s = squeeze(mean(in.MIM_s,2)); - MIM_pn(:,isb) = sum(MIM_s(:,1)< MIM_s(:,2:end),2)./(size(in.MIM_s,3)-1); -end - -%% Use Stouffer's method to aggregate p-values -nroi = size(MIM_pn,1); -for iroi = 1:nroi - MIM_pn_s(iroi) = fp_stouffer(squeeze(MIM_pn(iroi,:))); -end - -%% Use FDR-correction for multiple comparison's correction -[p_fdr, ~] = fdr(MIM_pn_s, 0.05); -MIM_pn_s(MIM_pn_s>p_fdr)=1; - -%% Plot -load cm17; -load('bs_results.mat'); % load cortex -MIM_pn_s(MIM_pn_s==0) = 0.001; % 1/nshuf -data = -log10(MIM_pn_s); -allplots_cortex_BS(cortex_highres,data, [35 52], cm17a ,'-log(p)', 0.3,... - [DIRFIG 'netMIM_right']) \ No newline at end of file diff --git a/libs/pellegrini/fp_tsdata_to_cpsd.m b/libs/pellegrini/fp_tsdata_to_cpsd.m deleted file mode 100644 index 8ef4dc4..0000000 --- a/libs/pellegrini/fp_tsdata_to_cpsd.m +++ /dev/null @@ -1,87 +0,0 @@ -function S = fp_tsdata_to_cpsd(X,fres,method,ind_1, ind_2, id_trials_1, id_trials_2, window,noverlap,nw,ntapers) -% Estimate cross-power spectral density from time series data between -% channels ind_1 and channels ind_2. -% Shuffle id_trials_2 for surrogate images. - -% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe - -%% -[n_chan,n_times,n_trials] = size(X); -X = demean(X); -X = permute(X,[2 1 3]); % transpose row, col (works for single-trial data too) -nfft = 2*fres; - -if ~exist('window','var') - window = min(n_times,nfft); % according to Chronux ... by default Matlab 'pwelch' splits data into 8 overlapping segments -end -assert(window <= n_times,'window cannot be longer than data'); - -if ~exist('overlap','var') - noverlap = round(window/2); -end -assert(noverlap < window,'overlap must be shorter than window'); - -if nargin < 3 || isempty(method) - method = 'WELCH'; -end - -if ~isequal(sort(ind_1),sort(unique(ind_1))) || ~isequal(sort(ind_2), sort(unique(ind_2))) - error('ind_1 and ind_2 must be unique') -end - -if strcmpi(method,'MT') - - nchunks = floor(((n_times-noverlap)/(window-noverlap))); % FFT chunks per channel - - if nargin < 10 || isempty(nw) - nw = 3; - end - - if nargin < 11 || isempty(ntapers) - ntapers = 2*nw -1; - end - - tapers = dpss(window,nw,ntapers,'calc'); % Slepian sequences: tapers is a matrix of size window x ntapers - taparray = tapers(:,:,ones(1,n_chan)); - - S = 0; - for r = 1:n_trials % works for single-trial too - - %trial shuffling - x_original = X(:,:,id_trials_1(r)); - x_perm = X(:,:,id_trials_2(r)); - - clear s - s = fp_cpsd_mt_matlab2019(x_original,x_perm,ind_1, ind_2,fres+1,window,noverlap,nchunks,taparray); - S = S+s; - - end - S = permute(S,[2 3 1])/n_trials; - - -elseif strcmpi(method,'WELCH') - - S = 0; - - for r = 1:n_trials - %trial shuffling - x_original = X(:,:,id_trials_1(r)); - x_perm = X(:,:,id_trials_2(r)); - - clear s - s = fp_cpsd_welch(x_original,x_perm,ind_1,ind_2,fres+1,window,noverlap); - S = S + s; - end - S = S/n_trials; - -else - error('unknown method ''%s''',method); -end - - - - - - - - diff --git a/libs/pellegrini/fp_unwrap_conn.m b/libs/pellegrini/fp_unwrap_conn.m deleted file mode 100644 index 007ad05..0000000 --- a/libs/pellegrini/fp_unwrap_conn.m +++ /dev/null @@ -1,50 +0,0 @@ -function [MIM_, MIC_, GC_, DIFFGC_, iCOH_, aCOH_] = fp_unwrap_conn(conn,nroi,filt1,PCA_inds) - -% Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe - -%% -%initialize output variables -MIM_ = []; -MIC_ = []; -GC_=[]; -DIFFGC_ = []; -iCOH_ =[]; -aCOH_ = []; - -% filt1 -% filt1.band_inds - -iinds = 0; -for iroi = 1:nroi - for jroi = (iroi+1):nroi - iinds = iinds + 1; - - if isfield(conn,'MIM') - MIM_(iroi, jroi) = mean(conn.MIM(filt1.band_inds, iinds)); - MIM_(jroi,iroi) = MIM_(iroi,jroi); - end - if isfield(conn,'MIC') - MIC_(iroi, jroi) = mean(conn.MIC(filt1.band_inds, iinds)); - MIC_(jroi,iroi) = MIC_(iroi,jroi); - end - if isfield(conn,'TRGC') - DIFFGC_(iroi,jroi) = mean(squeeze(conn.TRGC(filt1.band_inds,iinds,1) - conn.TRGC(filt1.band_inds,iinds,2))); - DIFFGC_(jroi,iroi) = -DIFFGC_(iroi,jroi); - end - if isfield(conn,'GC') - GC_(iroi,jroi) = mean(squeeze(conn.GC(filt1.band_inds,iinds,1) - conn.GC(filt1.band_inds,iinds,2))); - GC_(jroi,iroi) = -GC_(iroi,jroi); - end - end -end - -if isfield(conn,'COH') - for iroi = 1:nroi - for jroi = 1:nroi - iCOH_(iroi, jroi) = mean(mean(mean(abs(imag(conn.COH(filt1.band_inds, PCA_inds{iroi}, PCA_inds{jroi}))),1), 2), 3); - iCOH_(jroi,iroi) = iCOH_(iroi,jroi); - aCOH_(iroi, jroi) = mean(mean(mean(abs(conn.COH(filt1.band_inds, PCA_inds{iroi}, PCA_inds{jroi})),1), 2), 3); - aCOH_(jroi,iroi) = aCOH_(iroi,jroi); - end - end -end \ No newline at end of file diff --git a/pop_roi_statsplot.m b/pop_roi_statsplot.m index ae3c221..3942d4b 100644 --- a/pop_roi_statsplot.m +++ b/pop_roi_statsplot.m @@ -1,4 +1,5 @@ -% pop_roi_statsplot() - Generate p-values from FC null distributions and plots them. +% pop_roi_statsplot() - Generate p-values from FC null distributions and plots them. Based on Franziska Pellegrini's script +% fp_plot_fc_shuffletest. % % Inputs: % EEG - EEGLAB dataset with ROI activity computed. @@ -18,6 +19,8 @@ % 'MIC' : Maximized Imaginary Coherency for each ROI % 'freqrange' - [min max] frequency range or [integer] single frequency in Hz. Default is to plot broadband power. % 'alpha' - [integer] Significance level. Default is 0.05. +% +% Author: Tien Dung Nguyen, tien-dung.nguyen@charite.de function EEG = pop_roi_statsplot(EEG, varargin) diff --git a/roi_connstats.m b/roi_connstats.m index bf2b051..d8767f1 100644 --- a/roi_connstats.m +++ b/roi_connstats.m @@ -63,7 +63,7 @@ if ~isempty(tmpMethods1) npcs = repmat(EEG.roi.nPCA, 1, EEG.roi.nROI); - conn = shuffle_MIM(data, npcs, tmpMethods1, g.nshuf, 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'poolsize', g.poolsize); % (nfreq, nROI, nROI, nshuf) + conn = shuffle_CS(data, npcs, tmpMethods1, g.nshuf, 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'poolsize', g.poolsize); % (nfreq, nROI, nROI, nshuf) for iMethod = 1:length(tmpMethods1) EEG.roi.(tmpMethods1{iMethod}) = conn.(tmpMethods1{iMethod}); if strcmpi(tmpMethods1{iMethod}, 'MIM') || strcmpi(tmpMethods1{iMethod}, 'MIC') diff --git a/roi_pac.m b/roi_pac.m index 36e849b..b40522f 100644 --- a/roi_pac.m +++ b/roi_pac.m @@ -80,7 +80,7 @@ params.nROI = nROI; EEG.roi.PAC.roi_selection = roi_selection; end - [b_orig, b_anti, b_orig_norm, b_anti_norm] = pac_bispec(data, params); + [b_orig, b_anti, b_orig_norm, b_anti_norm] = data2bs_pac(data, params); % options which bispectral tensors to store switch bs_outopts diff --git a/test_pipes/test_pac.m b/test_pipes/test_pac.m index 3c9ef1c..76dd617 100644 --- a/test_pipes/test_pac.m +++ b/test_pipes/test_pac.m @@ -20,7 +20,7 @@ 'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1); -EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel); +EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3); %% Test bispectrum for single frequency inputs low = 10; @@ -44,6 +44,7 @@ tic EEG4 = pop_roi_connect(EEG, 'methods', {'PAC', 'MIM', 'COH'}, 'fcomb', fcomb); % test all 3 connectivity functions (data2spwctrgc, data2strgcmim, roi_pac)toc toc +EEG5 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'conn_stats', 'off', 'nshuf', 2); %% Test PAC plotting % Test for single frequency inputs diff --git a/utils/pac_bispec.m b/utils/data2bs_pac.m similarity index 81% rename from utils/pac_bispec.m rename to utils/data2bs_pac.m index c70c907..1ee62b1 100644 --- a/utils/pac_bispec.m +++ b/utils/data2bs_pac.m @@ -15,7 +15,7 @@ % Copyright (C) Franziska Pellegrini, franziska.pellegrini@charite.de, % Tien Dung Nguyen, tien-dung.nguyen@charite.de -function [b_orig, b_anti, b_orig_norm,b_anti_norm] = pac_bispec(data, params) +function [b_orig, b_anti, b_orig_norm,b_anti_norm] = data2bs_pac(data, params) % determine ROIs nroi = params.nROI; @@ -57,39 +57,43 @@ for proi = 1:nroi for aroi = proi:nroi - % upper freqs X = data([proi aroi],:,:); + + % upper freqs [bs_up,~] = data2bs_event(X(:,:)', segleng, segshift, epleng, freqinds_up); - biv_orig_up = ([abs(bs_up(1, 2, 2, :)) abs(bs_up(2, 1, 1, :))]); + % call function bs2pac(bs_up), this function does everything below and can hopefully then also be called in shuffle_BS (need to include nshuf info at some point) + + biv_orig_up = squeeze(([mean(abs(bs_up(1, 2, 2, :))) mean(abs(bs_up(2, 1, 1, :)))])); % [Bkmm, Bmkk], average over frequency bands xx = bs_up - permute(bs_up, [2 1 3 4]); %Bkmm - Bmkm - biv_anti_up = ([abs(xx(1, 2, 2, :)) abs(xx(2, 1, 1, :))]); + biv_anti_up = squeeze(([abs(xx(1, 2, 2, :)) abs(xx(2, 1, 1, :))])); % normalized by threenorm [RTP_up,~] = data2bs_threenorm(X(:,:)', segleng, segshift, epleng, freqinds_up); bicoh_up = bs_up ./ RTP_up; - bicoh_up = mean(bicoh_up, 4); + bicoh_up = mean(bicoh_up, 4); % average over frequency bands biv_orig_up_norm = ([abs(bicoh_up(1, 2, 2)) abs(bicoh_up(2, 1, 1))]); xx = bicoh_up-permute(bicoh_up, [2 1 3]); biv_anti_up_norm = ([abs(xx(1, 2, 2)) abs(xx(2, 1, 1))]); % lower freqs [bs_low,~] = data2bs_event(X(:,:)', segleng, segshift, epleng, freqinds_low); - biv_orig_low = ([abs(bs_low(1, 2, 2, :)) abs(bs_low(2, 1, 1, :))]); + biv_orig_low = squeeze(([mean(abs(bs_low(1, 2, 2, :))) mean(abs(bs_low(2, 1, 1, :)))])); xx = bs_low - permute(bs_low, [2 1 3, 4]); - biv_anti_low = ([abs(xx(1, 2, 2, :)) abs(xx(2, 1, 1, :))]); + biv_anti_low = squeeze(([abs(xx(1, 2, 2, :)) abs(xx(2, 1, 1, :))])); % normalized by threenorm [RTP_low,~] = data2bs_threenorm(X(:,:)', segleng, segshift, epleng, freqinds_low); bicoh_low = bs_low ./ RTP_low; - bicoh_low = mean(bicoh_low, 4); + bicoh_low = mean(bicoh_low, 4); % average over frequency bands biv_orig_low_norm = ([abs(bicoh_low(1, 2, 2)) abs(bicoh_low(2, 1, 1))]); xx = bicoh_low-permute(bicoh_low, [2 1 3]); biv_anti_low_norm = ([abs(xx(1, 2, 2)) abs(xx(2, 1, 1))]); + % PAC_km(f1, f2) = 0.5 * |Bkmm(f1, f2-f1)| + 0.5 * |Bkmm(f1, f2)| b_orig(aroi,proi) = mean([biv_orig_up(1) biv_orig_low(1)]); b_orig(proi,aroi) = mean([biv_orig_up(2) biv_orig_low(2)]); b_anti(aroi,proi) = mean([biv_anti_up(1) biv_anti_low(1)]); - b_anti(proi,aroi) = mean([biv_anti_up(2) biv_anti_low(2)]); + b_anti(proi,aroi) = mean([biv_anti_up(2) biv_anti_low(2)]); % normalized versions b_orig_norm(aroi,proi) = mean([biv_orig_up_norm(1) biv_orig_low_norm(1)]); diff --git a/utils/shuffle_BS.m b/utils/shuffle_BS.m index 26d3283..499b72c 100644 --- a/utils/shuffle_BS.m +++ b/utils/shuffle_BS.m @@ -134,13 +134,15 @@ else fcomb = fcomb; end + [BS, ~] =fp_data2bs_event_uni(data(:, :)', ndat, floor(ndat/2), ndat, fcomb, nshuf); % nchan, nchan, nchan, npeaks, nshuf + % bs_up = BS(:, :, :, 1, :) + % bs_low = BS(:, :, :, 2, :) - % function [cs,nave]=fp_data2bs_event_uni(data, segleng, segshift, epleng, freqpairs, nshuf) - [CS, ~] =fp_data2bs_event_uni(data(:, :)', ndat, floor(ndat/2), ndat , fcomb, nshuf); - %[CS, ~] = fp_data2bs_event_uni(data(:, :)', ndat, floor(ndat/2), ndat, freqpairs, nshuf); - - nfreqs = size(CS, 3); - + for proi = 1:nROI + for aroi = proi:nROI + % do your stuff + end + end % PAC = zeros(nfreqs, ninds); diff --git a/utils/shuffle_MIM.m b/utils/shuffle_CS.m similarity index 97% rename from utils/shuffle_MIM.m rename to utils/shuffle_CS.m index eb143c2..40a6262 100644 --- a/utils/shuffle_MIM.m +++ b/utils/shuffle_CS.m @@ -25,7 +25,7 @@ % Stefan Haufe, haufe@tu-berlin.de % Tien Dung Nguyen, tien-dung.nguyen@charite.de -function conn = shuffle_MIM(data, npcs, output, nshuf, varargin) +function conn = shuffle_CS(data, npcs, output, nshuf, varargin) % Copyright (c) 2022 Franziska Pellegrini and Stefan Haufe % decode input parameters @@ -105,7 +105,6 @@ else shuf_inds = randperm(nepo); end - data_shuf = data(:, :, shuf_inds); % Starts here @@ -155,9 +154,6 @@ for iout = 1:length(output) eval(['conn.' output{iout} ' = ' output{iout} '_s;']) end - % shut down current parallel pool - %poolobj = gcp('nocreate'); - %delete(poolobj); % shut down current parallel pool only if the toolbox is available if license('test', 'Distrib_Computing_Toolbox') && ~isempty(ver('parallel'))