Merge branch 'master' of https://github.com/BUNPC/Homer3

% SYNTAX:

% data2 = hmrR_BandpassFilt( data, hpf, lpf )

% [data2, Aux2] = hmrR_BandpassFilt(data, Aux, hpf, lpf, turnon_dod, turnon_aux )

%

% UI NAME:

% Bandpass_Filter

% INPUT:

% data - SNIRF data type containing data time course to filter, time

%       vector, and channels.

% Aux - auxilliary data

% hpf - high pass filter frequency (Hz)

%       Typical value is 0 to 0.02.

% lpf - low pass filter frequency (Hz)

%       Typical value is 0.5 to 3.

% turnon_dod - apply filter on dod when set to 1

% turnon_aux - apply filter on aux when set to 1

%

% OUTPUT:

% data2 - SNIRF data type containing the filtered data time course, time

%        vector, and channels.

% Aux2 - filtered aux

%

% USAGE OPTIONS:

% Bandpass_Filter: dod = hmrR_BandpassFilt( dod, hpf, lpf )

% Bandpass_Filter: [dod, Aux2] = hmrR_BandpassFilt(dod, Aux, hpf, lpf, turnon_dod, turnon_aux )

%

% PARAMETERS:

% hpf: [0.010]

% lpf: [0.500]

%

% turnon_dod: 1

% turnon_aux: 1

function [data2, Aux2] = hmrR_BandpassFilt(data, Aux, hpf, lpf, turnon_dod, turnon_aux)

function [data2, ylpf] = hmrR_BandpassFilt( data, hpf, lpf )

if isa(data, 'DataClass')

if turnon_dod && turnon_aux

    foo{1} = data;

    foo{2} = Aux;

    data2 = DataClass().empty();

    Aux2 = AuxClass().empty();

elseif turnon_dod

    foo{1} = data;

    data2 = DataClass().empty();

elseif isa(data, 'AuxClass')

    data2 = AuxClass().empty();

    Aux2 = Aux;

elseif turnon_aux

    foo{1} = Aux;

    Aux2 = AuxClass().empty();

    data2 = data;

else

    data2 = data;

    Aux2 = Aux;

    return

end

for i = 1:size(foo,2)

    ylpf = [];

for ii=1:length(data)

    if isa(data, 'DataClass')

        data2(ii) = DataClass(data(ii));

    elseif isa(data, 'AuxClass')

        data2(ii) = AuxClass(data(ii));

    end

    for ii=1:length(foo{i})

        if isa(foo{i}, 'DataClass')

            data2(ii) = DataClass(foo{i}(ii));

            y = data2(ii).GetDataTimeSeries();

            fs = data2(ii).GetTime();

        elseif isa(foo{i}, 'AuxClass')

            Aux2(ii) = AuxClass(foo{i}(ii));

            y = Aux2(ii).GetDataTimeSeries();

            fs = Aux2(ii).GetTime();

        end

        % convert t to fs

        % assume fs is a time vector if length>1

        else

            y2 = ylpf;

        end

        if isa(foo{i}, 'DataClass')

            data2(ii).SetDataTimeSeries(y2);

        elseif isa(foo{i}, 'AuxClass')

            Aux2(ii).SetDataTimeSeries(y2);

        end

    end

end

 No newline at end of file

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

%

% UI NAME:

% User_Friendly_Name_For_hmrR_tCCA

% hmrR_tCCA

%

% DESCRIPTION:

% This script generates regressors using the regularized temporally embedded

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

%

% USAGE OPTIONS:

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

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

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

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

%

%

% PARAMETERS:

        end

        %% Select and prepare aux channels

        % lowpass filter aux signals from SNIRF aux argument 

        aux = hmrR_BandpassFilt(aux, 0, 0.5);

        % Extract variables from SNIRF aux

        kk = 1;

        for ii = 1:length(aux)

        % zscore AUX signals

        AUX = zscore(AUX);

        param.NumOfEmb = ceil(timelag*fq / tCCAparams(2));

        %% DO THE TCCA WORK

        filterFilename = sprintf('./tCCAfilter_%d.txt', iBlk);

        tCCAexists = isfile(sprintf('./tCCAfilter_%d.txt', iBlk));

        isTrainingRun = runIdxResting == runIdx;

        if ~tCCAexists && isTrainingRun

            dotCCA = 'train';

        elseif tCCAexists && ~isTrainingRun

            dotCCA = 'apply';

        else

            dotCCA = 'skip';

        end

        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

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

                %       number of embeddings

        param.NumOfEmb = ceil(timelag*fq / param.tau);

        if runIdx == runIdxResting

            %% if the tCCAfilter variable is empty, learn and put it out (this is the training/resting run)

                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);    %%%%%%% JD: rtcca generates error, "Undefined function or variable 'cshrink'". %%%%%%% 

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

                        % 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)

                    %% perform tCCA with shrinkage for all fNIRS channels

                [REG,  ADD] = rtcca_dummy(d_long,AUX,param,flags);    %%%%%%% JD: rtcca generates error, "Undefined function or variable 'cshrink'". %%%%%%% 

                %% save and output tCCA filter.

                    [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

                %%%%%%% JD: Undefined variable "ADD_trn" or class "ADD_trn.ccac". %%%%%%% 

                    if ctr < 1

                        % use only auxiliary tcca components that have correlation > ct

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

                    compindex=1:10;

                        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)

                    % if numel(ADD_trn.ccac) >= ctr

                    %     compindex = 1:ctr;

                    % else % not as many components available as regressors requested, provide all

                    %     compindex = 1:numel(ADD_trn.ccac);

                    % end                    

                    compindex=1:10;

                        if numel(ADD_trn.ccac) >= ctr

                            compindex = 1:ctr;

                        else % not as many components available as regressors requested, provide all

                            compindex = 1:numel(ADD_trn.ccac);

                        end

                    end

                    % return the reduced number of available regressors

                    Aaux{iBlk} = REG(:,compindex);

                    % return reduced mapping matrix Av, this is the tCCA filter

                tCCAfilter = ADD.Av(:,compindex);

                    tCCAfilter = ADD_trn.Av(:,compindex);

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

                rcMap{iBlk} = [];

                    rcMap{iBlk} = 'all';

                end

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

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

            print_filter(tCCAfilter);

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

        elseif exist(filterFilename,'file')

            %% if the tCCAfilter variable exists, apply the filtering and generate the tCCA regressors

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

            case 'apply'

                %% if the tCCAfilter variable exists, load it, apply the filtering and generate the tCCA regressors

                % Load the filter for the iBlk data block

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

                tCCAfilter = load(filterFilename,'-ascii');

            print_filter(tCCAfilter);

                % Temporal embedding and zscoring of auxiliary data

                aux_sigs = AUX;

                    aux_emb = [aux_emb aux];

                end

                %zscore

            % aux_emb = zscore(aux_emb);

            aux_emb = 1;     % JD: Dummy value 

                aux_emb = zscore(aux_emb);

                % Apply tCCA filter (generate regressors)

            Aaux{iBlk} = aux_emb * tCCAfilter;

                Aaux = aux_emb * tCCAfilter;

                % if the individual channel regressor flag was set, indicate by providing rc mapping

                if flagICRegressors

                    for cc = 1:size(d_long,2)

                        % 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

                rcMap{iBlk} = [];

                    rcMap{iBlk} = 'all';

                end

        else

            Aaux{iBlk} = [];

            rcMap{iBlk} = [];

            % tCCAfilter = []; if uncommented, retraining necessary whenever tcca is temporarily flagged out of processing stream

            case 'skip'

                Aaux = [];

                rcMap = [];

        end

    end

else

    Aaux = [];

    rcMap = [];

    % tCCAfilter = []; if uncommented, retraining necessary whenever tcca is temporarily flagged out of processing stream

end