diff --git a/README.md b/README.md
index d6c4e75..903ba08 100644
--- a/README.md
+++ b/README.md
@@ -2,10 +2,10 @@
 
 # What is ROIconnect?
 
-ROIconnect is a freely available open-source plugin to [EEGLAB](https://github.com/sccn/eeglab) for EEG data analysis. It allows you to perform functional connectivity analysis between regions of interests (ROIs) on source level.  The results can be visualized in 2-D and 3-D. ROIs are defined based on popular fMRI atlases, and source localization can be performed through LCMV beamforming or eLORETA. Connectivity analysis can be performed between all pairs of brain regions using Granger Causality, Time-reversed Granger Causality, Multivariate Interaction Measure, Maximized Imaginary Coherency, and other methods. This plugin is compatible with Fieldtrip, Brainstorm and NFT head models.
+ROIconnect is a freely available open-source plugin to [EEGLAB](https://github.com/sccn/eeglab) for EEG data analysis. It allows you to perform functional connectivity analysis between regions of interest (ROIs) on source level.  The results can be visualized in 2-D and 3-D. ROIs are defined based on popular fMRI atlases, and source localization can be performed through LCMV beamforming or eLORETA. Connectivity analysis can be performed between all pairs of brain regions using Granger Causality, Time-reversed Granger Causality, Multivariate Interaction Measure, Maximized Imaginary Coherency, and other methods. This plugin is compatible with Fieldtrip, Brainstorm and NFT head models.
 
 📚 Check out the following paper to learn about recommended methods and pipelines for connectivity experiments:
-> Pellegrini, F., Delorme, A., Nikulin, V. & Haufe, S., 2022. Identifying best practices for detecting inter-regional functional connectivity from EEG. bioRxiv 2022.10.05.510753. https://doi.org/10.1101/2022.10.05.510753
+> Pellegrini, F., Delorme, A., Nikulin, V., & Haufe, S. (2023). Identifying good practices for detecting inter-regional linear functional connectivity from EEG. NeuroImage, 120218. [doi: 10.1016/j.neuroimage.2023.120218](https://doi.org/10.1016/j.neuroimage.2023.120218)
 
 ⚠️ Disclaimer: This plugin implements the best-practice pipeline that we identified for our studied setting. We believe it can be used in such environments without hesitation. Additional code to reproduce our experiments entirely is provided in a [separate repository](https://github.com/fpellegrini/FCsim). In the medium term, we intend to extend this plugin for other use cases, which will be backed up by respective validation studies.
 
diff --git a/libs/pellegrini/shuffle_MIM.m b/libs/pellegrini/shuffle_MIM.m
index 132b6ce..6a2463b 100644
--- a/libs/pellegrini/shuffle_MIM.m
+++ b/libs/pellegrini/shuffle_MIM.m
@@ -8,6 +8,7 @@
     ninds = length(inds);
     
     warning('One iteration takes about 90 seconds.')
+    fprintf('Generating null distribution using %d shuffles...\n', nshuf)
     fprintf('Progress of %d:', nshuf);
     for ishuf = 1:nshuf %one iteration takes ~90 sec on my local laptop
         if mod(ishuf, 10) == 0
@@ -52,6 +53,7 @@
         MIM_s(:, :, :, ishuf) = get_connect_mat(conn.MIM, nroi, +1);
         % [MIM_s(:,:,ishuf), ~, ~, ~, ~, ~] = fp_unwrap_conn(conn,nroi,filt1,PCA_inds);
     end
+    fprintf('\n');
 end
 
 function measure = get_connect_mat( measureOri, nROI, signVal)
diff --git a/plotconnectivity.m b/plotconnectivity.m
index 1533e72..924d4a0 100644
--- a/plotconnectivity.m
+++ b/plotconnectivity.m
@@ -193,6 +193,7 @@
 
 % rename labels
 % -------------
+if ~strcmpi(g.brainimg, 'off') % CTC added .... to branch around code if not plotting
 if isempty(g.labels)
     for iPnt = 1:length(anglesInit)
         g.labels{iPnt} = sprintf('  Area %d', iPnt);
@@ -205,7 +206,7 @@
         g.labels{iPnt} = [ ' ' g.labels{iPnt} ];
     end
 end
-
+end % CTC added .... to branch around above code if not plotting
 warning off;
 if isempty(limits)
     limits(2) = max(array(:));
diff --git a/pop_roi_activity.m b/pop_roi_activity.m
index 8255c37..ae38514 100644
--- a/pop_roi_activity.m
+++ b/pop_roi_activity.m
@@ -17,7 +17,12 @@
 %  'resample'         - [integer] resample to the desired sampling rate. Default
 %                       is 100. Adujst the model order accordingly. ROIconnect
 %                       has only be tested with 100 Hz sampling rates.
-%  'fooof'            - ['on'|'off'] enable FOOOF analysis (this method can be used to parameterize neural power spectra and is described here: https://fooof-tools.github.io/fooof/). Default is 'off'.
+%  'effectchanges'    - ['on'|'off'] apply resampling and epoching to the
+%                       dataset given as input. Otherwise, only update
+%                       EEG.roi structure. Default is 'off'.
+%  'fooof'            - ['on'|'off'] enable FOOOF analysis (this method can be 
+%                       used to parameterize neural power spectra and is described here: 
+%                       https://fooof-tools.github.io/fooof/). Default is 'off'.
 %  'fooof_frange'     - [ ] FOOOF fitting range. Default is [1 30] like in the MATLAB example: 
 %                       https://github.com/fooof-tools/fooof_mat/blob/main/examples/fooof_example_one_spectrum.m.
 %  'freqresolution'   - [integer] Desired frequency resolution (in number of frequencies). If
@@ -244,6 +249,7 @@
     'atlas'           'string'              {}               '';
     'resample'        'real'                {}               100;
     'regepochs'       'string'              { 'on' 'off'}    'off'; % ignored
+    'effectchanges'   'string'              { 'on' 'off'}    'off';
     'nPCA'            'real'                {}               3;
     'epochlen'        'real'                {}               2;
     'fooof'           'string'              { 'on' 'off'}    'off';
@@ -251,12 +257,13 @@
     'fooof_frange'     ''                   {}               [1 30]}, 'pop_roi_activity', 'ignore');
 if ischar(g), error(g); end
 
+EEGOUT = EEG;
 if ~isempty(g.resample) && ~isequal(EEG.srate, g.resample)
-    EEG = pop_resample(EEG, g.resample);
+    EEGOUT = pop_resample(EEGOUT, g.resample);
 end
-if EEG.trials == 1
+if EEGOUT.trials == 1
     recurrence = 2;
-    EEG = eeg_regepochs(EEG, recurrence, [0 2]);
+    EEGOUT = eeg_regepochs(EEGOUT, recurrence, [0 2]);
 end
 
 chansel = {};
@@ -276,11 +283,17 @@
     g.leadfield = EEG.dipfit.sourcemodel;
 end    
 
-EEG = roi_activity(EEG, 'leadfield', g.leadfield, 'headmodel', EEG.dipfit.hdmfile, ...
+EEGOUT = roi_activity(EEGOUT, 'leadfield', g.leadfield, 'headmodel', EEG.dipfit.hdmfile, ...
     'model', g.model, 'modelparams', g.modelparams, 'sourcemodel', sourceModelFile, ...
     'sourcemodel2mni', sourceModel2MNI, 'nPCA', g.nPCA,'fooof', g.fooof, 'fooof_frange', g.fooof_frange, ...
     'freqresolution', g.freqresolution, 'sourcemodelatlas', g.atlas, 'chansel', chansel, moreargs{:});
 
+if strcmpi(g.effectchanges, 'on')
+    EEG = EEGOUT;
+else
+    EEG.roi = EEGOUT.roi;
+end
+
 if nargout > 1
     for iOption = 1:2:length(options)
         if strcmpi(options{iOption}, 'leadfield') && isequal(options{iOption+1}, EEG.dipfit.sourcemodel)
diff --git a/pop_roi_connect.m b/pop_roi_connect.m
index 985210c..92dcf31 100644
--- a/pop_roi_connect.m
+++ b/pop_roi_connect.m
@@ -31,7 +31,9 @@
 %                       'MIC'   : Maximized Imaginary Coherency for each ROI
 %                       'PAC'   : Phase-amplitude coupling between ROIs
 %  'snippet'        - ['on'|off]  Option to compute connectivity over snippets. Default is 'off'. 
-%  'snip_length'    - ['on'|'off']  Length of the snippets. Default is 60 seconds.
+%  'firstsnippet'   - ['on'|off]  Only use the first snippet (useful for fast computation). Default is 'off'. 
+%  'snip_length'    - ['on'|'off'] Length of the snippets. Default is 60 seconds.
+%  'erroronsnippet' - ['on'|'off'] Error if snippet too short. Default 'on'.
 %  'fcsave_format'  - ['mean_snips'|'all_snips']  Option to save mean over snippets 
 %                     (shape: 101,68,68) or all snippets (shape: n_snips,101,68,68). Default is 'mean_snips.'
 %  'freqresolution' - [integer] Desired frequency resolution (in number of frequencies). 
@@ -170,7 +172,9 @@
       'naccu'          'integer'  { }                           0;
       'methods'        'cell'     { }                           { };
       'snippet'        'string'   { 'on', 'off' }               'off';
-      'nepochs'        'real'                {}                 [ ];
+      'firstsnippet'   'string'   { 'on', 'off' }               'off'; 
+      'erroronsnippet' 'string'   { 'on', 'off' }               'off'; 
+      'nepochs'        'real'                {}               [];
       'snip_length'    'integer'  { }                           60; 
       'fcsave_format'  'string'   { 'mean_snips', 'all_snips'}  'mean_snips';
       'freqresolution' 'integer'  { }                           0; 
@@ -237,7 +241,12 @@
     snip_eps = snippet_length/(size(EEG.data,2)/EEG.srate); % n epochs in snippet
     nsnips = floor(EEG.trials/snip_eps);
     if nsnips < 1
-        error('Snippet length cannot exceed data length.')
+        if strcmpi(opt.erroronsnippet, 'on')
+            error('Snippet length cannot exceed data length.')
+        else
+            frpintf(2, 'Snippet length cannot exceed data length, using the whole data')
+            nsnips = 1;
+        end
     end
     diff = (EEG.trials * EEG.pnts/EEG.srate) - (nsnips * EEG.pnts/EEG.srate * snip_eps);
     if diff ~= 0
@@ -246,7 +255,9 @@
 
     source_roi_data_save = EEG.roi.source_roi_data;
     for isnip = 1:nsnips
-        roi_snip = source_roi_data_save(:,:,(isnip-1)* snip_eps + 1 : (isnip-1)* snip_eps + snip_eps); % cut source data into snippets
+        begSnip = (isnip-1)* snip_eps + 1;
+        endSnip = min((isnip-1)* snip_eps + snip_eps, size(source_roi_data_save,3));
+        roi_snip = source_roi_data_save(:,:, begSnip:endSnip ); % cut source data into snippets
         EEG.roi.source_roi_data = single(roi_snip);
         EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection); % compute connectivity over one snippet
         for fc = 1:n_conn_metrics 
@@ -254,26 +265,42 @@
             fc_matrix = EEG.roi.(fc_name);
             conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure
         end
-    end
-    
-    % compute mean over connectivity of each snippet
-    for fc = 1:n_conn_metrics
-        fc_name = g.methods{fc};
-        [first_dim, second_dim, third_dim] = size(conn_matrices_snips{1,fc});
 
-        conn_cell = conn_matrices_snips(:,fc); % store all matrices of one metric in a cell
-        mat = cell2mat(conn_cell);
-        reshaped = reshape(mat, first_dim, nsnips, second_dim, third_dim);
-        reshaped = squeeze(permute(reshaped, [2,1,3,4]));
-        if strcmpi(g.fcsave_format, 'all_snips')
-            EEG.roi.(fc_name) = reshaped;
-        else
-            if nsnips > 1
-                mean_conn = squeeze(mean(reshaped, 1)); 
+        if strcmpi(g.firstsnippet, 'on')
+            nsnips = 1;
+        end
+    
+        source_roi_data_save = EEG.roi.source_roi_data;
+        for isnip = 1:nsnips
+            roi_snip = source_roi_data_save(:,:,(isnip-1)* snip_eps + 1 : (isnip-1)* snip_eps + snip_eps); % cut source data into snippets
+            EEG.roi.source_roi_data = single(roi_snip);
+            EEG = roi_connect(EEG, 'morder', g.morder, 'naccu', g.naccu, 'methods', g.methods,'freqresolution',g.freqresolution); % compute connectivity over one snippet
+            for fc = 1:n_conn_metrics 
+                fc_name = options{2}{fc};
+                fc_matrix = EEG.roi.(fc_name);
+                conn_matrices_snips{isnip,fc} = fc_matrix; % store each connectivity metric for each snippet in separate structure
+            end
+        end
+        
+        % compute mean over connectivity of each snippet
+        for fc = 1:n_conn_metrics
+            fc_name = options{2}{fc};
+            [first_dim, second_dim, third_dim] = size(conn_matrices_snips{1,fc});
+    
+            conn_cell = conn_matrices_snips(:,fc); % store all matrices of one metric in a cell
+            mat = cell2mat(conn_cell);
+            reshaped = reshape(mat, first_dim, nsnips, second_dim, third_dim);
+            reshaped = squeeze(permute(reshaped, [2,1,3,4]));
+            if strcmpi(g.fcsave_format, 'all_snips')
+                EEG.roi.(fc_name) = reshaped;
             else
-                mean_conn = reshaped;
+                if nsnips > 1
+                    mean_conn = squeeze(mean(reshaped, 1)); 
+                else
+                    mean_conn = reshaped;
+                end
+                EEG.roi.(fc_name) = mean_conn; % store mean connectivity in EEG struct
             end
-            EEG.roi.(fc_name) = mean_conn; % store mean connectivity in EEG struct
         end
     end
 elseif strcmpi(g.conn_stats, 'on')
diff --git a/pop_roi_connectplot.m b/pop_roi_connectplot.m
index 3499481..0b47ac8 100644
--- a/pop_roi_connectplot.m
+++ b/pop_roi_connectplot.m
@@ -763,10 +763,15 @@ function roi_plotcoloredlobes( EEG, matrix, titleStr, measure, hemisphere, group
         imagesc(matrix); colormap(cmap);
     end
     cb = colorbar;
-    tf = isMATLABReleaseOlderThan("R2022a");
-    if tf
+%     tf = isMATLABReleaseOlderThan("R2022a");
+%     if tf
+%         caxis([clim_min clim_max])
+%     else
+%         clim([clim_min clim_max])
+%     end
+    try 
         caxis([clim_min clim_max])
-    else
+    catch
         clim([clim_min clim_max])
     end
     if isDKatlas == 1
diff --git a/pop_roi_statsplot.m b/pop_roi_statsplot.m
new file mode 100644
index 0000000..ef3864e
--- /dev/null
+++ b/pop_roi_statsplot.m
@@ -0,0 +1,87 @@
+% pop_roi_statsplot() - Generate p-values from FC null distributions and plots them.
+%
+% Inputs:
+%   EEG        -  EEGLAB dataset with ROI activity computed.
+%
+% Optional inputs:
+%   'measure'  - [cell] Cell of strings corresponding to methods.
+%                       'CS'    : Cross spectrum
+%                       'aCOH'  : Coherence
+%                       'cCOH'  : (complex-valued) coherency
+%                       'iCOH'  : absolute value of the imaginary part of coherency
+%                       'wPLI'  : Weighted Phase Lag Index
+%                       'PDC'   : Partial directed coherence
+%                       'TRPDC' : Time-reversed partial directed coherence
+%                       'DTF'   : Directed transfer entropy
+%                       'TRDTF' : Time-reversed directed transfer entropy
+%                       'MIM'   : Multivariate Interaction Measure for each ROI
+%                       'MIC'   : Maximized Imaginary Coherency for each ROI
+%  'freqrange' - [min max] frequency range in Hz. Default is to plot broadband power.
+
+function EEG = pop_roi_statsplot(EEG, varargin)
+
+    if nargin < 2
+        help roi_connstats;
+        return
+    end
+
+    if ~isfield(EEG, 'roi') || ~isfield(EEG.roi, 'source_roi_data')
+        error('Cannot find ROI data - compute ROI data first');
+    end
+
+    % decode input parameters
+    % -----------------------
+    g = finputcheck(varargin,  { 
+        'measure'        'string'   { }    '';
+        'freqrange'      'real'     { }    []}, 'pop_roi_statsplot');
+    if ischar(g), error(g); end
+    S = EEG.roi;
+
+    if isempty(g.measure)
+        error('You must define a measure to plot');
+    end
+    
+    % extract frequency indices
+    if ~isempty(g.freqrange)
+        try
+            frq_inds = find(S.freqs >= g.freqrange(1) & S.freqs < g.freqrange(2));
+            title = sprintf('%1.1f-%1.1f Hz frequency band', g.freqrange(1), g.freqrange(2));
+        catch
+            frq_inds = find(S.freqs == g.freqrange(1)); % if a single frequency is passed
+            title = sprintf('%1.1f Hz', g.freqrange(1));
+        end
+    else
+        frq_inds = 1:length(S.freqs);
+        title = 'broadband';
+    end
+    
+    % select frequency or frequency band
+    if length(frq_inds) > 1
+        matrix = squeeze(mean(S.(g.measure)(frq_inds, :, :, :)));
+    else
+        matrix = squeeze(S.(g.measure)(frq_inds, :, :, :));
+    end
+    
+    % generate p-values by comparing the true FC (first shuffle) to null distribution
+    netFC = squeeze(mean(matrix, 2));
+    FC_pn = sum(netFC(:, 1) < netFC(:, 2:end), 2)./(size(matrix, 3) - 1);
+
+    % use FDR-correction for multiple comparison's correction
+    alpha = 0.05; % add this as a parameter?
+    [p_fdr, ~] = fdr(FC_pn, alpha);
+    FC_pn(FC_pn > p_fdr) = 1;
+
+    % plot 
+    load cm17;
+    load cortex; 
+    FC_pn(FC_pn==0) = 1 / (size(netFC, 2) - 1);  % 1 / nshuf
+    data = -log10(FC_pn);
+    try
+        allplots_cortex_BS(cortex_highres, data, [min(data) max(data)], cm17a ,'-log(p)', 0.3);
+        h = textsc(title, 'title');
+        set(h, 'fontsize', 20);
+    catch
+        warning('There are no "significant" p-values to be plotted. These are the p-values:')
+        disp(FC_pn)
+    end
+end
\ No newline at end of file
diff --git a/roi_activity.m b/roi_activity.m
index c43b0e2..5c277e2 100644
--- a/roi_activity.m
+++ b/roi_activity.m
@@ -220,7 +220,11 @@
     error('There must be the same number of vertices/voxels in the leadfield and source model');
 end
 if isempty(g.chansel)
-    g.chansel = [1:EEG.nbchan];
+    if isfield(EEG.dipfit, 'chansel')
+        g.chansel = EEG.dipfit.chansel;
+    else
+        g.chansel = 1:EEG.nbchan;
+    end
 else
     g.chansel = eeg_decodechan(EEG.chanlocs, g.chansel);
 end
@@ -328,11 +332,12 @@
     
     % compute power using the Welch method
     disp('Computing ROI activity...');
-    [tmpWelch,ftmp] = pwelch(tmpData(:,:), data_pnts, floor(data_pnts/2), data_pnts, EEG.srate); % ftmp should be equal frqs 
+    [tmpWelch,ftmp] = pwelch(tmpData, data_pnts, floor(data_pnts/2), data_pnts, EEG.srate); % ftmp should be equal frqs 
     tmpWelch = reshape(tmpWelch, size(tmpWelch,1), EEG.trials, size(source_voxel_data,2), size(source_voxel_data,3));
     tmpWelch = squeeze(mean(tmpWelch,2)); % remove trials size freqs x voxels x 3
     tmpWelch = squeeze(mean(tmpWelch,3)); % remove 3rd dim size freqs x voxels
-    
+    source_voxel_power = tmpWelch;
+
     % fooof settings
     if strcmpi(g.fooof, 'on')
         f_range = g.fooof_frange; % freq range where 1/f should be fitted 
@@ -377,12 +382,15 @@
 if strcmpi(g.exportvoxact, 'on')
     EEG.roi.source_voxel_data     = source_voxel_data; % large (takes lots of RAM)
 end
+EEG.roi.source_voxel_power    = single(source_voxel_power);
 EEG.roi.source_roi_data       = single(source_roi_data);
 EEG.roi.source_roi_power      = source_roi_power; % used for plotting
 EEG.roi.source_roi_power_norm = source_roi_power_norm; % used for cross-sprectum
 EEG.roi.freqs     = frqs;
 EEG.roi.nPCA      = g.nPCA;
 EEG.roi.nROI      = nROI;
+EEG.roi.pnts      = EEG.pnts;
+EEG.roi.srate     = EEG.srate;
 EEG.roi.atlas     = cortex.Atlas;
 EEG.roi.srate     = EEG.srate;
 EEG.roi.leadfield = g.leadfield;
diff --git a/roi_connect.m b/roi_connect.m
index 2d6bfda..e554d6c 100644
--- a/roi_connect.m
+++ b/roi_connect.m
@@ -114,9 +114,9 @@
     tmpMethods1 = intersect(g.methods, methodset1);
     if ~isempty(tmpMethods1)
         if isempty(g.roi_selection)
-            conn_mult = data2sctrgcmim(source_roi_data, EEG.pnts, g.morder, 0, g.naccu, [], inds, tmpMethods1, [], 'freqresolution', g.freqresolution);
+            conn_mult = data2sctrgcmim(source_roi_data, EEG.roi.pnts, g.morder, 0, g.naccu, [], inds, tmpMethods1, [], 'freqresolution', g.freqresolution);
         elseif ~isempty(g.roi_selection)
-            conn_mult = data2sctrgcmim(source_roi_data, EEG.pnts, g.morder, 0, g.naccu, [], inds, tmpMethods1, [], 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'nPCA', nPCA);
+            conn_mult = data2sctrgcmim(source_roi_data, EEG.roi.pnts, g.morder, 0, g.naccu, [], inds, tmpMethods1, [], 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'nPCA', nPCA);
             EEG.roi.roi_selection = g.roi_selection;
         end
         fields = fieldnames(conn_mult);
@@ -137,7 +137,7 @@
             elseif  strcmpi(tmpMethods1{iMethods}, 'GC') || strcmpi(tmpMethods1{iMethods}, 'TRGC')
                 TRGCnet = EEG.roi.(tmpMethods1{iMethods})(:, :, 1) - EEG.roi.(tmpMethods1{iMethods})(:, :, 2);
                 EEG.roi.(tmpMethods1{iMethods}) = get_connect_mat( TRGCnet, nROI, -1); 
-            else
+            else % wPLI
                 measure = rm_components(EEG.roi.(tmpMethods1{iMethods}), EEG.roi.nPCA, tmpMethods1{iMethods}); % only keep the first principal component
                 EEG.roi.(tmpMethods1{iMethods}) = measure;
             end
@@ -147,9 +147,9 @@
     tmpMethods2 = intersect(g.methods, methodset2);
     if ~isempty(tmpMethods2)
         if isempty(g.roi_selection)
-            conn_mult = data2spwctrgc(source_roi_data, EEG.pnts, g.morder, 0, g.naccu, [], tmpMethods2, [], 'freqresolution', g.freqresolution);
+            conn_mult = data2spwctrgc(source_roi_data, EEG.roi.pnts, g.morder, 0, g.naccu, [], tmpMethods2, [], 'freqresolution', g.freqresolution);
         else
-            conn_mult = data2spwctrgc(source_roi_data, EEG.pnts, g.morder, 0, g.naccu, [], tmpMethods2, [], 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'nPCA', nPCA);
+            conn_mult = data2spwctrgc(source_roi_data, EEG.roi.pnts, g.morder, 0, g.naccu, [], tmpMethods2, [], 'freqresolution', g.freqresolution, 'roi_selection', g.roi_selection, 'nPCA', nPCA);
         end
         fields = fieldnames(conn_mult);
         for iField = 1:length(fields)
diff --git a/roi_connstats.m b/roi_connstats.m
index 436fb15..b98ac47 100644
--- a/roi_connstats.m
+++ b/roi_connstats.m
@@ -10,8 +10,6 @@
 %                       'aCOH'  : Coherence
 %                       'cCOH'  : (complex-valued) coherency
 %                       'iCOH'  : absolute value of the imaginary part of coherency
-%                       'GC'    : Granger Causality
-%                       'TRGC'  : Time-reversed Granger Causality
 %                       'wPLI'  : Weighted Phase Lag Index
 %                       'PDC'   : Partial directed coherence
 %                       'TRPDC' : Time-reversed partial directed coherence
@@ -19,7 +17,7 @@
 %                       'TRDTF' : Time-reversed directed transfer entropy
 %                       'MIM'   : Multivariate Interaction Measure for each ROI
 %                       'MIC'   : Maximized Imaginary Coherency for each ROI
-%  nshuf'        - [integer] number of shuffles for statistical significance testing. The first shuffle is the true value. Default is 101. 
+%  nshuf         - [integer] number of shuffles for statistical significance testing. The first shuffle is the true value. Default is 1001. 
 %  roi_selection - [cell array of integers] Cell array of ROI indices {1, 2, 3, ...} indicating for which regions (ROIs) connectivity should be computed. 
 %                  Default is all (set to EEG.roi.nROI).
 
diff --git a/roi_lookupregions.m b/roi_lookupregions.m
index 56d5ab6..d46ac09 100644
--- a/roi_lookupregions.m
+++ b/roi_lookupregions.m
@@ -75,15 +75,18 @@
 end
 
 indNet = [];
+len    = 0;
 for iNet = 1:length(allnetworknames)
     if contains(networkname, allnetworknames{iNet})
-        indNet = iNet;
-        break
+        if length(  allnetworknames{iNet} ) > len
+            indNet = iNet;
+            len    = length( allnetworknames{iNet} );
+        end
     end
 end
 
 if isempty(indNet)
-    fprintf('Network %s not found\n', networkname);
+    fprintf('The brain region "%s" was not found (it is ok if it is not a brain region)\n', networkname);
 elseif ~isstruct(networkdefs)
     regions = roiTable.(allnetworknames{iNet});
 else
diff --git a/roi_pac.m b/roi_pac.m
index 372f67a..36e849b 100644
--- a/roi_pac.m
+++ b/roi_pac.m
@@ -52,7 +52,7 @@
     segleng = ndat;
     segshift = floor(ndat/2);
     epleng = ndat;
-    fs = EEG.srate;
+    fs = EEG.roi.srate;
 
     params = [];
     params.nROI = nROI;
diff --git a/roi_sourceplot.m b/roi_sourceplot.m
index e33c23c..5a0bc67 100644
--- a/roi_sourceplot.m
+++ b/roi_sourceplot.m
@@ -12,9 +12,13 @@
 %                MNI locations.
 %
 % Required inputs:
-%  'freqselect' - [real] frequency of interest or frequency range of interest.
+%  'freqrange' - [real] frequency of interest or frequency range of interest.
 %                 Defaut is all frequencies.
 %
+% Example:
+%  % run ROI activity using ROI connect to compute activity in ROI (and voxels)
+%  roi_sourceplot(EEG.roi.freqs, EEG.roi.source_voxel_power, EEG.roi.cortex );
+%
 % Author:  Arnaud Delorme
 
 % Copyright (C) Arnaud Delorme, arnodelorme@gmail.com
@@ -51,13 +55,17 @@
 g = finputcheck(varargin, { ...
     'freqrange'       {'cell' 'real'}  { [] [] }       {};
     'saveasfile'      'string'         { }             '';
+    'mask'            'string'         { }             '';
     'noplot'          'string'         { 'on' 'off' }  'off';
     'precomputed'     'struct'         { }             struct([]);
     'limits'          ''               []              [];
-    'slice'           'integer'        []              [5 8 10 15 18]}, 'roi_sourceplot');
+    'slice'           'integer'        []              []}, 'roi_sourceplot');
 if ischar(g)
     error(g);
 end
+if isempty(g.freqrange)
+    g.freqrange = [freqs(1) freqs(end)];
+end
 if ~iscell(g.freqrange)
     g.freqrange = { g.freqrange };
 end
@@ -68,23 +76,7 @@
 alldata = [];
 for iFreq = 1:length(g.freqrange)
     
-    g.freqselect = g.freqrange{iFreq};
-    if isempty(g.freqselect)
-        indFreq = 1:length(freqs);
-    elseif length(g.freqselect) == 1
-        [~,indFreq] = min(abs(freqs-g.freqselect));
-        indFreq(2) = indFreq;
-    elseif length(g.freqselect) == 2
-        [~,indFreq1] = min(abs(freqs-g.freqselect(1)));
-        [~,indFreq2] = min(abs(freqs-g.freqselect(2)));
-        indFreq = indFreq1:indFreq2;
-    else
-        error('Frequency selection must be an array of 1 or 2 elements');
-    end
-    if isempty(indFreq)
-        error('No frequency found');
-    end
-    
+   
     % transform to volume
     if ischar(sourcemodel)
         sourceProjtmp = load('-mat', sourcemodel);
@@ -94,62 +86,96 @@
     if ~isfield(sourceProjtmp, 'Vertices')
         sourceProjtmp.Vertices = sourceProjtmp.pos;
     end
-    xl = [min(sourceProjtmp.Vertices(:,1)) max(sourceProjtmp.Vertices(:,1)) ];
-    yl = [min(sourceProjtmp.Vertices(:,2)) max(sourceProjtmp.Vertices(:,2)) ];
-    zl = [min(sourceProjtmp.Vertices(:,3)) max(sourceProjtmp.Vertices(:,3)) ];
-    downscale = 5;
-    volMat = zeros(ceil(diff(xl)/downscale+1), ceil(diff(yl)/downscale+1), ceil(diff(zl)/downscale+1));
-    clear sourcemodelout;
-    sourcemodelout.dim = size(volMat);
-    [r,c,v] = ind2sub(size(volMat),find(volMat == 0));
-    sourcemodelout.pos = [r,c,v];
-    sourcemodelout.inside    = zeros(size(sourcemodelout.pos,1),1);
-    sourcemodelout.unit      = 'mm';
-    %sourcemodelout.transform = traditionaldipfit([xl(1)-5 yl(1)-5 zl(1)-5 0 0 0 5 5 5]);
-    sourcemodelout.transform = [5 0 0 -75;0 5 0 -105;0 0 5 -45;0 0 0 1]; % 5mm grid
-    allInds = zeros(size(sourceProjtmp.Vertices));
-    allIndVolume = zeros(length(sourceProjtmp.Vertices),1);
-    for iVert = 1:length(sourceProjtmp.Vertices)
-        xVal = round((sourceProjtmp.Vertices(iVert,1)-xl(1))/downscale)+1;
-        yVal = round((sourceProjtmp.Vertices(iVert,2)-yl(1))/downscale)+1;
-        zVal = round((sourceProjtmp.Vertices(iVert,3)-zl(1))/downscale)+1;
-        ind = sub2ind(size(volMat), xVal, yVal, zVal);
-        volMat(xVal, yVal, zVal) = mean(sourceact(iVert,indFreq), 2);
-        allIndVolume(iVert) = ind;
-        allInds(iVert,:) = [xVal yVal zVal];
-        sourcemodelout.inside(ind) = true;
+
+    stride = unique(abs(diff(sourceProjtmp.Vertices)));
+    if iscell(sourceact)
+        % select precomputed activity matrix
+        volMat = loreta2volume(sourceact{iFreq}, sourceProjtmp.Vertices, [stride(2) stride(2) stride(2)]);
+    else
+        % select frequency range in activity matrix
+        g.freqselect = g.freqrange{iFreq};
+        if isempty(g.freqselect)
+            indFreq = 1:length(freqs);
+        elseif length(g.freqselect) == 1
+            [~,indFreq] = min(abs(freqs-g.freqselect));
+            indFreq(2) = indFreq;
+        elseif length(g.freqselect) == 2
+            [~,indFreq1] = min(abs(freqs-g.freqselect(1)));
+            [~,indFreq2] = min(abs(freqs-g.freqselect(2)));
+            indFreq = indFreq1:indFreq2;
+        else
+            error('Frequency selection must be an array of 1 or 2 elements');
+        end
+        if isempty(indFreq)
+            error('No frequency found');
+        end
+        volMat = loreta2volume(mean(sourceact(:,indFreq), 2), sourceProjtmp.Vertices, [stride(2) stride(2) stride(2)]);
     end
+    if stride(2) ~= floor(stride(2))
+        error('Non-integer stride')
+    end
+
+%     xl = [min(sourceProjtmp.Vertices(:,1)) max(sourceProjtmp.Vertices(:,1)) ];
+%     yl = [min(sourceProjtmp.Vertices(:,2)) max(sourceProjtmp.Vertices(:,2)) ];
+%     zl = [min(sourceProjtmp.Vertices(:,3)) max(sourceProjtmp.Vertices(:,3)) ];
+%     downscale = 10;
+%     volMat = zeros(ceil(diff(xl)/downscale+1), ceil(diff(yl)/downscale+1), ceil(diff(zl)/downscale+1));
+%     clear sourcemodelout;
+%     sourcemodelout.dim = size(volMat);
+%     [r,c,v] = ind2sub(size(volMat),find(volMat == 0));
+%     sourcemodelout.pos = [r,c,v];
+%     sourcemodelout.inside    = zeros(size(sourcemodelout.pos,1),1);
+%     sourcemodelout.unit      = 'mm';
+%     %sourcemodelout.transform = traditionaldipfit([xl(1)-5 yl(1)-5 zl(1)-5 0 0 0 5 5 5]);
+%     sourcemodelout.transform = [5 0 0 -75;0 5 0 -105;0 0 5 -45;0 0 0 1]; % 5mm grid
+%     allInds = zeros(size(sourceProjtmp.Vertices));
+%     allIndVolume = zeros(length(sourceProjtmp.Vertices),1);
+%     for iVert = 1:length(sourceProjtmp.Vertices)
+%         xVal = round((sourceProjtmp.Vertices(iVert,1)-xl(1))/downscale)+1;
+%         yVal = round((sourceProjtmp.Vertices(iVert,2)-yl(1))/downscale)+1;
+%         zVal = round((sourceProjtmp.Vertices(iVert,3)-zl(1))/downscale)+1;
+%         ind = sub2ind(size(volMat), xVal, yVal, zVal);
+%         volMat(xVal, yVal, zVal) = mean(sourceact(iVert,indFreq), 2);
+%         allIndVolume(iVert) = ind;
+%         allInds(iVert,:) = [xVal yVal zVal];
+%         sourcemodelout.inside(ind) = true;
+%     end
     
     % put precomputed data in VolMat
-    for iSlice = 1:length(g.slice)
-        fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), g.slice(iSlice));
+    if ~isempty(g.slice)
+        tmpSlice = g.slice;
+    else
+        tmpSlice = ceil(linspace(3, size(volMat,3)-2, 5));
+    end
+    for iSlice = 1:length(tmpSlice)
+        fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), tmpSlice(iSlice));
         if isfield(g.precomputed, fieldVal)
             res = g.precomputed.(fieldVal);
             res(isnan(res)) = 0;
             res(:,:) = res(end:-1:1,:);
             res = res';
-            volMat(:,:,g.slice(iSlice)) = res;
+            volMat(:,:,tmpSlice(iSlice)) = res;
         end
     end
     
     % plot
-    res = squeeze(volMat(:,:,g.slice));
+    res = squeeze(volMat(:,:,tmpSlice));
     if isempty(g.limits)
-        mi = min(res(:));
-        mx = max(res(:));
+        mi = min(res(:),[],'omitnan');
+        mx = max(res(:),[],'omitnan');
         fprintf('Loreta limits from data: %1.2f to %1.2f\n', mi, mx);
     else
         mi = 0;
         mx = g.limits(iFreq);
     end
-    cmap = colormap('jet');
-    for iSlice = 1:length(g.slice)
-        res = squeeze(volMat(:,:,g.slice(iSlice)));
+    cmap = colormap('turbo');
+    for iSlice = 1:length(tmpSlice)
+        res = squeeze(volMat(:,:,tmpSlice(iSlice)));
         res = res';
         res(:,:) = res(end:-1:1,:);
         
         % save and retreive data
-        fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), g.slice(iSlice));
+        fieldVal = sprintf('loreta%1.0fto%1.0fHz_slice%d', g.freqrange{iFreq}(1), g.freqrange{iFreq}(2), tmpSlice(iSlice));
         alldata.(fieldVal) = res;
         
         if strcmpi(g.noplot, 'off')
@@ -157,15 +183,19 @@
             for iPix1 = 1:size(res,1)
                 for iPix2 = 1:size(res,2)
                     if res(iPix1,iPix2) ~= 0
-                        ind = ceil((res(iPix1,iPix2)-mi)/(mx-mi)*(size(cmap,1)-1))+1;
-                        ind = max(1, ind);
-                        ind = min(size(cmap,1), ind);
-                        resrgb(iPix1,iPix2,:) = cmap(ind,:);
+                        if isnan(res(iPix1,iPix2))
+                            resrgb(iPix1,iPix2,:) = [0.9 0.9 0.9];
+                        else
+                            ind = ceil((res(iPix1,iPix2)-mi)/(mx-mi)*(size(cmap,1)-1))+1;
+                            ind = max(1, ind);
+                            ind = min(size(cmap,1), ind);
+                            resrgb(iPix1,iPix2,:) = cmap(ind,:);
+                        end
                     end
                 end
             end
             
-            subplot(length(g.slice), length(g.freqrange), iFreq + length(g.freqrange)*(iSlice-1));
+            subplot(length(tmpSlice), length(g.freqrange), iFreq + length(g.freqrange)*(iSlice-1));
             imagesc(resrgb); axis equal; axis off;
             if iSlice == 1
                 if length(g.freqrange{iFreq}) == 2
@@ -216,4 +246,24 @@
 cfg.location     = [26 8 10];
 cfg.latency      = 10;
 cfg.slicepos     = [8 9];
-ft_sourceplot(cfg, sourceProj);
\ No newline at end of file
+ft_sourceplot(cfg, sourceProj);
+
+% convert coordinates to volume
+% -----------------------------
+function vol = loreta2volume(act, pos, stride)
+
+minx = min(pos(:,1)); maxx = max(pos(:,1));
+miny = min(pos(:,2)); maxy = max(pos(:,2));
+minz = min(pos(:,3)); maxz = max(pos(:,3));
+
+vol = zeros((maxx-minx)/stride(1), (maxy-miny)/stride(2), (maxz-minz)/stride(3));
+
+for iPos = 1:size(pos,1)
+    vol( (pos(iPos,1)-minx)/stride(1)+1, ...
+         (pos(iPos,2)-miny)/stride(2)+1, ...
+         (pos(iPos,3)-minz)/stride(3)+1) = act(iPos);
+    % in case using the sourcemodel with shrunk coordinates
+    %     vol( floor((pos(iPos,1)-minx)/stride(1)+1), ...
+    %          floor((pos(iPos,2)-miny)/stride(2)+1), ...
+    %          floor((pos(iPos,3)-minz)/stride(3)+1)) = act(iPos);
+end
diff --git a/test_pipes/pipeline_connectivity.m b/test_pipes/pipeline_connectivity.m
index 6e9ea48..ddce455 100644
--- a/test_pipes/pipeline_connectivity.m
+++ b/test_pipes/pipeline_connectivity.m
@@ -25,7 +25,7 @@
 EEG = pop_dipfit_settings( EEG, 'hdmfile',fullfile(eeglabp, 'plugins','dipfit','standard_BEM','standard_vol.mat'), ...
     'coordformat','MNI','mrifile',fullfile(eeglabp, 'plugins','dipfit','standard_BEM','standard_mri.mat'),...
     'chanfile',fullfile(eeglabp, 'plugins','dipfit','standard_BEM','elec', 'standard_1005.elc'),...
-    'coord_transform',[0.83215 -15.6287 2.4114 0.081214 0.00093739 -1.5732 1.1742 1.0601 1.1485] ,'chansel',[1:32] );
+    'coord_transform',[0.83215 -15.6287 2.4114 0.081214 0.00093739 -1.5732 1.1742 1.0601 1.1485] );
 
 EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'plugins','dipfit','LORETA-Talairach-BAs.mat'), ...
     'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1);
diff --git a/test_pipes/test_mim_stats.m b/test_pipes/test_mim_stats.m
index a4e037e..4660e30 100644
--- a/test_pipes/test_mim_stats.m
+++ b/test_pipes/test_mim_stats.m
@@ -8,7 +8,7 @@
 EEG = pop_resample( EEG, 100);
 EEG = pop_epoch( EEG, { }, [-0.5 1.5], 'newname', 'EEG Data epochs epochs', 'epochinfo', 'yes');
 EEG = pop_select( EEG, 'trial',1:30);
-[ALLEEG, EEG, CURRENTSET] = eeg_store(ALLEEG, EEG);
+% [ALLEEG, EEG, CURRENTSET] = eeg_store(ALLEEG, EEG);
 eeglab redraw;
 
 EEG = pop_dipfit_settings( EEG, 'hdmfile',fullfile(eeglabp, 'plugins','dipfit','standard_BEM','standard_vol.mat'), ...
@@ -22,21 +22,6 @@
 EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3);
 
 %% Create null distribution
-EEG = pop_roi_connect(EEG, 'methods', {'COH', 'MIM'}, 'conn_stats', 'on', 'nshuf', 3);
-
-% ROI selection
-MIM = squeeze(mean(EEG.roi.MIM, 1)); % broadband for now
-
-% generate p-values by comparing true MIM to null distribution
-netMIM = squeeze(mean(MIM, 2));
-MIM_pn = sum(netMIM(:,1) < netMIM(:,2:end),2)./(size(MIM,3)-1);
-
-% plot 
-load cm17;
-load cortex; 
-MIM_pn(MIM_pn==0) = 1 / (size(netMIM, 2)-1);  % 1 / nshuf
-data = -log10(MIM_pn);
-% allplots_cortex_BS(cortex_highres, data, [min(data) max(data)], cm17a ,'-log(p)', 0.3);
-allplots_cortex_BS(cortex_highres, data, [0 5], cm17a ,'-log(p)', 0.3);
-
-
+% EEG = pop_roi_connect(EEG, 'methods', {'COH', 'MIM'}, 'conn_stats', 'on', 'nshuf', 3);
+load('test_pipes/MIM_shuf.mat')
+pop_roi_statsplot(EEG, 'measure', 'MIM', 'freqrange', [8 13]);
\ No newline at end of file