Merge on behalf of Meryem Ayse Yucel and Yuanyuan Gao of new user functions to...

% SYNTAX:

% [yavg, yavgstd, tHRF, nTrials, ynew, yresid, ysum2, beta, R] = hmrR_GLM(data, stim, probe, mlActAuto, Aaux, tIncAuto, trange, rcMap, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder)

% [yavg, yavgstd, tHRF, nTrials, ynew, yresid, ysum2, beta, R, hmrstats] = hmrR_GLM(data, stim, probe, mlActAuto, Aaux, tIncAuto, trange, rcMap, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, c_vector)

%

% UI NAME:

% GLM_HRF_Drift_SS

%            3. uses tCCA regressors for nuisance regression, in Aaux,

%            mapped by rcMap, provided by hmr_tCCA()

% driftOrder - Polynomial drift correction of this order

%

% c_vector - Contrast vector, has values 1, -1 or 0. E.g. to contrast cond

%           2 to cond 3 in an experimental paradigm with four conditions, c_vector is

%           [0 1 -1 0]

%

% OUTPUTS:

% yavg - the averaged results

%           (#coefficients x HbX x #Channels x #conditions)

% R - the correlation coefficient of the GLM fit to the data

%     (#Channels x HbX)

%

% hmrstats - outputs t and pvalues for GLM

%

% USAGE OPTIONS:

% GLM_HRF_Drift_SS_Concentration: [dcAvg, dcAvgStd, nTrials, dcNew, dcResid, dcSum2, beta, R] = hmrR_GLM(dc, stim, probe, mlActAuto, Aaux, tIncAuto, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder)

% GLM_HRF_Drift_SS_Concentration: [dcAvg, dcAvgStd, nTrials, dcNew, dcResid, dcSum2, beta, R, hmrstats] = hmrR_GLM(dc, stim, probe, mlActAuto, Aaux, tIncAuto, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, c_vector)

%

%

% PARAMETERS:

% rhoSD_ssThresh: 15.0

% flagNuisanceRMethod: 1

% driftOrder: 3

% c_vector: 0

%

%

function [data_yavg, data_yavgstd, nTrials, data_ynew, data_yresid, data_ysum2, beta_blks, yR_blks] = ...

    hmrR_GLM(data_y, stim, probe, mlActAuto, Aaux, tIncAuto, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder)

function [data_yavg, data_yavgstd, nTrials, data_ynew, data_yresid, data_ysum2, beta_blks, yR_blks, hmrstats] = ...

    hmrR_GLM(data_y, stim, probe, mlActAuto, Aaux, tIncAuto, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, c_vector)

% Init output

data_yavg     = DataClass().empty();

beta_blks     = cell(length(data_y),1);

yR_blks       = cell(length(data_y),1);

beta_label = [];

hmrstats = [];

% Check input args

if isempty(tIncAuto)

                    yest(:,lstML,conc) = At * foo(:,lstML,conc);

                    yvar(1,lstML,conc) = sum((squeeze(y(:,conc,lstML))-yest(:,lstML,conc)).^2)./(size(y,1)-1); % check this against eq(53) in Ye2009

                    for iCh = 1:length(lstML)

                        % GLM stats for each condition

                        bvar(:,lstML(iCh),conc) = yvar(1,lstML(iCh),conc) * pAinvAinvD;

                        tval(:,lstML(iCh),conc) =  foo(:,lstML(iCh),conc)./sqrt(bvar(:,lstML(iCh),conc));

                        pval(:,lstML(iCh),conc) = 1-tcdf(abs(tval(:,lstML(iCh),conc)),(size(y,1)-1));

                        %

                        % GLM stats for contrast between conditions, given a c_vector exists

                        if nCond > 1

                            if (sum(abs(c_vector)) ~= 0) && (size(c_vector,2) == nCond)

                                if ~exist('cv_extended') == 1

                                    cv_dummy = [];

                                    for m = 1:nCond

                                        cv_dummy = [cv_dummy ones(1,nB)*c_vector(m)];

                                    end

                                    cv_extended = [cv_dummy zeros(1,size(beta_label,2)-size(cv_dummy,2))];

                                end

                                tval_contrast(:,lstML(iCh),conc) = cv_extended * foo(:,lstML(iCh),conc)./sqrt(cv_extended * (pinvA*pinvA') * yvar(:,lstML(iCh),conc) * cv_extended');

                                pval_contrast(:,lstML(iCh),conc) = 1-tcdf(abs(tval_contrast(:,lstML(iCh),conc)),(size(y,1)-1));

                            end

                        end

                        %

                        for iCond=1:nCond

                            if size(tbasis,3)==1

                                yavgstd(:,lstML(iCh),conc,iCond) = diag(tbasis*diag(bvar([1:nB]+(iCond-1)*nB,lstML(iCh),conc))*tbasis').^0.5;

    yR_blks{iBlk}   = yR;

    % stats struct

    if glmSolveMethod == 1 %  for OLS only now

    if glmSolveMethod == 1 %  for OLS only for now

        % GLM stats for each condition 

        hmrstats.beta_label = beta_label;

        hmrstats.tval = tval;

        hmrstats.pval = pval;

        hmrstats.ml = ml;

        % GLM stats for contrast between conditions, if c_vector exists

        if (sum(abs(c_vector)) ~= 0) && (size(c_vector,2) == nCond) && nCond>1

            hmrstats.tval_contrast = tval_contrast;

            hmrstats.pval_contrast = pval_contrast;

            hmrstats.contrast = c_vector;

        end

    end

end

 No newline at end of file

for ii=1:length(intensity)

    dod(ii) = DataClass();

    d = intensity(ii).GetDataTimeSeries();

    % Optional (user prompt): Adding dc offset if intensity (d) has negative values

    if ~isempty(d(d<=0))

        quest = {'Intensity signal has negative values. If you would like to add a dc offset, please click YES. If you would like to proceed with negative values, hit CANCEL.'};

        dlgtitle = 'Warning';

        btn1 = 'YES';

        btn2 = 'CANCEL';

        detbtn = btn1;

        answer = questdlg(quest,dlgtitle,btn1,btn2,detbtn)

        switch answer

            case 'YES'

                for j = 1:size(d,2)

                    foo = d(:,j);

                    if ~isempty(foo<0)

                        d(:,j) = foo + abs(min(foo));

                        foo = d(:,j);

                    end

                    if ~isempty(foo == 0)

                        d(:,j) = foo + min(foo(foo > 0));

                    end

                end

        end

    end

    dm = mean(abs(d),1);

    nTpts = size(d,1);

    dod(ii).SetTime(intensity(ii).GetTime());

% SYNTAX:

% intensity = hmrR_PreprocessIntensity_MedianFilter( intensity )

%

% UI NAME:

% hmrR_PreprocessIntensity_MedianFilter

%

% DESCRIPTION:

% Applies a median filter to data to remove huge spikes.

%

% INPUT:

% intensity - SNIRF data type where the d matrix is intensity

%

% OUTPUT:

% intensity - SNIRF data type where the d matrix is intensity

%

% USAGE OPTIONS:

% Intensity_to_Intensity: d = hmrR_PreprocessIntensity_MedianFilter(data)

function d = hmrR_PreprocessIntensity_MedianFilter( intensity)

for ii=1:length(intensity)

    d = intensity(ii).GetDataTimeSeries();

    for j = 1:size(d,2)

        foo = d(:,j);

        xdata = (1:length( foo))';

        d(:,j) = interp1(xdata(~isnan( foo)), foo(~isnan( foo)),xdata,'spline');

        new_signal = zeros(size(foo));

        new_signal(1) = foo(1);

        new_signal(end) = foo(end);

        for tp = 2:length(foo)-1

            values = foo([tp-1,tp,tp+1]);

            median_value = median(values);

            new_signal(tp) = median_value;

        end

        d(:,j) = new_signal;

    end

    intensity(ii).SetDataTimeSeries(d);

end

% SYNTAX:

% intensity = hmrR_PreprocessIntensity_NAN( intensity )

%

% UI NAME:

% hmrR_PreprocessIntensity_NAN

%

% DESCRIPTION:

% replace NAN by spline interpolation of the nonnan values

%

% INPUT:

% intensity - SNIRF data type where the d matrix is intensity

%

% OUTPUT:

% intensity - SNIRF data type where the d matrix is intensity

%

% USAGE OPTIONS:

% Intensity_to_Intensity: d = hmrR_PreprocessIntensity_NAN(data)

function d = hmrR_PreprocessIntensity_NAN( intensity)

for ii=1:length(intensity)

    d = intensity(ii).GetDataTimeSeries();

    if ~isempty(find(isnan(d)))

        for j = 1:size(d,2)

            foo = d(:,j);

            if ~isempty(find(isnan(foo)))

                xdata = (1:length( foo))';

                d(:,j) = interp1(xdata(~isnan( foo)), foo(~isnan( foo)),xdata,'spline');

            end

        end

    end

    intensity(ii).SetDataTimeSeries(d);

end

% SYNTAX:

% intensity = hmrR_PreprocessIntensity_Negative( intensity )

%

% UI NAME:

% hmrR_PreprocessIntensity_Negative

%

% DESCRIPTION:

% Fix negative values in intensity by either adding an offset to make all

% numbers positive (+eps)(option 1) or just set values <=0 to +eps (option

% 2)

%

% INPUT:

% intensity - SNIRF data type where the d matrix is intensity

%

% OUTPUT:

% intensity - SNIRF data type where the d matrix is intensity

%

% USAGE OPTIONS:

% Intensity_to_Intensity: d = hmrR_PreprocessIntensity_Negative(data)

%

function d = hmrR_PreprocessIntensity_Negative( intensity )

for ii=1:length(intensity)

    d = intensity(ii).GetDataTimeSeries();

    if ~isempty(d(d<=0))

        quest = {'Intensity signal has negative values.'};

        dlgtitle = 'Warning';

        btn1 = 'OPTION1: Add a dc offset';

        btn2 = 'OPTION2: Set values <=0 to eps';

        btn3 = 'CANCEL';

        detbtn = btn3;

        answer = questdlg(quest,dlgtitle,btn1,btn2,btn3, detbtn);

        switch answer  

            case 'OPTION1: Add a dc offset'

                for j = 1:size(d,2)

                    foo = d(:,j);

                    if ~isempty(find(foo<0))

                        d(:,j) = foo + abs(min(foo)) + eps;

                    end

                end

            case 'OPTION2: Set values <=0 to eps'

                for j = 1:size(d,2)

                    foo = d(:,j);

                    if ~isempty(find(foo<0))

                        d(find(foo<=0),j) = eps;

                    end

                end

        end

    end

    intensity(ii).SetDataTimeSeries(d);

end