From 725563de49bf0c2b76cd8bc7531c311ec582c163 Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Thu, 14 Mar 2024 19:57:33 +0100 Subject: [PATCH 01/12] add TDE incl. simulation script --- .gitignore | 1 + eegplugin_roiconnect.m | 1 + libs/jurhar/bispectral_TD_est.m | 51 ++++++++++++++++++++ libs/jurhar/combine_sn.m | 31 ++++++++++++ libs/jurhar/mk_delay.m | 31 ++++++++++++ libs/jurhar/mk_series.m | 38 +++++++++++++++ libs/jurhar/segment_data.m | 35 ++++++++++++++ simulations/sim_tde_bispec.m | 83 +++++++++++++++++++++++++++++++++ 8 files changed, 271 insertions(+) create mode 100644 libs/jurhar/bispectral_TD_est.m create mode 100644 libs/jurhar/combine_sn.m create mode 100644 libs/jurhar/mk_delay.m create mode 100644 libs/jurhar/mk_series.m create mode 100644 libs/jurhar/segment_data.m create mode 100644 simulations/sim_tde_bispec.m diff --git a/.gitignore b/.gitignore index 3e3461a..567bad1 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ *.asv +.DS_Store 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/simulations/sim_tde_bispec.m b/simulations/sim_tde_bispec.m new file mode 100644 index 0000000..7c7d9cd --- /dev/null +++ b/simulations/sim_tde_bispec.m @@ -0,0 +1,83 @@ +% 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; +frqs = sfreqs(2 * fres, srate); + +[B2_xxx] = squeeze(data2bs_univar(Xlong', 2 * seglen, segshift, epleng, length(frqs)-1)); +para_xyx.chancomb = [1, 2, 1]; +[B2_xyx] = data2bs_univar([Xlong', Ylong'], 2 * seglen, segshift, epleng, length(frqs)-1, para_xyx); +[B2_yyy] = squeeze(data2bs_univar(Ylong', 2 * seglen, segshift, epleng, length(frqs)-1)); + +% required for antisymmetrization +para_yxx.chancomb = [2, 1, 1]; +[B2_yxx] = data2bs_univar([Xlong', Ylong'], 2 * seglen, segshift, epleng, length(frqs)-1, 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); + +% test band +band = [20 30]; +fmask = true(1, length(frqs)); +fmask(frqs < band(1) | frqs > band(2)) = 0; +[T_, I_] = bispectral_TD_est(B2_xxx, B2_xyx, B2_yyy, method, fmask(1:end-1), 1); + + +%% Plotting +delay_scale = (-seglen+1:seglen-1) / srate; +[peak_val, peak_idx] = max(aT); % compute estimated delay/peak +est_delay = delay_scale(peak_idx); + +figure; plot(delay_scale, aT) +xline(est_delay, '--r') +xlabel('Time (s)') +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 + error('The estimated delay does not match the true simulated delay') +else + disp('The estimated delay matches the true simulated delay.') +end From c013f6b3dd7418eba30bda2f145240fbd2cd4509 Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Thu, 14 Mar 2024 20:02:19 +0100 Subject: [PATCH 02/12] update comment --- simulations/sim_tde_bispec.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/simulations/sim_tde_bispec.m b/simulations/sim_tde_bispec.m index 7c7d9cd..feece70 100644 --- a/simulations/sim_tde_bispec.m +++ b/simulations/sim_tde_bispec.m @@ -56,7 +56,7 @@ [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); -% test band +% TDE on frequency bands band = [20 30]; fmask = true(1, length(frqs)); fmask(frqs < band(1) | frqs > band(2)) = 0; @@ -65,7 +65,7 @@ %% Plotting delay_scale = (-seglen+1:seglen-1) / srate; -[peak_val, peak_idx] = max(aT); % compute estimated delay/peak +[peak_val, peak_idx] = max(aT); % extract estimated delay/peak est_delay = delay_scale(peak_idx); figure; plot(delay_scale, aT) From 63aed0bb8b653a061ae3dbee28c48283e20353af Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Tue, 19 Mar 2024 17:12:29 +0100 Subject: [PATCH 03/12] typo --- roi_pac.m | 2 +- simulations/sim_tde_bispec.m | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/roi_pac.m b/roi_pac.m index 7a333c5..9f65e02 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: diff --git a/simulations/sim_tde_bispec.m b/simulations/sim_tde_bispec.m index feece70..06e4a43 100644 --- a/simulations/sim_tde_bispec.m +++ b/simulations/sim_tde_bispec.m @@ -68,9 +68,10 @@ [peak_val, peak_idx] = max(aT); % extract estimated delay/peak est_delay = delay_scale(peak_idx); -figure; plot(delay_scale, aT) +figure; plot(delay_scale, aT, '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 From 5a081c3db15b3283fdcc757aeb3b700b780bef23 Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Tue, 19 Mar 2024 23:10:12 +0100 Subject: [PATCH 04/12] update docs --- pop_roi_connect.m | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 2adc07b..9fd4525 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,10 @@ % 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_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. From 0664da7196ca37688a788d791e90e176dd028213 Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Tue, 19 Mar 2024 23:49:59 +0100 Subject: [PATCH 05/12] move methods check from roi_connect to pop_roi_connect --- pop_roi_connect.m | 12 ++++++++++++ roi_connect.m | 8 ++++---- roi_pac.m | 2 +- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index 9fd4525..b853748 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -55,6 +55,7 @@ % 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). +% 'tde_method' - [integer] TDE method, must be between 1:4, open bispectral_TD_est.m for details. Default is 1. % '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. @@ -186,6 +187,9 @@ 'fcomb' 'struct' { } struct; 'bs_outopts' 'integer' { } 1; 'roi_selection' 'cell' { } { }; + 'tde_regions' 'integer' { } [ ]; + 'tde_freqbands' 'integer' { } [ ]; + 'tde_method' 'integer' { 1:4 } 1; 'conn_stats' 'string' { } 'off'; ... 'nshuf' 'integer' { } 1001; ... 'poolsize' 'integer' { } 1}, 'pop_roi_connect'); @@ -202,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); @@ -345,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/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 9f65e02..f553097 100644 --- a/roi_pac.m +++ b/roi_pac.m @@ -38,7 +38,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 From 21ebb265d2d12ab05509ea66529d1006d0553fed Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Wed, 20 Mar 2024 00:10:32 +0100 Subject: [PATCH 06/12] fix simulation script, add TDE to pop_roi_connect --- pop_roi_connect.m | 4 ++-- simulations/sim_tde_bispec.m | 19 ++++++++++--------- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/pop_roi_connect.m b/pop_roi_connect.m index b853748..328457c 100644 --- a/pop_roi_connect.m +++ b/pop_roi_connect.m @@ -51,11 +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). -% 'tde_method' - [integer] TDE method, must be between 1:4, open bispectral_TD_est.m for details. Default is 1. % '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. @@ -187,9 +187,9 @@ 'fcomb' 'struct' { } struct; 'bs_outopts' 'integer' { } 1; 'roi_selection' 'cell' { } { }; + 'tde_method' 'integer' { 1:4 } 1; 'tde_regions' 'integer' { } [ ]; 'tde_freqbands' 'integer' { } [ ]; - 'tde_method' 'integer' { 1:4 } 1; 'conn_stats' 'string' { } 'off'; ... 'nshuf' 'integer' { } 1001; ... 'poolsize' 'integer' { } 1}, 'pop_roi_connect'); diff --git a/simulations/sim_tde_bispec.m b/simulations/sim_tde_bispec.m index 06e4a43..12c7f00 100644 --- a/simulations/sim_tde_bispec.m +++ b/simulations/sim_tde_bispec.m @@ -44,14 +44,14 @@ fres = srate; frqs = sfreqs(2 * fres, srate); -[B2_xxx] = squeeze(data2bs_univar(Xlong', 2 * seglen, segshift, epleng, length(frqs)-1)); +[B2_xxx] = squeeze(data2bs_univar(Xlong', ndat, segshift, epleng, length(frqs)-1)); para_xyx.chancomb = [1, 2, 1]; -[B2_xyx] = data2bs_univar([Xlong', Ylong'], 2 * seglen, segshift, epleng, length(frqs)-1, para_xyx); -[B2_yyy] = squeeze(data2bs_univar(Ylong', 2 * seglen, segshift, epleng, length(frqs)-1)); +[B2_xyx] = data2bs_univar([Xlong', Ylong'], ndat, segshift, epleng, length(frqs)-1, para_xyx); +[B2_yyy] = squeeze(data2bs_univar(Ylong', ndat, segshift, epleng, length(frqs)-1)); % required for antisymmetrization para_yxx.chancomb = [2, 1, 1]; -[B2_yxx] = data2bs_univar([Xlong', Ylong'], 2 * seglen, segshift, epleng, length(frqs)-1, para_yxx); +[B2_yxx] = data2bs_univar([Xlong', Ylong'], ndat, segshift, epleng, length(frqs)-1, 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); @@ -64,9 +64,10 @@ %% Plotting +% extract estimated delay/peak delay_scale = (-seglen+1:seglen-1) / srate; -[peak_val, peak_idx] = max(aT); % extract estimated delay/peak -est_delay = delay_scale(peak_idx); +[peak_val, peak_idx] = max(T); +est_delay = delay_scale(peak_idx); % in Hz figure; plot(delay_scale, aT, 'black') xline(est_delay, '--r') @@ -77,8 +78,8 @@ grid on %% Check if the estimated delay is in fact the true simulated delay -if ~est_delay == delay - error('The estimated delay does not match the true simulated delay') -else +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 From 3d457f613ee2a89ec016651431af90a5c47435dc Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Wed, 20 Mar 2024 13:54:46 +0100 Subject: [PATCH 07/12] update documentation --- roi_pac.m | 9 +++++---- simulations/sim_tde_bispec.m | 5 +++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/roi_pac.m b/roi_pac.m index f553097..43ad1a9 100644 --- a/roi_pac.m +++ b/roi_pac.m @@ -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 diff --git a/simulations/sim_tde_bispec.m b/simulations/sim_tde_bispec.m index 12c7f00..3ef6edb 100644 --- a/simulations/sim_tde_bispec.m +++ b/simulations/sim_tde_bispec.m @@ -27,7 +27,7 @@ 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); +noiseY = mk_series(t, nois50e_type, 0.7, 0); % Combine and segment [Xlong, Ylong] = combine_sn(Xraw, noiseX, noiseY, delay, srate, snr, theta, beta); @@ -43,6 +43,7 @@ epleng = ndat; fres = srate; frqs = sfreqs(2 * fres, srate); +maxfreqbins = floor(ndat/2); [B2_xxx] = squeeze(data2bs_univar(Xlong', ndat, segshift, epleng, length(frqs)-1)); para_xyx.chancomb = [1, 2, 1]; @@ -66,7 +67,7 @@ %% Plotting % extract estimated delay/peak delay_scale = (-seglen+1:seglen-1) / srate; -[peak_val, peak_idx] = max(T); +[peak_val, peak_idx] = max(aT); est_delay = delay_scale(peak_idx); % in Hz figure; plot(delay_scale, aT, 'black') From 27e6945ed7f9c37df534eaf2bee350f95776d33f Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Wed, 20 Mar 2024 15:03:11 +0100 Subject: [PATCH 08/12] flip y-axis for plotting --- pop_roi_connectplot.m | 43 ++++++++++++++++++------------------ simulations/sim_tde_bispec.m | 6 ++--- 2 files changed, 25 insertions(+), 24 deletions(-) diff --git a/pop_roi_connectplot.m b/pop_roi_connectplot.m index 1b720ee..f5222e4 100644 --- a/pop_roi_connectplot.m +++ b/pop_roi_connectplot.m @@ -572,27 +572,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,6 +832,7 @@ 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 + set(gca, 'YDir','normal') h = title([ 'ROI to ROI ' upper(replace_underscores(measure)) ' (' titleStr ')' ]); set(h, 'fontsize', 16); xtickangle(90) diff --git a/simulations/sim_tde_bispec.m b/simulations/sim_tde_bispec.m index 3ef6edb..da63533 100644 --- a/simulations/sim_tde_bispec.m +++ b/simulations/sim_tde_bispec.m @@ -66,11 +66,11 @@ %% Plotting % extract estimated delay/peak -delay_scale = (-seglen+1:seglen-1) / srate; +shift = (-seglen+1:seglen-1) / srate; [peak_val, peak_idx] = max(aT); -est_delay = delay_scale(peak_idx); % in Hz +est_delay = shift(peak_idx); % in Hz -figure; plot(delay_scale, aT, 'black') +figure; plot(shift, aT, 'black') xline(est_delay, '--r') xlabel('Time (s)') ylabel('a.u.') From 68c00e1b78b1d7927e65e89bd47e864507e53b71 Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Wed, 20 Mar 2024 17:00:38 +0100 Subject: [PATCH 09/12] add butterfly plot --- pop_roi_connectplot.m | 44 ++++++++++++++++++++++++------ test_pipes/pipeline_connectivity.m | 2 +- utils/get_labels.m | 16 +++++++++++ utils/replace_underscores.m | 4 +++ 4 files changed, 56 insertions(+), 10 deletions(-) create mode 100644 utils/get_labels.m create mode 100644 utils/replace_underscores.m diff --git a/pop_roi_connectplot.m b/pop_roi_connectplot.m index f5222e4..3e59580 100644 --- a/pop_roi_connectplot.m +++ b/pop_roi_connectplot.m @@ -32,6 +32,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' @@ -162,7 +163,7 @@ 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; @@ -327,6 +328,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 +404,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 +423,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 +436,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 +456,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 = []; @@ -519,15 +528,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 +548,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; @@ -833,7 +859,7 @@ function roi_plotcoloredlobes( EEG, matrix, titleStr, measure, hemisphere, group 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 set(gca, 'YDir','normal') - h = title([ 'ROI to ROI ' upper(replace_underscores(measure)) ' (' titleStr ')' ]); + h = title([ 'ROI to ROI ' replace_underscores(measure) ' (' titleStr ')' ]); set(h, 'fontsize', 16); xtickangle(90) 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/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 From 86b43934dd01387607bdabbe736e98856c4a6b65 Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Wed, 20 Mar 2024 17:48:12 +0100 Subject: [PATCH 10/12] add TDE --- pop_roi_connectplot.m | 65 ++++++++++++++++++++++++---- roi_tde.m | 98 +++++++++++++++++++++++++++++++++++++++++++ test_tde.m | 52 +++++++++++++++++++++++ 3 files changed, 208 insertions(+), 7 deletions(-) create mode 100644 roi_tde.m create mode 100644 test_tde.m diff --git a/pop_roi_connectplot.m b/pop_roi_connectplot.m index 3e59580..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' @@ -156,7 +158,7 @@ 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 @@ -167,7 +169,7 @@ 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 @@ -178,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 @@ -189,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 @@ -200,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 @@ -247,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 @@ -495,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 @@ -928,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_tde.m b/roi_tde.m new file mode 100644 index 0000000..2f28307 --- /dev/null +++ b/roi_tde.m @@ -0,0 +1,98 @@ +% 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(2 * 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 + fmask = true(1, length(frqs)); + fmask(frqs < tde_freqbands(1) | frqs > tde_freqbands(2)) = 0; + [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/test_tde.m b/test_tde.m new file mode 100644 index 0000000..b72bf26 --- /dev/null +++ b/test_tde.m @@ -0,0 +1,52 @@ +% 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 +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 From 035b4f3eff77bddf0b557b5362bd768aef596b37 Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Thu, 21 Mar 2024 11:04:15 +0100 Subject: [PATCH 11/12] remove tde on freq bands for now --- roi_tde.m | 13 ++++++++----- simulations/sim_tde_bispec.m | 27 ++++++++++++++------------- test_tde.m => test_pipes/test_tde.m | 9 +++++---- 3 files changed, 27 insertions(+), 22 deletions(-) rename test_tde.m => test_pipes/test_tde.m (88%) diff --git a/roi_tde.m b/roi_tde.m index 2f28307..4281119 100644 --- a/roi_tde.m +++ b/roi_tde.m @@ -44,7 +44,7 @@ segleng = ndat; epleng = ndat; fres = EEG.srate; - frqs = sfreqs(2 * fres, EEG.srate); + frqs = sfreqs(fres, EEG.srate); maxfreqbins = floor(segleng/2); % only keeep first PC @@ -90,9 +90,12 @@ [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 - fmask = true(1, length(frqs)); - fmask(frqs < tde_freqbands(1) | frqs > tde_freqbands(2)) = 0; - [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 + 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 index da63533..c3837d2 100644 --- a/simulations/sim_tde_bispec.m +++ b/simulations/sim_tde_bispec.m @@ -27,7 +27,7 @@ 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, nois50e_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); @@ -42,35 +42,36 @@ segshift = floor(ndat/2); epleng = ndat; fres = srate; -frqs = sfreqs(2 * fres, srate); maxfreqbins = floor(ndat/2); -[B2_xxx] = squeeze(data2bs_univar(Xlong', ndat, segshift, epleng, length(frqs)-1)); +[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, length(frqs)-1, para_xyx); -[B2_yyy] = squeeze(data2bs_univar(Ylong', ndat, segshift, epleng, length(frqs)-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, length(frqs)-1, para_yxx); +[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 -band = [20 30]; -fmask = true(1, length(frqs)); -fmask(frqs < band(1) | frqs > band(2)) = 0; -[T_, I_] = bispectral_TD_est(B2_xxx, B2_xyx, B2_yyy, method, fmask(1:end-1), 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(aT); +[peak_val, peak_idx] = max(T); est_delay = shift(peak_idx); % in Hz -figure; plot(shift, aT, 'black') +figure; plot(shift, T, 'black') xline(est_delay, '--r') xlabel('Time (s)') ylabel('a.u.') diff --git a/test_tde.m b/test_pipes/test_tde.m similarity index 88% rename from test_tde.m rename to test_pipes/test_tde.m index b72bf26..cd4517b 100644 --- a/test_tde.m +++ b/test_pipes/test_tde.m @@ -40,13 +40,14 @@ EEG1 = pop_roi_connect(EEG, 'methods', {'TDE'}, 'tde_regions', [roi1_idx roi2_idx], 'tde_method', 1); %% Test TDE on frequency bands -EEG2 = pop_roi_connect(EEG, 'methods', {'TDE'}, 'tde_regions', [roi1_idx roi2_idx], 'tde_method', 1, 'tde_freqbands', [11 13]); +% 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 +% % frequency band +% pop_roi_connectplot(EEG2, 'measure', 'tde'); +% pop_roi_connectplot(EEG2, 'measure', 'tde_anti'); \ No newline at end of file From 0e56d114d75ef6b65f1ab06aa0dcc74ee457adbd Mon Sep 17 00:00:00 2001 From: Tien Dung Nguyen Date: Thu, 21 Mar 2024 17:41:03 +0100 Subject: [PATCH 12/12] update unit test to avoid floating point errors --- .gitignore | 2 -- test_pipes/test_pac.m | 35 +++++++++++++++++++---------------- 2 files changed, 19 insertions(+), 18 deletions(-) delete mode 100644 .gitignore diff --git a/.gitignore b/.gitignore deleted file mode 100644 index 567bad1..0000000 --- a/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -*.asv -.DS_Store 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');