Merge pull request #76 from nguyen-td/master

Very nice. I reviewed the code, checked out the repository and run the test case. Everything looks fine.

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

libs/pellegrini/fp_data2bs_event_uni.m

@@ -1,4 +1,4 @@
-function [cs,nave]=fp_data2bs_event_uni(data,segleng,segshift,epleng,freqpairs,nshuf)
+function [cs,nave]=fp_data2bs_event_uni(data,segleng,segshift,epleng,freqpairs,ishuf)
 % Calculates bispectral-tensors  from data for event-related measurement
 % and their null distributions using a shuffling approach. 
 %
@@ -15,7 +15,7 @@
 %           e.g. segshift=segleng/2 makes overlapping segments
 % epleng: leng of each epoch
 % freqpairs: pairs of  frequency in bins
-% nshuf: required number of samples (shuffles) in the null distribution 
+% ishuf: one shuffle in the null distribution 
 % para: structure which is eventually used later
 %
 % output:
@@ -41,54 +41,52 @@
 nseg=floor((epleng-segleng)/segshift)+1; %total number of segments
 assert(nseg==1,'only possible with 1 segment')
 
-cs=zeros(nchan,nchan,nchan,2,nshuf);
+cs=zeros(nchan,nchan,nchan,2);
 
 %get Fourier coefficients
 coeffs = fp_fft_coeffs(data,segleng,segshift,epleng,freqpairs);
 
-for ishuf = 1:nshuf
-    nave=0;
-    csloc1=zeros(nchan,nchan,nchan);
-    csloc2=zeros(nchan,nchan,nchan);
-    cs1=zeros(nchan,nchan,nchan);
-    cs2=zeros(nchan,nchan,nchan);
+%pre-allocation
+nave=0;
+csloc1=zeros(nchan,nchan,nchan);
+csloc2=zeros(nchan,nchan,nchan);
+cs1=zeros(nchan,nchan,nchan);
+cs2=zeros(nchan,nchan,nchan);
 
-    if ishuf == 1 %the first shuffle contains the true order
-        inds = 1:nep;
-    else
-        inds = randperm(nep,nep); %indices for shuffling of epochs for channel 2 
-    end
+if ishuf == 1 %the first shuffle contains the true order
+    inds = 1:nep;
+else
+    inds = randperm(nep,nep); %indices for shuffling of epochs for channel 2 
+end
+
+for j=1:nep %loop over epochs 
+
+    %bispec of f1 (low peak), f2-f1 (left side lobe), f2 (high peak)
+    csloc1(1,1,1)=transpose(coeffs(1,1,j))        *coeffs(2,1,j)      *conj(coeffs(3,1,j));
+    csloc1(2,1,1)=transpose(coeffs(1,2,inds(j)))  *coeffs(2,1,j)      *conj(coeffs(3,1,j));
+    csloc1(1,2,1)=transpose(coeffs(1,1,j))        *coeffs(2,2,inds(j))*conj(coeffs(3,1,j));
+    csloc1(1,1,2)=transpose(coeffs(1,1,j))        *coeffs(2,1,j)      *conj(coeffs(3,2,inds(j)));
+    csloc1(2,2,1)=transpose(coeffs(1,2,inds(j)))  *coeffs(2,2,inds(j))*conj(coeffs(3,1,j));
+    csloc1(1,2,2)=transpose(coeffs(1,1,j))        *coeffs(2,2,inds(j))*conj(coeffs(3,2,inds(j)));
+    csloc1(2,2,2)=transpose(coeffs(1,2,inds(j)))  *coeffs(2,2,inds(j))*conj(coeffs(3,2,inds(j)));
+    csloc1(2,1,2)=transpose(coeffs(1,2,inds(j)))  *coeffs(2,1,j)      *conj(coeffs(3,2,inds(j)));
 
-    for j=1:nep %loop over epochs 
-
-        %bispec of f1 (low peak), f2-f1 (left side lobe), f2 (high peak)
-        csloc1(1,1,1)=transpose(coeffs(1,1,j))        *coeffs(2,1,j)      *conj(coeffs(3,1,j));
-        csloc1(2,1,1)=transpose(coeffs(1,2,inds(j)))  *coeffs(2,1,j)      *conj(coeffs(3,1,j));
-        csloc1(1,2,1)=transpose(coeffs(1,1,j))        *coeffs(2,2,inds(j))*conj(coeffs(3,1,j));
-        csloc1(1,1,2)=transpose(coeffs(1,1,j))        *coeffs(2,1,j)      *conj(coeffs(3,2,inds(j)));
-        csloc1(2,2,1)=transpose(coeffs(1,2,inds(j)))  *coeffs(2,2,inds(j))*conj(coeffs(3,1,j));
-        csloc1(1,2,2)=transpose(coeffs(1,1,j))        *coeffs(2,2,inds(j))*conj(coeffs(3,2,inds(j)));
-        csloc1(2,2,2)=transpose(coeffs(1,2,inds(j)))  *coeffs(2,2,inds(j))*conj(coeffs(3,2,inds(j)));
-        csloc1(2,1,2)=transpose(coeffs(1,2,inds(j)))  *coeffs(2,1,j)      *conj(coeffs(3,2,inds(j)));
-
-        %bispec of f1 (low peak), f2 high peak), f1+f2 (right side lobe)
-        csloc2(1,1,1)=transpose(coeffs(1,1,j))          *coeffs(3,1,j)      *conj(coeffs(4,1,j));
-        csloc2(2,1,1)=transpose(coeffs(1,2,inds(j)))    *coeffs(3,1,j)      *conj(coeffs(4,1,j));
-        csloc2(1,2,1)=transpose(coeffs(1,1,j))          *coeffs(3,2,inds(j))*conj(coeffs(4,1,j));
-        csloc2(1,1,2)=transpose(coeffs(1,1,j))          *coeffs(3,1,j)      *conj(coeffs(4,2,inds(j)));
-        csloc2(2,2,1)=transpose(coeffs(1,2,inds(j)))    *coeffs(3,2,inds(j))*conj(coeffs(4,1,j));
-        csloc2(1,2,2)=transpose(coeffs(1,1,j))          *coeffs(3,2,inds(j))*conj(coeffs(4,2,inds(j)));
-        csloc2(2,2,2)=transpose(coeffs(1,2,inds(j)))    *coeffs(3,2,inds(j))*conj(coeffs(4,2,inds(j)));
-        csloc2(2,1,2)=transpose(coeffs(1,2,inds(j)))    *coeffs(3,1,j)      *conj(coeffs(4,2,inds(j)));
-
-        cs1=cs1+csloc1;
-        cs2=cs2+csloc2;
-
-        nave=nave+1;
-    end
+    %bispec of f1 (low peak), f2 high peak), f1+f2 (right side lobe)
+    csloc2(1,1,1)=transpose(coeffs(1,1,j))          *coeffs(3,1,j)      *conj(coeffs(4,1,j));
+    csloc2(2,1,1)=transpose(coeffs(1,2,inds(j)))    *coeffs(3,1,j)      *conj(coeffs(4,1,j));
+    csloc2(1,2,1)=transpose(coeffs(1,1,j))          *coeffs(3,2,inds(j))*conj(coeffs(4,1,j));
+    csloc2(1,1,2)=transpose(coeffs(1,1,j))          *coeffs(3,1,j)      *conj(coeffs(4,2,inds(j)));
+    csloc2(2,2,1)=transpose(coeffs(1,2,inds(j)))    *coeffs(3,2,inds(j))*conj(coeffs(4,1,j));
+    csloc2(1,2,2)=transpose(coeffs(1,1,j))          *coeffs(3,2,inds(j))*conj(coeffs(4,2,inds(j)));
+    csloc2(2,2,2)=transpose(coeffs(1,2,inds(j)))    *coeffs(3,2,inds(j))*conj(coeffs(4,2,inds(j)));
+    csloc2(2,1,2)=transpose(coeffs(1,2,inds(j)))    *coeffs(3,1,j)      *conj(coeffs(4,2,inds(j)));
 
-    %shape cs: chan x chan x chan x peak_combination x shuffles 
-    cs(:,:,:,1,ishuf) = cs1./nave; 
-    cs(:,:,:,2,ishuf) = cs2./nave; 
+    cs1=cs1+csloc1;
+    cs2=cs2+csloc2;
 
-end
+    nave=nave+1;
+end
+
+%shape cs: chan x chan x chan x peak_combination x shuffle
+cs(:,:,:,1) = cs1./nave; 
+cs(:,:,:,2) = cs2./nave; 

pop_roi_connect.m

@@ -17,9 +17,9 @@
 %  'naccu'          - [integer]  Number of accumulation for stats. Default is 0.
 %  'methods'        - [cell] Cell of strings corresponding to methods.
 %                       'CS'    : Cross spectrum
-%                       'aCOH'  : Coherence
-%                       'cCOH'  : (complex-valued) coherency
-%                       'iCOH'  : absolute value of the imaginary part of coherency
+%                       'aCOH'  : Coherence 
+%                       'cCOH'  : (Complex-valued) Coherency
+%                       'iCOH'  : Absolute value of the imaginary part of Coherency
 %                       'GC'    : Granger Causality
 %                       'TRGC'  : Time-reversed Granger Causality
 %                       'wPLI'  : Weighted Phase Lag Index
@@ -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