diff --git a/.gitignore b/.gitignore deleted file mode 100644 index 3e3461a..0000000 --- a/.gitignore +++ /dev/null @@ -1 +0,0 @@ -*.asv diff --git a/eegplugin_roiconnect.m b/eegplugin_roiconnect.m index feb6f7c..ef393bc 100644 --- a/eegplugin_roiconnect.m +++ b/eegplugin_roiconnect.m @@ -48,6 +48,7 @@ addpath(fullfile(p, 'libs/Daniele_ARMA')); addpath(fullfile(p, 'libs/export_fig')); addpath(fullfile(p, 'libs/haufe')); + addpath(fullfile(p, 'libs/jurhar')); addpath(fullfile(p, 'libs/mvgc_v1.0')); addpath(fullfile(p, 'libs/mvgc_v1.0/core')); addpath(fullfile(p, 'libs/mvgc_v1.0/stats')); diff --git a/libs/jurhar/bispectral_TD_est.m b/libs/jurhar/bispectral_TD_est.m new file mode 100644 index 0000000..5a37140 --- /dev/null +++ b/libs/jurhar/bispectral_TD_est.m @@ -0,0 +1,51 @@ +% Implements Bispectral Time Delay Estimations from Nikias Paper +% (1988). Inputs are precomputed univariate bispectra. +% +% Inputs: +% B_xxx, B_xyx, B_yyy - (n_freq x n_freq) univariate bispectra +% method - determines the TDE method, must be between 1:4. +% frqs_idx - (n_freq x 1) vector of frequency bins? +% zeropad - zero-padding, used for frequency-domain bispectral estimates +% +% Outputs: +% T - (1 x 2 * n_freq - 1) time-domain time delay estimate (bispectral hologram) +% I - (n_freq x 2 * n_freq - 1) method-dependent coefficient to estimate bispectral delays + +function [T,I] = bispectral_TD_est(B_xxx, B_xyx, B_yyy, method, frqs_idx, zeropad) + + if nargin < 6 + zeropad = 1; + end + + assert(ismember(method,1:4), 'Method must be [1;4].'); + [m,n,boots] = size(B_xxx); + + switch method + case 1 + phi = angle(B_xyx)-angle(B_xxx); + I = exp(1i.*phi); + case 2 + phi = angle(B_xyx)- 0.5 * (angle(B_xxx) + angle(B_yyy)); + I = exp(1i.*phi); + case 3 + I = B_xyx ./ B_xxx; + case 4 + phi = angle(B_xyx)- 0.5 * (angle(B_xxx) + angle(B_yyy)); + I = abs(B_xyx).* exp(1i.*phi) ./ sqrt(abs(B_xxx).*abs(B_yyy)); + end + + % translate back to time domain. + if nargin > 4 & ~isempty(frqs_idx) + frqs_idx = [frqs_idx zeros(1,n-length(frqs_idx))]; + I = frqs_idx.' * frqs_idx .* I; + end + + % use for frequency-domain bispec estimates (not time domain) + if zeropad + I = [I zeros(m,n-1,boots)]; + end + + T = sum(I,1); + T = abs(fftshift(ifft(T),2)); + +end \ No newline at end of file diff --git a/libs/jurhar/combine_sn.m b/libs/jurhar/combine_sn.m new file mode 100644 index 0000000..f0b4a87 --- /dev/null +++ b/libs/jurhar/combine_sn.m @@ -0,0 +1,31 @@ +function [X,Y] = combine_sn(rawX,noiseX,noiseY,delay,samplefreq,snr,theta,beta) +% Combines signals and noise vectors with specified SNR (alpha), theta, +% beta. + + [rawX,rawY] = mk_delay(delay,rawX); + + % zero-center signal and noise + rawX = rawX-mean(rawX); + rawY = rawY-mean(rawY); + noiseX = noiseX-mean(noiseX); + noiseY = noiseY-mean(noiseY); + + [~,~,~,pX] = obw(rawX,samplefreq); + [~,~,~,pY] = obw(rawY,samplefreq); + [~,~,~,pnoiseX] = obw(noiseX,samplefreq); + [~,~,~,pnoiseY] = obw(noiseY,samplefreq); + + % Normalize signal and noise powers + rawX = rawX / sqrt(pX); + rawY = rawY /sqrt(pY); + noiseX = noiseX / sqrt(pnoiseX); + noiseY = noiseY /sqrt(pnoiseY); + + if nargin < 8 | isempty(beta) + beta =1; + end + + % Combine signal and noise with specified SNR + X = snr * rawX + (1-snr) * (noiseX + theta(1) * noiseY); + Y = snr * beta * rawY + (1-snr) * (noiseY + theta(2) * noiseX); +end \ No newline at end of file diff --git a/libs/jurhar/mk_delay.m b/libs/jurhar/mk_delay.m new file mode 100644 index 0000000..fd54004 --- /dev/null +++ b/libs/jurhar/mk_delay.m @@ -0,0 +1,31 @@ +function [X,Y] = mk_delay(delay,sig1,sig2) +%% Syntax +% [X,Y] = mk_delay(delay,sig1,sig2) +%% Arguments +% _input_ +% delay integer, introduced delay +% sig1 (1,N1) time series 1 +% sig2 (1,N2) time series 2 (optional) +% _output_ +% X (1,N1-delay), sig1 derived time series +% Y (1,N2-delay), sig2 derived time series +%% Description +% Takes time series sig1 and sig2 and time shifts them by delay bins. delay +% can take positive and negative values. If only sig1 is provided, then X,Y +% are time-shifted instances of sig1. +% +% (C) Tin Jurhar, 2022 + + if nargin < 3 + sig2 = sig1; + end + + if delay < 0 + X = sig1(1:end-abs(delay)); + Y = sig2(max(abs(delay))+1:end); + else + X = sig1(max(abs(delay))+1:end); + Y = sig2(1:end-abs(delay)); + end +end + diff --git a/libs/jurhar/mk_series.m b/libs/jurhar/mk_series.m new file mode 100644 index 0000000..0197c35 --- /dev/null +++ b/libs/jurhar/mk_series.m @@ -0,0 +1,38 @@ +% Function to generate simulated signals. Signals are generated either as +% as samples from an exponential distribution, a Gaussian distribution or +% from a Gaussian pink noise process. The option to model the signal as +% spikes is disabled here. +% +% Inputs: +% t - [integer] length of signal in datapoints +% source_type - [String] underlying model of the signal: 'exp', 'G', 'pink', 'pinksq' +% param - [integer] noise parameter +% noise_mask - [boolean] if 1, additional Gaussian noise is added to the signal +% +% Output: +% signal - (1 x t) simulated signal of length t + +function [signal] = mk_series(t,source_type,param,noise_mask) + switch source_type + case 'exp' + signal= exprnd(param,1,t); + case 'G' % Gaussian Noise + signal = randn(1,t); + case 'pink' + signal = mkpinknoise(t,1,param).'; + case 'pinksq' + signal = mkpinknoise(t,1,param).'; + signal = signal.^2; +% case 'spiky' +% load('spike_signal.mat'); +% signal = spikes_signal(randperm(t)).'; + end + + if nargin >3 & noise_mask % add Gaussian noise for numerical stability + signal = signal./norm(signal,'fro'); + noise= randn(1,t); + noise = noise./ norm(noise,'fro'); + + signal = 0.95 * signal + 0.05 * noise; + end +end diff --git a/libs/jurhar/segment_data.m b/libs/jurhar/segment_data.m new file mode 100644 index 0000000..705dcd7 --- /dev/null +++ b/libs/jurhar/segment_data.m @@ -0,0 +1,35 @@ +function Xs = segment_data(X,seglen,shift) + +%% Syntax +% Xs = segment_data(X,seglen,shift) +%% Arguments +% _input_ +% X (1,N) Time series of length N +% seglen integer, desired segment length +% shift integer, number of bins shifted for next segment,(optional) +% _output_ +% Xs (nseg,seglen) output matrix of segmeneted timeseries +%% Description +% Segments input time series Xs (row vector) into segments of length seglen +% with desired shift. +% +% (C) Tin Jurhar, 2022 + if nargin < 3 | isempty(shift) + shift = seglen; + end + + assert(shift<=seglen,"shift must be geq seglen.") + assert(length(size(X))==2 & (size(X,1)==1 | size(X,2)==1),"X must be 1-dimensional.") + if size(X,2)==1 + X = X.'; + end + assert(seglen <= length(X),'seglen must be geq than the input vector.') + + [~,lenX] = size(X); + Xs = []; + i = 1; + while i+seglen-1 <= lenX + Xs = [Xs; X(i:i+seglen-1)]; + i = i+shift; + end +end \ No newline at end of file diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 2adc07b..328457c 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -29,7 +29,8 @@ % 'TRDTF' : Time-reversed directed transfer entropy % 'MIM' : Multivariate Interaction Measure for each ROI % 'MIC' : Maximized Imaginary Coherency for each ROI -% 'PAC' : Phase-amplitude coupling between ROIs +% 'PAC' : Phase-amplitude coupling between ROIs (see the 'fcomb' and 'bs_outopts' input parameters) +% 'TDE' : Time-delay estimation between two selected ROIs (see the 'tde_regions' and 'tde_freqbands' input parameters) % 'snippet' - ['on'|off] Option to compute connectivity over snippets. Default is 'off'. % '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. @@ -50,6 +51,11 @@ % 5 - only store: b_anti, b_anti_norm % '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 empty (in this case, connectivity will be computed for all ROIs). +% 'tde_method' - [integer] TDE method, must be between 1:4, open bispectral_TD_est.m for details. Default is 1. +% 'tde_regions' - [seed target] Array containing the seed and target region for time-delay estimation. Regions need to be passed as region *indices*, +% e.g., [2 10] will compute time-delays from region 2 -> 10 and 10 -> 2, corresponding to the information flow in both directions separately. +% Default is [] (will throw an error). +% 'tde_freqbands' - [f1 f2] Array containing the frequency band for which bispectral time-delays will be estimated. Default is [] (broadband). % 'conn_stats' - ['on'|'off'] Run statistics on connectivity metrics. Default is 'off'. % 'nshuf' - [integer] number of shuffles for statistical significance testing. The first shuffle is the true value. Default is 1001. % 'freqrange' - [min max] frequency range in Hz. This is used to compute and plot p-values. Default is to plot broadband power. @@ -181,6 +187,9 @@ 'fcomb' 'struct' { } struct; 'bs_outopts' 'integer' { } 1; 'roi_selection' 'cell' { } { }; + 'tde_method' 'integer' { 1:4 } 1; + 'tde_regions' 'integer' { } [ ]; + 'tde_freqbands' 'integer' { } [ ]; 'conn_stats' 'string' { } 'off'; ... 'nshuf' 'integer' { } 1001; ... 'poolsize' 'integer' { } 1}, 'pop_roi_connect'); @@ -197,6 +206,11 @@ return end +tmpMethods = setdiff(g.methods, { 'CS' 'COH' 'cCOH' 'aCOH' 'iCOH' 'GC' 'TRGC' 'wPLI' 'PDC' 'TRPDC' 'DTF' 'TRDTF' 'MIM' 'MIC' 'PAC' 'TDE'}); +if ~isempty(tmpMethods) + error('Unknown methods %s', vararg2str(tmpMethods)) +end + % compute connectivity over snippets if strcmpi(g.snippet, 'on') && strcmpi(g.conn_stats, 'off') % n_conn_metrics = length(g.methods); @@ -340,6 +354,9 @@ if strcmpi(g.snippet, 'off') && ~isempty(intersect(g.methods, {'PAC'})) EEG = roi_pac(EEG, g.fcomb, g.bs_outopts, g.roi_selection); end + if strcmpi(g.snippet, 'off') && ~isempty(intersect(g.methods, {'TDE'})) + EEG = roi_tde(EEG, g.tde_method, g.tde_regions, g.tde_freqbands); + end end % TO-DO: add snippet option for stats mode diff --git a/pop_roi_connectplot.m b/pop_roi_connectplot.m index 1b720ee..86a2279 100644 --- a/pop_roi_connectplot.m +++ b/pop_roi_connectplot.m @@ -22,8 +22,10 @@ % 'iCOH': Absolute value of the imaginary part of Coherency % 'MIC' : Maximized Imaginary Coherency for each ROI % 'MIM' : Multivariate Interaction Measure for each ROI -% 'pac' : Phase-amplitude coupling for a certain frequency (band) combination based on bicoherence -% 'pac_anti': Phase-amplitude coupling for a certain frequency (band) combination based on the antisymmetrized bicoherence +% 'PAC' : Phase-amplitude coupling for a certain frequency (band) combination based on bicoherence +% 'PAC_anti': Phase-amplitude coupling for a certain frequency (band) combination based on the antisymmetrized bicoherence +% 'TDE': Bispectral time-delay estimation +% 'TDE_anti': Antisymmetrized bispectral time-delay estimation % 'freqrange' - [min max] frequency range or [integer] single frequency in Hz. Default is to plot broadband power. % 'smooth' - [float] smoothing factor for cortex surface plotting % 'plotcortex' - ['on'|'off'] plot results on smooth cortex. Default is 'on' @@ -32,6 +34,7 @@ % 'plot3d' - ['on'|'off'] ... Default is 'off' % 'plot3dparams' - [cell] optional parameters for the generation of brain movies. Check the related documentation in roi_plotbrainmovie.m % 'plotmatrix' - ['on'|'off'] plot results as ROI to ROI matrix. Default is 'off' +% 'plotbutterfly' - ['on'|'off'] plot results as a frequency x connectivity plot (butterfly plot). Default is 'off' % 'plotbarplot' - ['on'|'off'] plot ROI based power spectrum as barplot. Default is 'off' % 'hemisphere' - ['all'|'left'|'right'] hemisphere options for ROI to ROI matrix. Default is 'all' % 'grouphemispheres' - ['on'|'off'] group ROIs by hemispheres (left hemisphere, then right hemisphere). Default is 'off' @@ -155,18 +158,18 @@ splot(end ).unit = 'aCOH'; splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = -1; + splot(end ).psd = 0; splot(end ).plot3d = plot3dFlag; end if isfield(EEG.roi, 'cCOH') splot(end+1).label = 'ROI to ROI coherency'; splot(end ).labelshort = 'Coherency'; - splot(end ).acronym = 'cCoh'; + splot(end ).acronym = 'cCOH'; splot(end ).unit = 'cCOH'; splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = -1; + splot(end ).psd = 0; splot(end ).plot3d = plot3dFlag; end @@ -177,7 +180,7 @@ splot(end ).unit = 'iCOH'; splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = -1; + splot(end ).psd = 0; splot(end ).plot3d = plot3dFlag; end @@ -188,7 +191,7 @@ splot(end ).unit = 'GC'; splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = -1; + splot(end ).psd = 0; splot(end ).plot3d = plot3dFlag; end @@ -199,7 +202,7 @@ splot(end ).unit = 'TRGC'; % not used yet splot(end ).cortex = cortexFlag; splot(end ).matrix = 1; - splot(end ).psd = -1; + splot(end ).psd = 0; splot(end ).plot3d = plot3dFlag; end @@ -246,6 +249,28 @@ splot(end ).psd = 0; splot(end ).plot3d = plot3dFlag; end + + if isfield(EEG.roi, 'TDE') + splot(end+1).label = 'ROI to ROI Time-delay estimation'; + splot(end ).labelshort = 'Time-delay estimation'; + splot(end ).acronym = 'TDE'; % TDE based on bispectrum + splot(end ).unit = 'TDE'; % not used yet + splot(end ).cortex = cortexFlag; + splot(end ).matrix = 1; + splot(end ).psd = -1; + splot(end ).plot3d = plot3dFlag; + end + + if isfield(EEG.roi, 'TDE') + splot(end+1).label = 'ROI to ROI Time-delay estimation'; + splot(end ).labelshort = 'Time-delay estimation'; + splot(end ).acronym = 'TDE_anti'; % TDE based on antisymmetrized bispectrum + splot(end ).unit = 'TDE'; % not used yet + splot(end ).cortex = cortexFlag; + splot(end ).matrix = 1; + splot(end ).psd = -1; + splot(end ).plot3d = plot3dFlag; + end if nargin < 2 @@ -327,6 +352,7 @@ 'plot3d' 'string' { 'on' 'off' } 'off'; 'plot3dparams' 'cell' { } {}; 'plotmatrix' 'string' { 'on' 'off' } 'off'; + 'plotbutterfly' 'string' { 'on' 'off' } 'off'; 'noplot' 'string' { 'on' 'off' } 'off'; 'plotbarplot' 'string' { 'on' 'off'} 'off'; 'hemisphere' 'string' {'all' 'left' 'right'} 'all'; @@ -402,7 +428,7 @@ error('This option is obsolete'); end - if strcmpi(g.plotcortex, 'on') && strcmpi(lower(g.measure), 'roipsd') + if strcmpi(g.plotcortex, 'on') && strcmpi(g.measure, 'roipsd') cortexPlot = 10*log10( mean(EEG.roi.source_roi_power(frq_inds,:)) ); end @@ -421,9 +447,10 @@ % TRGCnet = S.TRGC(:, :, 1) - S.TRGC(:, :, 2); TRGC = S.TRGC; end - matrix = squeeze(mean(TRGC(frq_inds, :, :))); + butterflyplot = squeeze(mean(TRGC, 2)); + matrix = squeeze(mean(TRGC(frq_inds, :, :), 1)); cortexPlot = mean(matrix, 2); - cortexTitle = [ upper(g.measure) ' (' titleStr '); Red = net sender; Blue = net receiver' ]; + cortexTitle = [ g.measure ' (' titleStr '); Red = net sender; Blue = net receiver' ]; case { 'mim' 'mic' } if strcmpi(g.measure, 'MIC') @@ -433,16 +460,19 @@ % MI = S.MIM(:, :); MI = S.MIM; end - matrix = squeeze(mean(MI(frq_inds, :, :),1)); + butterflyplot = squeeze(mean(MI, 2)); + matrix = squeeze(mean(MI(frq_inds, :, :), 1)); cortexPlot = mean(matrix, 2); case { 'acoh' 'ccoh' 'icoh'} if strcmpi(g.measure, 'aCOH') + butterflyplot = squeeze(mean(S.aCOH, 2)); matrix = squeeze(mean(S.aCOH(frq_inds, :, :), 1)); elseif strcmpi(g.measure, 'cCOH') error(['Complex values are not supported. To plot the absolute values, compute "aCOH", ' ... 'to plot the imaginary part, compute "iCOH".']) else + butterflyplot = squeeze(mean(S.iCOH, 2)); matrix = squeeze(mean(S.iCOH(frq_inds, :, :), 1)); end cortexPlot = mean(matrix, 2); @@ -450,14 +480,17 @@ case { 'crossspecpow' 'coh' 'crossspecimag' } if strcmpi(g.measure, 'coh') PS = abs(S.COH); % do not know what to do here + butterflyplot = squeeze(mean(PS, 2)); PSarea2area = squeeze(mean(PS(frq_inds, :, :))); cortexPlot = mean(PSarea2area, 2); elseif strcmpi(g.measure, 'crossspecimag') PS = abs(imag(cs2coh(S.CS))); + butterflyplot = squeeze(mean(PS, 2)); PSarea2area = squeeze(mean(PS(frq_inds, :, :))); cortexPlot = mean(PSarea2area, 2); else PS = cs2psd(S.CS); + butterflyplot = squeeze(mean(PS, 2)); apow = squeeze(sum(sum(reshape(PS(frq_inds, :), [], S.nROI), 1), 2)).*S.source_roi_power_norm'; cortexPlot = 10*log10(apow); PSarea2area = []; @@ -486,6 +519,19 @@ end cortexPlot = mean(matrix, 2); + case {'tde' 'tde_anti'} + if strcmpi(g.measure, 'tde') + t_xy = S.TDE.T_XY; + t_yx = S.TDE.T_YX; + else + t_xy = S.TDE.aT_XY; + t_yx = S.TDE.aT_YX; + end + plot_tde(t_xy, S.TDE.shift, S.TDE.region_X, S.TDE.region_Y, S.TDE.method) + plot_tde(t_yx, S.TDE.shift, S.TDE.region_Y, S.TDE.region_X, S.TDE.method) + + % disable cortex plot + g.plotcortex = 'off'; end % get seed @@ -519,15 +565,17 @@ % plot on cortical surface if strcmpi(g.plotcortex, 'on') && cortexFlag ~= -1 - cortexTitle = [ plotOpt.labelshort ' (' titleStr ')' ]; + if ~strcmpi(g.measure, 'gc') || ~strcmpi(g.measure, 'trgc') + cortexTitle = [ plotOpt.labelshort ' (' titleStr ')' ]; + end if isempty(g.plotcortexseedregion) - allplots_cortex_BS(cortex_highres, cortexPlot, [min(cortexPlot) max(cortexPlot)], cm17a, upper(g.measure), g.smooth); + allplots_cortex_BS(cortex_highres, cortexPlot, [min(cortexPlot) max(cortexPlot)], cm17a, g.measure, g.smooth); % allplots_cortex_BS(cortex_highres, cortexPlot, [min(cortexPlot) max(cortexPlot)], cm17a, upper(splot.unit), g.smooth); % allplots_cortex_BS(S.cortex, cortexPlot, [min(cortexPlot) max(cortexPlot)], cm17a, plotOpt.unit, g.smooth); else cortexTitle = [ cortexTitle ' for area ' int2str(seed_idx)]; cortexPlot = squeeze(matrix(seed_idx,:)); - allplots_cortex_BS(S.cortex, cortexPlot, [min(cortexPlot) max(cortexPlot)], cm17a, upper(g.measure), g.smooth, [], {coordinate}); + allplots_cortex_BS(S.cortex, cortexPlot, [min(cortexPlot) max(cortexPlot)], cm17a, g.measure, g.smooth, [], {coordinate}); % allplots_cortex_BS(cortex_highres, cortexPlot, [min(cortexPlot) max(cortexPlot)], cm17a, upper(g.measure), g.smooth, [], {coordinate}); % allplots_cortex_BS(S.cortex, cortexPlot, [min(cortexPlot) max(cortexPlot)], cm17a, plotOpt.unit, g.smooth, [], {coordinate}); end @@ -537,6 +585,21 @@ warning('EEG.roi.cortex does not contain the field "Faces" required to plot surface topographies.') end + % make butterfly plot + if strcmpi(g.plotbutterfly, 'on') && ~isempty(matrix) + if strcmpi(g.measure, 'pac') || strcmpi(g.measure, 'pac_anti') + warning('Butterfly plots (frequency x connectivity) cannot be computed for PAC because frequencies have already been specified for the computation.') + else + figure; plot(EEG.roi.freqs, butterflyplot, 'LineWidth', 1) + h = title([ 'ROI to ROI ' replace_underscores(g.measure) ' (' titleStr ')' ]); + set(h, 'fontsize', 16); + xlabel('Frequency (Hz)') + ylabel([replace_underscores(g.measure) ' (a.u.)']) + grid on + end + + end + % plot 3D if strcmpi(g.plot3d, 'on') && ~isempty(matrix) PSarea2area = matrix.*seedMask; @@ -572,27 +635,27 @@ end end -function labels = get_labels(EEG) - % retrieve labels from atlas - labels = strings(1,length(EEG.roi.atlas.Scouts)); - for i = 1:length(labels) - scout = struct2cell(EEG.roi.atlas.Scouts(i)); - labels(i) = char(scout(1)); - end - labels = cellstr(labels); - - % remove region labels that were not selected - if isfield(EEG.roi, 'roi_selection') - if ~isempty(EEG.roi.roi_selection) - labels = labels(cell2mat(EEG.roi.roi_selection)); - end - end -end - -function new_labels = replace_underscores(labels) - % remove underscores in label names to avoid bug - new_labels = strrep(labels, '_', ' '); -end +% function labels = get_labels(EEG) +% % retrieve labels from atlas +% labels = strings(1,length(EEG.roi.atlas.Scouts)); +% for i = 1:length(labels) +% scout = struct2cell(EEG.roi.atlas.Scouts(i)); +% labels(i) = char(scout(1)); +% end +% labels = cellstr(labels); +% +% % remove region labels that were not selected +% if isfield(EEG.roi, 'roi_selection') +% if ~isempty(EEG.roi.roi_selection) +% labels = labels(cell2mat(EEG.roi.roi_selection)); +% end +% end +% end +% +% function new_labels = replace_underscores(labels) +% % remove underscores in label names to avoid bug +% new_labels = strrep(labels, '_', ' '); +% end function [colors, color_idxx, roi_idxx, labels_sorted, roi_loc] = get_colored_labels(EEG) labels = get_labels(EEG); @@ -832,7 +895,8 @@ function roi_plotcoloredlobes( EEG, matrix, titleStr, measure, hemisphere, group else set(gca,'ytick',1:n_roi_labels,'yticklabel',labels(hem_idx{1}:hem_idx{2}:n_roi_labels), 'fontsize', 7, 'TickLength',[0.015, 0.02], 'LineWidth',0.75); end - h = title([ 'ROI to ROI ' upper(replace_underscores(measure)) ' (' titleStr ')' ]); + set(gca, 'YDir','normal') + h = title([ 'ROI to ROI ' replace_underscores(measure) ' (' titleStr ')' ]); set(h, 'fontsize', 16); xtickangle(90) end @@ -901,3 +965,17 @@ function roi_largeplot(EEG, mim, trgc, roipsd, titleStr) set(lgd, 'Position', [0.44 0.06 0.25 0.25]); end end + +function plot_tde(T, shift, region_X, region_Y, method) + [~, peak_idx] = max(T); + est_delay = shift(peak_idx); % in Hz + + figure; plot(shift, T, 'LineWidth', 1) + xline(est_delay, '--r') + xlabel('Time (s)') + ylabel('a.u.') + h = title(sprintf('%s -> %s TDE | Method %d', region_X, region_Y, method)); + set(h, 'fontsize', 16); + subtitle("\tau = " + num2str(est_delay) + " (s)") + grid on +end diff --git a/roi_connect.m b/roi_connect.m index 236bc78..8edd3e3 100644 --- a/roi_connect.m +++ b/roi_connect.m @@ -80,10 +80,10 @@ 'roi_selection' 'cell' { } { } }, 'roi_connect'); if ischar(g), error(g); end if isempty(g.naccu), g.naccu = 0; end - tmpMethods = setdiff(g.methods, { 'CS' 'COH' 'cCOH' 'aCOH' 'iCOH' 'GC' 'TRGC' 'wPLI' 'PDC' 'TRPDC' 'DTF' 'TRDTF' 'MIM' 'MIC' 'PAC'}); - if ~isempty(tmpMethods) - error('Unknown methods %s', vararg2str(tmpMethods)) - end +% tmpMethods = setdiff(g.methods, { 'CS' 'COH' 'cCOH' 'aCOH' 'iCOH' 'GC' 'TRGC' 'wPLI' 'PDC' 'TRPDC' 'DTF' 'TRDTF' 'MIM' 'MIC' 'PAC'}); +% if ~isempty(tmpMethods) +% error('Unknown methods %s', vararg2str(tmpMethods)) +% end inds = {}; ninds = 0; if isempty(g.roi_selection) diff --git a/roi_pac.m b/roi_pac.m index 7a333c5..43ad1a9 100644 --- a/roi_pac.m +++ b/roi_pac.m @@ -1,4 +1,4 @@ -% roi_pac() - Compute phase-amplitude coupling (PAC) between ROIs (cf. Zandvoort and Nolte, 2021). +% roi_pac() - Computes phase-amplitude coupling (PAC) between ROIs (cf. Zandvoort and Nolte, 2021). % Wrapper function for pac_bispec(). % % Usage: @@ -7,10 +7,11 @@ % Inputs: % EEG - EEGLAB dataset with ROI activity computed. % fcomb - [struct] Frequency combination for which PAC is computed (in Hz). Must have fields 'low' and -% 'high' with fcomb.low < fcomb.high. For example, fcomb.low = 10 and fcomb.high = 50 if single -% frequencies are used. fcomb.low = [4 8] and fcomb.high = [48 50] if frequency bands are used -% (might take a long time to compute, so use with caution). Default is {} (this will cause an error). -% bs_outopts - [integer] Option which bispectral tensors should be stored in EEG.roi.PAC. Default is 1. +% 'high' with fcomb.low < fcomb.high. For example, fcomb.low = 10 and fcomb.high = 50 if single +% frequencies are used. fcomb.low = [4 8] and fcomb.high = [48 50] if frequency bands are used +% (might take a long time to compute, so use with caution). Default is {} (this will cause an error). +% bs_outopts - [integer] Option which bispectral tensors should be stored in EEG.roi.PAC. Default is 1. "orig" means +% original bispectrum, "anti" means antisymmetrized bispectrum, "norm" means normalized bispectrum. % 1 - store all tensors: b_orig, b_anti, b_orig_norm, b_anti_norm % 2 - only store: b_orig, b_anti % 3 - only store: b_orig_norm, b_anti_norm @@ -38,7 +39,7 @@ if ~isfield(fcomb, 'low') || ~isfield(fcomb, 'high') help roi_pac; - error('Frequency pair cannot be found - check the documentation for the "fcomb" input parameter in "roi_pac".') + error('Frequency pair cannot be found - check the documentation for the "fcomb" input parameter in "roi_pac.m".') end if fcomb.high < fcomb.low diff --git a/roi_tde.m b/roi_tde.m new file mode 100644 index 0000000..4281119 --- /dev/null +++ b/roi_tde.m @@ -0,0 +1,101 @@ +% roi_tde() - Performs time-delay estimation (TDE) based on antisymmetrized bispectra +% between two channels. (cf. Jurhar et al., 2024, Nikias and Pan, 1988). +% +% Usage: +% EEG = roi_pac(EEG, method, tde_regions, freqbands); +% +% Inputs: +% EEG - EEGLAB dataset with ROI activity computed. Will contain a new EEG.roi.TDE field with the following subfields: +% T_XY - bispectral TDE from region X -> Y +% aT_XY - antisyemmtrized bispectral TDE from region X -> Y +% T_YX - bispectral TDE from region Y -> X +% aT_YX - bispectral TDE from region Y -> X +% shift - time shift (x-axis) +% method - same as 'tde_method' +% region_X - region X with its actual label according to EEG.roi.atlas.Scouts +% region_Y - region Y with its actual label according to EEG.roi.atlas.Scouts +% tde_method - [integer] TDE method, must be between 1:4, open bispectral_TD_est.m for details. Default is 1. +% tde_regions - [seed target] Array containing the seed and target region for time-delay estimation. Regions need to be passed as region *indices*, +% e.g., [2 10] will compute time-delays from region 2 -> 10 and 10 -> 2, corresponding to the information flow in both directions separately. +% Default is [] (will throw an error). +% tde_freqbands - [f1 f2] Array containing the frequency band for which bispectral time-delays will be estimated. Default is [] (broadband). + +function EEG = roi_tde(EEG, tde_method, tde_regions, tde_freqbands) + + if nargin < 4 + help roi_tde; + return + end + + if ~isfield(EEG, 'roi') || ~isfield(EEG.roi, 'source_roi_data') + error('Cannot find ROI data - compute ROI data first'); + else + data = EEG.roi.source_roi_data; + end + + if isempty(tde_regions) + help roi_tde + error('No seed and target regions selected - check the the documentation for the "tde_regions" input parameter in "roi_tde.m"') + end + + % set parameters + ndat = EEG.pnts * 2; + segshift = floor(ndat/2); + segleng = ndat; + epleng = ndat; + fres = EEG.srate; + frqs = sfreqs(fres, EEG.srate); + maxfreqbins = floor(segleng/2); + + % only keeep first PC + nPCA = EEG.roi.nPCA; + if nPCA > 1 + warning('Only the first principal component will be used to determine TDE.') + data = data(1:nPCA:end, :, :); + end + + % get seed and target region + X = data(tde_regions(1), :)'; + Y = data(tde_regions(2), :)'; + + % compute time-delays + [T_XY, ~, aT_XY, ~] = compute_TDE(X, Y, tde_method, tde_freqbands, segleng, segshift, epleng, maxfreqbins, frqs); + [T_YX, ~, aT_YX, ~] = compute_TDE(Y, X, tde_method, tde_freqbands, segleng, segshift, epleng, maxfreqbins, frqs); + + % store delays in EEG struct + EEG.roi.TDE.T_XY = T_XY; + EEG.roi.TDE.aT_XY = aT_XY; + EEG.roi.TDE.T_YX = T_YX; + EEG.roi.TDE.aT_YX = aT_YX; + EEG.roi.TDE.shift = (-floor(segleng/2)+1:floor(segleng/2)-1) / EEG.srate; % x-axis + EEG.roi.TDE.method = tde_method; + EEG.roi.TDE.region_X = EEG.roi.atlas.Scouts(tde_regions(1)).Label; + EEG.roi.TDE.region_Y = EEG.roi.atlas.Scouts(tde_regions(2)).Label; +end + +function [T, I, aT, aI] = compute_TDE(X, Y, tde_method, tde_freqbands, segleng, segshift, epleng, maxfreqbins, frqs) + + % compute univariate bispectra + [B2_xxx] = squeeze(data2bs_univar(X, segleng, segshift, epleng, maxfreqbins)); + para_xyx.chancomb = [1, 2, 1]; + [B2_xyx] = data2bs_univar([X, Y], segleng, segshift, epleng, maxfreqbins, para_xyx); + [B2_yyy] = squeeze(data2bs_univar(Y, segleng, segshift, epleng, maxfreqbins)); + + % required for antisymmetrization + para_yxx.chancomb = [2, 1, 1]; + [B2_yxx] = data2bs_univar([X, Y], segleng, segshift, epleng, maxfreqbins, para_yxx); + + % compute time-delays + if isempty(tde_freqbands) + [T, I] = bispectral_TD_est(B2_xxx, B2_xyx, B2_yyy, tde_method, [], 1); + [aT, aI] = bispectral_TD_est(B2_xxx, B2_xyx - B2_yxx, B2_yyy, tde_method, [], 1); % antisymmetrization + else + error('TDE on frequency bands is not implemented yet.') +% neg_frqs = cat(1, -flip(frqs(2:end)), frqs(1:end-1)); % create new freq vector with negative frequencies +% fmask = false(1, maxfreqbins); +% fmask(find(neg_frqs == tde_freqbands(1)):find(neg_frqs == tde_freqbands(2))) = 1; +% % fmask(find(neg_frqs == -tde_freqbands(2)):find(neg_frqs == -tde_freqbands(1))) = 1; +% [T, I] = bispectral_TD_est(B2_xxx, B2_xyx, B2_yyy, tde_method, fmask(1:end-1), 1); +% [aT, aI] = bispectral_TD_est(B2_xxx, B2_xyx - B2_yxx, B2_yyy, tde_method, fmask(1:end-1), 1); % antisymmetrization + end +end \ No newline at end of file diff --git a/simulations/sim_tde_bispec.m b/simulations/sim_tde_bispec.m new file mode 100644 index 0000000..c3837d2 --- /dev/null +++ b/simulations/sim_tde_bispec.m @@ -0,0 +1,87 @@ +% Simulate artificial two-channel dataset to estimate time delays using +% antisymmetrized bispectra. The times series consists of a signal and noise +% component that can be specified. + +%% Setup +clear +clc + +srate = 100; % Hz +t = 130 * srate; % length of signal in datapoints (130 seconds) +seglen = 2 * srate; % segment length in datapoints (2 seconds) +delay = 50; % delay in time bins + +param = 0.7 ; % parameter used to define distribution. lambda for +theta = [0.7 0.7]; + +signal_type = 'exp'; % ["exp", "G", "pink", "pinksq", "spiky"] +noise_type = 'G'; % ["exp", "G", "pink", "pinksq", "spiky"] + +snr = 0.8; +beta = 1; + +nboot = []; % leave to allow plotting +method = 1; % 1:4, chose which bispectral TDE should be displayed, check documentation in bispectral_TD_est + +%% Create data +eeglab +Xraw = mk_series(t+abs(delay), signal_type, 0.7,0); +noiseX = mk_series(t, noise_type, 0.7, 0); +noiseY = mk_series(t, noise_type, 0.7, 0); + +% Combine and segment +[Xlong, Ylong] = combine_sn(Xraw, noiseX, noiseY, delay, srate, snr, theta, beta); +X = segment_data(Xlong, seglen); +Y = segment_data(Ylong, seglen); + +% Unmixed noise comparison +[XlongU, YlongU] = combine_sn(Xraw, noiseX, noiseY, delay, srate, snr,[0 0], beta); + +%% Estimate bispectral TDE +ndat = seglen * 2; +segshift = floor(ndat/2); +epleng = ndat; +fres = srate; +maxfreqbins = floor(ndat/2); + +[B2_xxx] = squeeze(data2bs_univar(Xlong', ndat, segshift, epleng, maxfreqbins)); +para_xyx.chancomb = [1, 2, 1]; +[B2_xyx] = data2bs_univar([Xlong', Ylong'], ndat, segshift, epleng, maxfreqbins, para_xyx); +[B2_yyy] = squeeze(data2bs_univar(Ylong', ndat, segshift, epleng, maxfreqbins)); + +% required for antisymmetrization +para_yxx.chancomb = [2, 1, 1]; +[B2_yxx] = data2bs_univar([Xlong', Ylong'], ndat, segshift, epleng, maxfreqbins, para_yxx); + +[T, I] = bispectral_TD_est(B2_xxx, B2_xyx, B2_yyy, method, [], 1); +[aT, aI] = bispectral_TD_est(B2_xxx, B2_xyx - B2_yxx, B2_yyy, method, [], 1); + +% % TDE on frequency bands +% frqs = sfreqs(fres, srate); +% neg_frqs = cat(1, -flip(frqs(2:end)), frqs(1:end-1)); % create new freq vector with negative frequencies +% band = [20 30]; +% fmask = false(1, maxfreqbins); +% fmask(find(neg_frqs == band(1)):find(neg_frqs == band(2))) = 1; fmask(find(neg_frqs == -band(2)):find(neg_frqs == -band(1))) = 1; +% [T_, I_] = bispectral_TD_est(B2_xxx, B2_xyx, B2_yyy, method, fmask(1:end-1), 1); + + +%% Plotting +% extract estimated delay/peak +shift = (-seglen+1:seglen-1) / srate; +[peak_val, peak_idx] = max(T); +est_delay = shift(peak_idx); % in Hz + +figure; plot(shift, T, 'black') +xline(est_delay, '--r') +xlabel('Time (s)') +ylabel('a.u.') +title(sprintf('TDE | Method %d', method)) +subtitle("\tau = " + num2str(est_delay) + " (s)") +grid on + +%% Check if the estimated delay is in fact the true simulated delay +if est_delay == (delay / srate) + disp('The estimated delay matches the true simulated delay.') +else + error('The estimated delay does not match the true simulated delay') +end diff --git a/test_pipes/pipeline_connectivity.m b/test_pipes/pipeline_connectivity.m index 965d46a..5956c3a 100644 --- a/test_pipes/pipeline_connectivity.m +++ b/test_pipes/pipeline_connectivity.m @@ -49,6 +49,6 @@ EEG = pop_roi_connect(EEG, 'methods', measures(iMeasure)); t(iMeasure) = toc; end -pop_roi_connectplot(EEG, 'measure', 'MIM', 'plotcortex', 'on', 'plotmatrix', 'on'); +pop_roi_connectplot(EEG, 'measure', 'MIM', 'plotcortex', 'on', 'plotmatrix', 'on', 'plotbutterfly', 'on'); pop_roi_connectplot(EEG, 'measure', 'MIM', 'plot3d', 'on', 'plotcortex', 'off', 'plot3dparams', {'thresholdper', 0.05, 'brainmovieopt' {'nodeColorDataRange', [], 'nodeSizeLimits', [0 0.2]}}); diff --git a/test_pipes/test_pac.m b/test_pipes/test_pac.m index 994ff94..9ab4874 100644 --- a/test_pipes/test_pac.m +++ b/test_pipes/test_pac.m @@ -34,19 +34,20 @@ EEG1 = pop_roi_connect(EEG, 'methods', {'PAC', 'MIM', 'COH'}, 'fcomb', fcomb); % test all 3 connectivity functions (data2spwctrgc, data2strgcmim, roi_pac) tic -EEG2 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'conn_stats', 'on', 'nshuf', 3, 'poolsize', 12); % compute only b_anti, b_anti_norm +EEG2 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'conn_stats', 'on', 'nshuf', 3, 'poolsize', 2); % 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 +EEG3 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'bs_outopts', 5); -if ~isequal(squeeze(EEG2.roi.PAC.b_orig(:, :, 1)), EEG3.roi.PAC.b_orig) +% test if first shuffle equals the true value +tolerance = 1e-7; +if ~all(ismembertol(EEG1.roi.PAC.b_orig, squeeze(EEG2.roi.PAC.b_orig(:, :, 1)), tolerance), 'all') 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) +if ~all(ismembertol(EEG1.roi.PAC.b_anti, squeeze(EEG2.roi.PAC.b_anti(:, :, 1)), tolerance), 'all') error 'The first shuffle in the surrogate connectivity array is not the true matrix.' end - %% Test bispectrum for frequency band inputs low = [4 8]; high = [48 50]; @@ -58,27 +59,29 @@ 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', 3, 'bs_outopts', 1); +EEG5 = pop_roi_connect(EEG, 'methods', {'PAC'}, 'fcomb', fcomb, 'conn_stats', 'off', 'nshuf', 3, 'bs_outopts', 5); toc -if ~isequal(squeeze(EEG5.roi.PAC.b_anti(:, :, 1)), EEG4.roi.PAC.b_anti) +% test if first shuffle equals the true value +tolerance = 1e-7; +if ~all(ismembertol(EEG4.roi.PAC.b_anti, squeeze(EEG5.roi.PAC.b_anti(:, :, 1)), tolerance), 'all') 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(EEG3, 'measure', 'pac', 'plotmatrix', 'on'); -pop_roi_connectplot(EEG3, 'measure', 'pac_anti', 'plotmatrix', 'on'); +pop_roi_connectplot(EEG1, 'measure', 'PAC', 'plotmatrix', 'on'); +pop_roi_connectplot(EEG1, 'measure', 'PAC_anti', 'plotmatrix', 'on'); % Provoke errors by plotting bispectral tensors that do not exist -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 +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(EEG5, 'measure', 'pac', 'plotmatrix', 'on'); -pop_roi_connectplot(EEG5, 'measure', 'pac_anti', 'plotmatrix', 'on'); +pop_roi_connectplot(EEG4, 'measure', 'PAC', 'plotmatrix', 'on'); +pop_roi_connectplot(EEG4, 'measure', 'PAC_anti', 'plotmatrix', 'on'); -% plot MIM and COH as a sanity check -pop_roi_connectplot(EEG1, 'measure', 'mim', 'plotmatrix', 'on'); -pop_roi_connectplot(EEG1, 'measure', 'coh', 'plotmatrix', 'on'); +% plot MIM and aCOH as a sanity check +pop_roi_connectplot(EEG1, 'measure', 'MIM', 'plotmatrix', 'on'); +pop_roi_connectplot(EEG1, 'measure', 'aCOH', 'plotmatrix', 'on'); diff --git a/test_pipes/test_tde.m b/test_pipes/test_tde.m new file mode 100644 index 0000000..cd4517b --- /dev/null +++ b/test_pipes/test_tde.m @@ -0,0 +1,53 @@ +% Test bispectral time-delay estimation. +%% Run pipeline +clear +eeglab + +eeglabp = fileparts(which('eeglab.m')); +EEG = pop_loadset('filename','eeglab_data_epochs_ica.set','filepath',fullfile(eeglabp, 'sample_data/')); +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); +eeglab redraw; + +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] ); + +% 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); + +EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'functions','supportfiles','head_modelColin27_5003_Standard-10-5-Cap339.mat'), ... + '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); + +%% Compute iCOH to extract a promising region combination that shows connectivity +EEG = pop_roi_connect(EEG, 'methods', {'iCOH'}); +pop_roi_connectplot(EEG, 'measure', 'iCOH', 'plotcortex', 'off', 'plotmatrix', 'on', 'plotbutterfly', 'on'); % look at broadband iCOH +pop_roi_connectplot(EEG, 'measure', 'iCOH', 'freqrange', [11 13], 'plotcortex', 'off', 'plotmatrix', 'on', 'plotbutterfly', 'on'); + +%% Test TDE on broadband +% find regions from "visual analysis" above +roi1 = 'lingual L'; +roi2 = 'lateraloccipital L'; +labels = replace_underscores(get_labels(EEG)); +roi1_idx = find(ismember(labels, roi1)); +roi2_idx = find(ismember(labels, roi2)); + +EEG1 = pop_roi_connect(EEG, 'methods', {'TDE'}, 'tde_regions', [roi1_idx roi2_idx], 'tde_method', 1); + +%% Test TDE on frequency bands +% not fully implemented yet, needs further testing +% EEG2 = pop_roi_connect(EEG, 'methods', {'TDE'}, 'tde_regions', [roi1_idx roi2_idx], 'tde_method', 1, 'tde_freqbands', [11 13]); + +%% Plotting +% broadband +pop_roi_connectplot(EEG1, 'measure', 'tde'); +pop_roi_connectplot(EEG1, 'measure', 'tde_anti'); + +% % frequency band +% pop_roi_connectplot(EEG2, 'measure', 'tde'); +% pop_roi_connectplot(EEG2, 'measure', 'tde_anti'); \ No newline at end of file diff --git a/utils/get_labels.m b/utils/get_labels.m new file mode 100644 index 0000000..a0e6579 --- /dev/null +++ b/utils/get_labels.m @@ -0,0 +1,16 @@ +function labels = get_labels(EEG) + % retrieve labels from atlas + labels = strings(1,length(EEG.roi.atlas.Scouts)); + for i = 1:length(labels) + scout = struct2cell(EEG.roi.atlas.Scouts(i)); + labels(i) = char(scout(1)); + end + labels = cellstr(labels); + + % remove region labels that were not selected + if isfield(EEG.roi, 'roi_selection') + if ~isempty(EEG.roi.roi_selection) + labels = labels(cell2mat(EEG.roi.roi_selection)); + end + end +end diff --git a/utils/replace_underscores.m b/utils/replace_underscores.m new file mode 100644 index 0000000..55e7018 --- /dev/null +++ b/utils/replace_underscores.m @@ -0,0 +1,4 @@ +function new_labels = replace_underscores(labels) + % remove underscores in label names to avoid bug + new_labels = strrep(labels, '_', ' '); +end \ No newline at end of file