debugged, commented rtCCA and disabled flagICRegressors. Removed debug plot from GLM function

    % Set other data blocks

    beta_blks{iBlk} = beta;

    yR_blks{iBlk}   = yR;

    %

    % debug

    figure

    plot (squeeze(yavg(:,2,:,2)))

end

% SYNTAX:

% [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, runIdx, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

% [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, runIdx, flagtCCA, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

%

% UI NAME:

% hmrR_tCCA

%

% DESCRIPTION:

% This script generates regressors using the regularized temporally embedded

% Canonical Correlation Anlaysis described in detail in von Lhmann, NeuroImage,

% 2020. If you use this method, please use the following CITATION:

% Canonical Correlation Anlaysis described in detail in von Lhmann et al. (2020), NeuroImage. 

% If you use this method, please use the following CITATION:

% von Lhmann, Alexander, et al. "Improved physiological noise regression

% in fNIRS: A multimodal extension of the General Linear Model using temporally

% embedded Canonical Correlation Analysis." NeuroImage 208 (2020): 116472.

% probe - SNIRF probe type containing source/detector geometry data (See SNIRF Spec for more details)

% runIdx - the index of the run in a multi-run session

% flagtCCA - turns the function on / off

% flagICRegressors - selects regressor generation strategy. (0/false) chooses a common set

%           of regressors for all fNIRS channels. (1/true) generates one

%           individual regressor per fNIRS channel.

% tCCAparams - These are the parameters for tCCA function

%            1 - 1 timelag [s]

%            2 - 2 step size [samples]

%

% OUTPUTS:

% Aaux - A matrix of auxilliary regressors (#time points x #Aux regressors)

% rcMap - An array of cells (1 x #fNIRS channels) containing

%           aux regressor indices for individual regressor-channel mapping.

%           Only relevant when flagICRegressors = 1, otherwise rcMap is empty.

% rcMap - Currently always empty or "all". 

%           Under development: Will also provide an array of cells (1 x #fNIRS channels) 

%           containing aux regressor indices for individual regressor-channel mapping.

%           Only relevant when flagICRegressors = 1.

%

% USAGE OPTIONS:

% hmrR_tCCA_Concentration_Data: [Aaux, rcMap] = hmrR_tCCA(dc, aux, probe, iRun, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

% hmrR_tCCA_OD_Data: [Aaux, rcMap] = hmrR_tCCA(dod, aux, probe, iRun, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

% hmrR_tCCA_Concentration_Data: [Aaux, rcMap] = hmrR_tCCA(dc, aux, probe, iRun, flagtCCA, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

% hmrR_tCCA_OD_Data: [Aaux, rcMap] = hmrR_tCCA(dod, aux, probe, iRun, flagtCCA, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

%

%

% PARAMETERS:

% flagtCCA: 1

% flagICRegressors: 0

% tCCAparams: [3 2 0.3]

% tCCAaux_inx: [1 2 3 4 5 6 7 8]

% rhoSD_ssThresh: 15.0

% runIdxResting: 1

% tResting: [30 210]

%

function [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, runIdx, flagtCCA, flagICRegressors, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

function [Aaux, rcMap] = hmrR_tCCA(data, aux, probe, runIdx, flagtCCA, tCCAparams, tCCAaux_inx, rhoSD_ssThresh, runIdxResting, tResting)

%% COMMENTS/THOUGHTS/QUESTIONS ALEX

% 2) Output canonical correlation coefficients as quality metric?

% 3) Output fNIRS signal(s)/Aux regressors for visualization/quality control?

% 4) single channel (regressor) feature: regressor channel map currently assumes first half HbO, secon half HbR. CHECK HOW NEEDS TO BE ADAPTED FOR GLM?

% 5) Implement the variable low/bandpass filter coefficients from previous processing stream !!!

%%

%% flags and tCCA settings

flags.pcaf =  [0 0]; % no pca of X or AUX

% flags.shrink = true; % perform shrinkage in the CCA

flags.shrink = false;  % JD: set to false because rtcca generate error, "Undefined function or variable 'cshrink'".

% flagICRegressors - selects regressor generation strategy. (0/false) chooses a common set

% of regressors for all fNIRS channels. (1/true) generates one individual

% regressor per fNIRS channel. This feature is currently under development,

% please do not use it.

flagICRegressors = 0;

% tCCA parameters

param.tau = tCCAparams(2); %stepsize for embedding in samples (tune to sample frequency!)

timelag = tCCAparams(1);

        switch dotCCA

            case 'train'

                %% if the tCCAfilter variable is not existing, learn and save it (this is the training/resting run)

                %       If flagICRegressors = 1,

                %       then each column corresponds to an indiviudal channel

                %       regressor. Otherwise, columns correspond to common regressors

                %       Columns of the tCCA filter matrix correspond to common regressors

                %       for all channels in descending order ranked by canonical correlation coefficient

                %       number of embeddings

                %       number of embeddings. If flagICRegressors = 1 (currently N/A), 

                %       each column corresponds to an indiviudal channel regressor. 

                % cut data to selected time window before training

                % warning if resting segment is shorter than 1min

                d_long = d_long(cIdx,:);

                AUX = AUX(cIdx,:);

                if flagICRegressors %% learn channel specific filters and regressors

                    for cc = 1:size(d_long,2)

                        %% perform tCCA with shrinkage for each fNIRS channel

                        [REG,  ADD] = rtcca_dummy(d_long(:,cc),AUX,param,flags);

                        %% save individual regressor and channel map

                        % return the reduced number of available regressors

                        Aaux{iBlk}(:,cc) = REG(:,1);

                        % assemble tCCA filter for channel specific regressors from individual reduced mapping matrices Av

                        tCCAfilter(:,cc) = ADD.Av(:,1);

                        % save channel-regressor map. First half HbO, second half HbR

                        rcMap{iBlk}{cc} = cc;

                    end

                    % reshape (HbO first row, HbR second row)

                    rcMap = reshape(rcMap,[2,numel(rcMap)/2]);

                else %% learn common set of regressors for all channels (default)

                if ~flagICRegressors %% learn common set of regressors for all channels (default)

                    %% perform tCCA with shrinkage for all fNIRS channels

                    [REG,  ADD_trn] = rtcca(d_long,AUX,param,flags);                 %% save and output tCCA filter.

                    %reduce filter matrix with the help of correlation threshold or max number of regressors

                    if ctr < 1

                        % use only auxiliary tcca components that have correlation > ct

                        compindex=find(ADD_trn.ccac>param.ct);

                        %% JD had a comment here as Undefined variable "ADD_trn" or class "ADD_trn.ccac"

                    else

                        % use only the first ctr auxiliary tcca components (fixed number of

                        % regressors = ctr)

                    tCCAfilter = ADD_trn.Av(:,compindex);

                    % set channel-regressor map to empty (GLM will use all available regressors for all channels)

                    rcMap{iBlk} = 'all';

                else %% learn channel specific filters and regressors. UNDER DEVELOPMENT - DO NOT USE

                    for cc = 1:size(d_long,2)/2

                        %% perform tCCA with shrinkage for each fNIRS channel (always HbO and HbR together)

                        [REG,  ADD] = rtcca(d_long(:,cc),AUX,param,flags);

                        %% save individual regressor and channel map

                        % return the reduced number of available regressors

                        Aaux{iBlk}(:,cc) = REG(:,1);

                        % assemble tCCA filter for channel specific regressors from individual reduced mapping matrices Av

                        tCCAfilter(:,cc) = ADD.Av(:,1);

                        % save channel-regressor map. First half HbO, second half HbR

                        rcMap{iBlk}{cc} = cc;

                    end

                    % reshape (HbO first row, HbR second row)

                    rcMap{iBlk} = reshape(rcMap{iBlk},[2,numel(rcMap{iBlk})/2]);

                end

                Aaux = []; % we do not want to use regressors for the resting run

                Aaux = []; % we do not want to use regressors for the resting run (available only for investigative purposes)

                %% save tCCAfilter matrix that was learned from resting state data.

                fprintf('hmrR_tCCA: run idx = %d. Generated and Saved tCCAfilter\n', runIdx)

                save(filterFilename, '-ascii', 'tCCAfilter');

function [ REG, ADD ] = rtcca_dummy( X, AUX, param, flags )

% rtcca uses multimodal data, here fNIRS, acceleration and auxiliary signals, to

% extract regressors for GLM-based noise reduction using temporally embedded

% CCA regularized with shrinkage of covariance matrices with 'optimal' parameter from Ledoit & Wolf 2004

% This function is an adaptation of the BLISSARD artifact rejection code by

% Alexander von Luhmann, written in 2018 at TUB

% Main difference: CCA is performed with shrinked covariance matrices and

% is therefore implemented as generalized eigenvalue problem (in previous

% versions Matlab's canoncorr was used)

%

% INPUTS:

% X:    input nirs data with dimensions |time x channels|

% AUX:   input accel data with dimensions |time x channels|

%                    feeding in orthogonalized (PCA) data is advised.

% param.tau:  temporal embedding parameter (lag in samples)

% param.NumOfEmb: Number of temporally embedded copies

% param.ct:   correlation threshold. returns only those regressors that

%               have a canonical correlation greater than the threshold.

% flags.pcaf:  flags for performing pca of [AUX X] as preprocessing step.

%       default: [0 0]

% flags.shrink: regularize CCA with automatic shrinkage of covariance

%       matrices. default: true

%

% OUTPUTS:

% REG: Regressors found by temp. embedded CCA

% ADD.ccac: CCA correlation coefficients between projected fNIRS sources

%       and projected aux data

% ADD.aux_emb:   Temp. embedded accelerometer signals

% ADD.U,V:  sources in CCA space found by CCA

% ADD.gammaX,gammaAUX:  selected shrinkage parameters

% ADD.Av_red:   return (cca threshold) reduced mapping matrix Av

REG = round(100*rand(5,10));

ADD.Av = round(100*rand(5,10));