Merge pull request #48 from nguyen-td/feature

Feature

libs/pellegrini/fp_cpsd_mt.m

@@ -0,0 +1,61 @@
+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
+

libs/pellegrini/fp_cpsd_mt_matlab2015.m

@@ -0,0 +1,61 @@
+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
+

libs/pellegrini/fp_cpsd_mt_matlab2019.m

@@ -0,0 +1,58 @@
+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
+

libs/pellegrini/fp_cpsd_welch.m

@@ -0,0 +1,26 @@
+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       
+

libs/pellegrini/fp_fc_shuffletest.m

@@ -0,0 +1,23 @@
+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')

libs/pellegrini/fp_npcs2inds.m

@@ -0,0 +1,19 @@
+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

libs/pellegrini/fp_plot_fc_shuffletest.m

@@ -0,0 +1,36 @@
+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'])

libs/pellegrini/fp_tsdata_to_cpsd.m

@@ -0,0 +1,87 @@
+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
+
+
+
+
+
+
+
+