add TDE incl. simulation script

.gitignore

@@ -1 +1,2 @@
 *.asv
+.DS_Store

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

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