Documentation update

FuncRegistry/UserFunctions/Archive/hmrS_SessAvgStd2_oldversion_nirs.m

@@ -5,13 +5,13 @@
 % Run_Average_Standard_Deviation
 %
 % DESCRIPTION:
-% Calculate avearge HRF standard deviation of all runs for one subject. 
+% Calculate average HRF standard deviation  of all runs for one subject. 
 %
 % INPUTS:
-% yAvgRuns:
-% ySum2Runs:
-% mlActRuns:
-% nTrialsRuns:
+% yAvgRuns: all runs for one subject, cell array of size [Number of runs X Number of data blocks]
+% ySum2Runs: ???
+% mlActRuns: logical cell array of size [Number of runs X Number of data blocks x Data channels*], determines active trials
+% nTrialsRuns: the number of trials averaged for each condition across all runs
 %
 % OUTPUTS:
 % yAvgStdOut: the standard deviation across runs

FuncRegistry/UserFunctions/hmrG_SubjAvg.m

@@ -8,12 +8,14 @@
 % Calculate the block average for all subjects, for all common stimuli accross subjects.
 %
 % INPUTS:
-% yAvgSubjs:
-% nTrialsSubjs:
+% yAvgSubjs:    cell array of length nSubj containing DataClass objects
+%               with the averaged results per subject
+% nTrialsSubjs: nSubj x nDataBlks cell array containing the number of 
+%               trials per subject per block
 %
 % OUTPUTS:
 % yAvgOut: the averaged results
-% nTrials: 
+% nTrials: amount of trials per block over which the average was computed
 %
 % USAGE OPTIONS:
 % Subj_Average_on_Concentration_Data: [dcAvg, nTrials] = hmrG_SubjAvg(dcAvgSubjs, nTrialsSubjs)
@@ -33,32 +35,42 @@
     subjCh = [];
     nStim = 0;
     grp1 = [];
-    nT = [];    
+    nT = []; % number of trials
 
     for iSubj = 1:nSubj
 
         yAvgOut(iBlk) = DataClass();        
 
+        % get current subject input for the current data block
         yAvg      = yAvgSubjs{iSubj}(iBlk).GetDataTimeSeries('reshape');
-        if isempty(yAvg)
+        if isempty(yAvg) % if no input, create flag and skip subject
             err(iBlk, iSubj) = -1;
             continue;
         end
 
         tHRF      = yAvgSubjs{iSubj}(iBlk).GetTime();
-        nT        = nTrialsSubjs{iSubj}{iBlk};
-        datatype  = yAvgSubjs{iSubj}(iBlk).GetDataTypeLabel();
+        nT        = nTrialsSubjs{iSubj}{iBlk};  % get number of trials in this block
+        datatype  = yAvgSubjs{iSubj}(iBlk).GetDataTypeLabel();  % check if Hb or OD data
+
+        % check if the trial contains Hb or OD data
         if strncmp(datatype{1}, 'HRF Hb', length('HRF Hb'))
+            % get source-detector pairs (renders nChans x 2 matrix; 
+            % 1st column contains source indeces, 2nd column contains dectector indices)
             ml    = yAvgSubjs{iSubj}(iBlk).GetMeasListSrcDetPairs('reshape');
         elseif strcmp(datatype{1}, 'HRF dOD')
             ml    = yAvgSubjs{iSubj}(iBlk).GetMeasList('reshape');
         end
 
-        nCond = size(nT,2);
+        nCond = size(nT,2);  % number of conditions to average over
         yAvgOut(iBlk).SetTime(tHRF);
 
+
         if strncmp(datatype{1}, 'HRF Hb', length('HRF Hb'))            
             if iSubj==1
+                % initialize data array for group 1 with HbO, HbR and 
+                % HbT data per condition
+                % shape: time x 3 x nChan x nCond (second dim includes the
+                % three Hb data typeS)                
                 grp1 = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), nCond);
             end
 
@@ -68,7 +80,7 @@
             end
 
             for iC = 1:nCond
-                if sum(nT(:,iC))==0
+                if sum(nT(:,iC))==0  % if no trials for iC condition, skip
                     continue;
                 end
 
@@ -103,15 +115,20 @@
                         end
                     end
                 end
+                % keep count of number of conditions (i.e. number of times
+                % the averages were added) to compute mean afterwards 
                 subjCh(lstPass,iC) = subjCh(lstPass,iC) + 1; %#ok<*AGROW>
             end
 
             yAvg = [];
             if iSubj == length(yAvgSubjs)
                 for iC = 1:size(grp1,4)
                     for iCh = 1:size(grp1,3)
+                        % compute group mean
                         yAvg(:,:,iCh,iC) = grp1(:,:,iCh,iC) / subjCh(iCh,iC);
                         if iSubj == nSubj
+                            % initialize data array for each HbO type and
+                            % condition, per channel
                             yAvgOut(iBlk).AddChannelHbO(ml(iCh,1), ml(iCh,2), iC);
                             yAvgOut(iBlk).AddChannelHbR(ml(iCh,1), ml(iCh,2), iC);
                             yAvgOut(iBlk).AddChannelHbT(ml(iCh,1), ml(iCh,2), iC);                            
@@ -177,6 +194,7 @@
                         yAvgOut(iBlk).AddChannelDod(ml(iCh,1), ml(iCh,2), ml(iCh,4), iC);
                     end
                 end
+                % save the average over subjects in output data array
                 yAvgOut(iBlk).AppendDataTimeSeries(yAvg);
             end
 

FuncRegistry/UserFunctions/hmrG_SubjAvgStd.m

@@ -8,7 +8,7 @@
 % Calculate avearge HRF standard deviation and standard error of all subjects in a group.
 %
 % INPUTS:
-% ySubjAvg:
+% ySubjAvg: NIRS data cell array of cize [Number of subjects; Number of data blocks ]
 %
 % OUTPUTS:
 % yAvgStdOut: the standard deviation across subjects.

FuncRegistry/UserFunctions/hmrR_BandpassFilt.m

@@ -6,6 +6,8 @@
 %
 % DESCRIPTION:
 % Perform a bandpass filter on time course data.
+% The type of filter used is a Butterworth filter of order 3 for the
+% lowpass case and 5 for the highpass. 
 %
 % INPUT:
 % data - SNIRF data type containing data time course to filter, time vector, and channels.
@@ -21,19 +23,22 @@
 % Bandpass_Filter_OpticalDensity: dod = hmrR_BandpassFilt(dod, hpf, lpf)
 % Bandpass_Filter_Auxiliary: aux = hmrR_BandpassFilt(aux, hpf, lpf)
 %
-% PARAMETERS:
+% DEFAULT PARAMETERS:
 % hpf: [0.000]
 % lpf: [0.500]
 %
 % PREREQUISITES:
 % Intensity_to_Delta_OD: dod = hmrR_Intensity2OD( intensity )
 
 function [data2, ylpf] = hmrR_BandpassFilt( data, hpf, lpf )
+
+% instantiate output data object where computation will be performed
 if isa(data, 'DataClass')
     data2 = DataClass().empty();
 elseif isa(data, 'AuxClass')
     data2 = AuxClass().empty();
 end
+
 ylpf = [];
 for ii=1:length(data)
     if isa(data, 'DataClass')
@@ -70,22 +75,23 @@
 
     % low pass filter
     lpf_norm = lpf / (fs / 2);
-    if lpf_norm > 0  % No lowpass if filter is 
+    if lpf_norm > 0  % No lowpass if cutoff is <= 0
         FilterOrder = 3;
         [z, p, k] = butter(FilterOrder, lpf_norm, 'low');
-        [sos, g] = zp2sos(z, p, k);
+        [sos, g] = zp2sos(z, p, k); % zero-pole-gain to second-order sections model conversion
         y2 = filtfilt(sos, g, double(y)); 
     end
 
     % high pass filter
     hpf_norm = hpf / (fs / 2);
-    if hpf_norm > 0
+    if hpf_norm > 0 % No highpass if cutoff is <= 0
         FilterOrder = 5;
         [z, p, k] = butter(FilterOrder, hpf_norm, 'high');
-        [sos, g] = zp2sos(z, p, k);
+        [sos, g] = zp2sos(z, p, k); % zero-pole-gain to second-order sections model conversion
         y2 = filtfilt(sos, g, y2);
     end
 
+    % store filtered time series in output data object
     data2(ii).SetDataTimeSeries(y2);
 
 end

FuncRegistry/UserFunctions/hmrR_BlockAvg.m

@@ -12,22 +12,22 @@
 % the data range, then the trial is excluded from the average.
 %
 % INPUT:
-% data: SNIRF.data container with the delta OD or delat concentration data
+% data: SNIRF.data container with the delta OD or delta concentration data
 % stim: SNIRF.stim container with the stimulus condition data
 % trange: defines the range for the block average [tPre tPost]
 %
 % OUTPUT:
 % data_avg: SNIRF.data container with the averaged results
 % data_std: SNIRF.data container with the standard deviation across trials
 % nTrials: the number of trials averaged for each condition
-% data_sum2: SNIRF.data container ...
+% data_sum2: SNIRF.data container of the squared sum of the trials
 % yTrials: a structure containing the individual trial responses
 %
 % USAGE OPTIONS:
 % Block_Average_on_Concentration_Data: [dcAvg, dcAvgStd, nTrials, dcSum2] = hmrR_BlockAvg( dc, stim, trange )
 % Block_Average_on_Delta_OD_Data: [dodAvg, dodAvgStd, nTrials, dodSum2] = hmrR_BlockAvg( dod, stim, trange )
 %
-% PARAMETERS:
+% DEFAULT PARAMETERS:
 % trange: [-2.0, 20.0]
 %
 % PREREQUISITES:
@@ -64,27 +64,32 @@
     y = data(kk).GetDataTimeSeries('reshape');    % Get the data vector 
     t = data(kk).GetTime();    % Get the time vector 
     dt = t(2)-t(1);
+    % convert pre and post stimulus time to samples
     nPre = round(trange(1)/dt);
     nPost = round(trange(2)/dt);
     nTpts = size(y,1);
-    tHRF = nPre*dt:dt:nPost*dt;
+    tHRF = nPre*dt:dt:nPost*dt; % time axis for averaged response (HRF)
     if strncmp(datatype{1}, 'Hb', 2)
         ml = data(kk).GetMeasListSrcDetPairs('reshape');
+        % initialize array containing all Hb signals (HbO, HbR and Hbt)
+        % shape: time x nHb x nTrials x nConditions
         yblk = zeros(nPost-nPre+1,size(y,2),size(y,3),size(s,2));
     elseif strcmp(datatype{1}, 'dOD')
         ml = data(kk).GetMeasList('reshape');
+        % initialize array containing all OD signals 
+        % shape: time x nOD x nTrials x nConditions 
         yblk = zeros(nPost-nPre+1,size(y,2),size(s,2));
     else
         return;
     end
-    
-    for iC = 1:size(s,2)
-        lstS = find(s(:,iC)==1);
-        nBlk = 0;
-        for iT = 1:length(lstS)
+
+    for iC = 1:size(s,2) % iterate over experimental conditions
+        lstS = find(s(:,iC)==1); % get stimuli onsets
+        nBlk = 0; % valid trials counter
+        for iT = 1:length(lstS) % iterate over block trials
             if (lstS(iT)+nPre)>=1 && (lstS(iT)+nPost)<=nTpts
                 if strncmp(datatype{1}, 'Hb', 2)
-                    nBlk = nBlk + 1;
+                    nBlk = nBlk + 1; 
                     yblk(:,:,:,nBlk) = y(lstS(iT)+[nPre:nPost],:,:); %changed from yblk(:,:,:,end+1)
                 elseif strcmp(datatype{1}, 'dOD')
                     nBlk = nBlk + 1;
@@ -103,12 +108,16 @@
 
             % Loop over all channels
             for ii=1:size(yavg,3)
-                foom = ones(size(yavg,1),1)*mean(yavg(1:-nPre,:,ii,iC),1);
+                % set the baseline of the average to zero by subtracting
+                % the mean of the average for t<0.
+                foom = ones(size(yavg,1),1)*mean(yavg(1:-nPre,:,ii,iC),1); % mean of the average for t<0.
                 yavg(:,:,ii,iC) = yavg(:,:,ii,iC) - foom;
 
+                % similarly, per block
                 for iBlk = 1:nBlk
                     yTrials(iC).yblk(:,:,ii,iBlk) = yTrials(iC).yblk(:,:,ii,iBlk) - foom;
                 end
+                % compute squared sum of the trials
                 ysum2(:,:,ii,iC) = sum( yTrials(iC).yblk(:,:,ii,1:nBlk).^2 ,4);
 
                 % Snirf stuff: set channel descriptors
@@ -122,7 +131,7 @@
                 data_std(kk).AddChannelHbR(ml(ii,1), ml(ii,2), iC);
                 data_std(kk).AddChannelHbT(ml(ii,1), ml(ii,2), iC);
 
-                % 
+                % squared sum of the tirals
                 data_sum2(kk).AddChannelHbO(ml(ii,1), ml(ii,2), iC);
                 data_sum2(kk).AddChannelHbR(ml(ii,1), ml(ii,2), iC);
                 data_sum2(kk).AddChannelHbT(ml(ii,1), ml(ii,2), iC);
@@ -140,12 +149,16 @@
 
             % Loop over all wavelengths
             for ii=1:size(yavg,2)
-                foom = ones(size(yavg,1),1)*mean(yavg(1:-nPre,ii,iC),1);
+                % set the baseline of the average to zero by subtracting
+                % the mean of the average for t<0.
+                foom = ones(size(yavg,1),1)*mean(yavg(1:-nPre,ii,iC),1); % mean of the average for t<0.
                 yavg(:,ii,iC) = yavg(:,ii,iC) - foom;
 
+                % similarly, per block
                 for iBlk = 1:nBlk
                     yTrials(iC).yblk(:,ii,iBlk) = yTrials(iC).yblk(:,ii,iBlk) - foom;
                 end
+                % compute squared sum of the trials
                 ysum2(:,ii,iC) = sum( yTrials(iC).yblk(:,ii,1:nBlk).^2 ,3);
 
                 % Snirf stuff: set channel descriptors

FuncRegistry/UserFunctions/hmrR_GLM.m

@@ -1,6 +1,7 @@
 % SYNTAX:
 % [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)
+% hmrR_GLM(data_y, stim, probe, mlActAuto, Aaux, tIncAuto, rcMap, trange, glmSolveMethod, idxBasis,...
+% paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, c_vector)
 %
 % UI NAME:
 % GLM_HRF_Drift_SS
@@ -17,11 +18,12 @@
 % data - this is the concentration data with dimensions #time points x [HbO/HbR/HbT] x #channels
 % stim - stimulation vector (# time points x #conditions)=1 at stim onset otherwise =0
 % probe - source detector stucture (units should be consistent with rhoSD_ssThresh)
-% mlActAuto -
-% Aaux - A matrix of auxilliary regressors (#time points x #Aux channels)
+% mlActAuto - cell array of vectors, one for each time base in data, specifying 
+%             active/inactive channels with 1 meaning active, 0 meaning inactive
+% Aaux - a matrix of auxilliary regressors (#time points x #Aux channels)
 % tIncAuto - a vector (#time points x 1) indicating which data time points
 %            are motion (=0) or not (=1)
-% rcMap - An array of cells (1 x #fNIRS channels) containing aux regressor
+% rcMap - an array of cells (1 x #fNIRS channels) containing aux regressor
 %           indices for individual regressor-channel mapping. Currently
 %           only relevant when flagNuisanceRMethod = 3 (tCCA regressors).
 % trange - defines the range for the block average [tPre tPost dt]. If dt
@@ -45,7 +47,7 @@
 % 			         (t/(p*q))^p * exp(p-t/q)
 %                Defaults: p=8.6 q=0.547
 %                The peak is at time p*q.  The FWHM is about 2.3*sqrt(p)*q.
-% paramsBasis - Parameters for the basis function (chosen via idxBasis)
+% paramsBasis - parameters for the basis function (chosen via idxBasis)
 %               if idxBasis=1 [stdev step ~ ~ ~ ~] where stdev is the width of the
 %                  gaussian and step is the temporal spacing between
 %                  consecutive gaussians
@@ -69,8 +71,8 @@
 %            2. if performed with average of all short separation channels.
 %            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
+% 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]
 %
@@ -84,17 +86,19 @@
 % yresid - the residual between the data y and the GLM fit
 % ysum2 - an intermediate matrix for calculating stdev across runs
 % beta - the coefficients of the temporal basis function fit for the HRF
-%           (#coefficients x HbX x #Channels x #conditions)
+%        (#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 p values for GLM and the corresponding beta_label and ml
 %     (#Betas x #Channels x HbX) for conditions
 %     (#Channels x HbX) for contrasts 
 %
 % USAGE OPTIONS:
-% 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)
+% 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:
+% DEFAULT PARAMETERS:
 % trange: [-2.0, 20.0]
 % glmSolveMethod: 1
 % idxBasis: 1
@@ -124,6 +128,7 @@
 
 % Check input args
 if isempty(tIncAuto)
+    % assume there are no motion artefacts, include all data points
     tIncAuto = cell(length(data_y),1);
 end
 if isempty(mlActAuto)