Merge pull request #89 from nguyen-td/feature

Feature

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'));

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

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

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
+

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

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

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