New function hmrR_TaskNorm

FuncRegistry/UserFunctions/hmrR_TaskNorm.m

@@ -0,0 +1,273 @@
+% SYNTAX:
+% [dcAvg, dcAvgStd, nTrials] = hmrR_TaskNorm(dc, stim, deltat, show )
+%
+% UI NAME:
+% Task_Normalization
+%
+% DESCRIPTION:
+% Calculate the average of the HRF over the duration of the task.
+% copyright  � 2023, Medea.
+% Authors: Anna Falivene, Charlotte Johnson
+
+%
+% INPUT:
+% dc: SNIRF.data container with the concentration data
+% stim: SNIRF.stim container with the stimulus condition data
+% deltat: defines the range before the onset and after the duration
+% [tPre_onset tPost_duration]. Both positive (or null) values
+% show: parameter to able/disable figure(). if show=0 no figure will
+% appear, if =1 figures with the HBO and HBR concentration for each
+% channel will appear
+%
+% OUTPUT:
+% dcAvg: averaged results
+% dcAvgStd: standard deviation across trials
+% nTrials: the number of trials averaged
+%
+% USAGE OPTIONS:
+% Task_Normalization_on_Concentration_Data: [dcAvg, dcAvgStd, nTrials] = hmrR_TaskNorm(dc, stim, deltat, show )
+%
+% PARAMETERS:
+% deltat: [0, 0]
+% show: 1
+%
+% PREREQUISITES:
+% Delta_OD_to_Conc: dc = hmrR_OD2Conc( dod, probe, ppf )
+%
+function [dcAvg, dcAvgStd, nTrials] = hmrR_TaskNorm(dc, stim, deltat, show )
+
+snirf = SnirfClass(dc, stim); %read the concentration data
+
+delta1=deltat(1); %baseline duration
+delta2=deltat(2); %post task duration
+dcAvg= DataClass();
+dcAvgStd = DataClass();
+d = dc.GetDataTimeSeries();
+
+%to compute the mean duration of the stimuli/tasks considered
+m_durations1=[];
+for ic=1:size(stim,2)
+    if size(stim(1, ic).data) ~=0
+
+        s1= stim(:,ic).GetData();
+        m_durations1=[m_durations1;mean(s1(:,2),1)];
+    end
+end
+
+for ic=1:size(stim,2)
+    if size(stim(1, ic).data) ~=0
+
+        s1= stim(:,ic).GetData();
+    else
+        continue
+    end
+    t = snirf.GetTimeCombined();
+    Fs=1/(t(2)-t(1)); %sampling rate
+    nTrials=size(s1,1);
+    namestr=[];
+
+
+    for i=1:nTrials
+        my_field = string(['Repetition',num2str(i)]);
+         namestr=[namestr;my_field];
+        str.(my_field) = [];
+
+%data segmentation of each repetition
+        T1=find(t>=s1(i,1),1,'first');%finding onset
+        T2=find(t<=s1(i,1)+s1(i,2),1,'last');%finding end of trial (onset+dur)
+        newtime=t(T1:T2);
+        newdata=d(T1:T2,:); 
+        rsdata=(round(mean(m_durations1)))*50;%number of samples you want to have when interpolating/normalizing 
+        normtime=linspace(newtime(1),newtime(end),rsdata);
+        %data interpolation for the task execution phase of the trial
+        norm_data=interp1(newtime, newdata, normtime, 'spline');
+%data interpolation for the baseline and post-task phases of the trial
+        if delta1~=0 | delta2~=0
+
+            %data segmentation of baseline and post task phases
+            T1d=T1-round(delta1*Fs);
+            T2d=T2+round(delta2*Fs);
+
+            rs1=(delta1*50);%defines numer of samples of delta1
+            rs2=(delta2*50);%defines numer of samples of delta2
+            n1=linspace(t(T1d),t(T1),rs1);%=Generate linearly spaced vector.
+            n2=linspace(t(T2),t(T2d),rs2);% =Generate linearly spaced vector.
+            if delta1~=0
+                d1=interp1(t(T1d:T1), d(T1d:T1,:), n1, 'spline');% interpolation of data in delta1
+            else
+                d1=[];
+            end
+            if delta2~=0
+                d2=interp1(t(T2:T2d), d(T2:T2d,:), n2, 'spline');% interpolation of data in delta2
+            else
+                d2=[];
+            end
+
+            norm_data=[ d1(1:end-1,:); norm_data; d2(2:end,:)];
+
+        end
+
+        str.(namestr(i,:))=norm_data;%saves new data in a structure
+    end
+
+    meanvalues=[]; stdvalues=[]; %forming a matrix where the MEAN and STD values are stored
+
+    ml = dc.GetMeasListSrcDetPairs('reshape'); %opening the measlist useful to calculate mean and stds 
+
+    %compute mean and std over the trials
+    for k=1:(3*size(ml,1))
+        ALLM=[];
+
+        for nt=1:nTrials
+            ALLM=[ALLM str.(namestr(nt,:))(:,k)];
+        end
+        meanvalues=[meanvalues mean(ALLM,2)];
+        stdvalues=[stdvalues std(ALLM,1,2)];
+
+    end
+
+    % '%cycle' vector
+    if delta1~=0 & delta2~=0
+        x1=linspace(0,100,rsdata+rs1+rs2-2);%=Generate linearly spaced vector from 0-100%.
+        disp (['onset: ', num2str(x1(rs1-1)), '%', ' - ', 'end of stimulus: ' ,num2str(x1(rs1-1+rsdata)), '%'])
+    elseif delta1~=0 & delta2==0 %only baseline
+        x1=linspace(0,100,rsdata+rs1-1);
+        disp (['onset: ', num2str(x1(rs1-1)), '%'])
+    elseif delta1==0 & delta2~=0 %no baseline but only resting period after the stim
+        x1=linspace(0,100,rsdata+rs2-1);
+        disp (['end of stimulus: ' ,num2str(x1(rs1-1+rsdata)), '%'])
+    else
+        x1=linspace(0,100,rsdata); %only exercise task
+    end
+
+    avghrf=[]; figs=[];
+    tot_avg=[];
+
+    %based on BlockAVG, Computing mean and std per channel in the different tasks.
+    for ll=1:3:size(meanvalues,2)-2
+
+        yHBO=meanvalues(:,ll); %HbO
+        yHBR=meanvalues(:,ll+1); %hbR
+        yHBT=meanvalues(:,ll+2); %hbt
+        sdHBO=stdvalues(:,ll);
+        sdHBR=stdvalues(:,ll+1);
+        sdHBT=stdvalues(:,ll+2);
+
+
+        ALLHBO=[];ALLHBR=[];%creating a matrix for All HbO and HbR data
+
+        for nt=1:nTrials
+%data concentration of each trial
+            ALLHBO=[ALLHBO str.(namestr(nt,:))(:,ll)];
+            ALLHBR=[ALLHBR str.(namestr(nt,:))(:,ll+1)];
+        end
+
+%baseline correction (as in Homer BlockAvg function) can be performed only if delta1 is different from 0
+        if delta1~=0
+            foomHBO = ones(size(yHBO,1),1)*mean(yHBO(1:rs1-1),1);
+            foomHBR = ones(size(yHBR,1),1)*mean(yHBR(1:rs1-1),1);
+            foomHBT = ones(size(yHBT,1),1)*mean(yHBT(1:rs1-1),1);
+            avghrf(:,ll)=yHBO-foomHBO;
+            avghrf(:,ll+1)=yHBR-foomHBR;
+            avghrf(:,ll+2)=yHBT-foomHBT;
+
+        else
+            %no correction
+            avghrf(:,ll)=yHBO;
+            avghrf(:,ll+1)=yHBR;
+            avghrf(:,ll+2)=yHBT;
+        end
+
+%%% Plot section 
+        f=figure(Visible=show);
+        ax1=subplot('Position',[0.35 0.58 0.33 0.34 ]);
+        txt=['Source ',num2str(dc.measurementList(1,ll).sourceIndex),' detector ', num2str(dc.measurementList(1,ll).detectorIndex)];
+
+        patch([x1 fliplr(x1)], [avghrf(:,ll)-sdHBO ;flipud(avghrf(:,ll)+sdHBO)], [0.6350 0.0780 0.1840] , 'FaceAlpha', 0.5 )
+        hold on
+        plot(x1,avghrf(:,ll),'r', 'LineWidth',2)
+
+        patch([x1 fliplr(x1)], [avghrf(:,ll+1)-sdHBR ;flipud(avghrf(:,ll+1)+sdHBR)], [0.3 0.75 0.9], 'FaceAlpha', 0.5)
+        plot(x1,avghrf(:,ll+1),'Color', '#0072BD', 'LineWidth',2)
+
+        hold off
+        lgd=legend('HbO','', 'HbR','');
+        title(txt )
+
+        ax2=subplot(223);
+        plot(x1,ALLHBO, 'LineWidth',1)
+        legend(namestr(:,:))
+        title('HbO_2 in each repetition')
+
+
+        ax3=subplot(224);
+        plot(x1,ALLHBR, 'LineWidth',1)
+        legend(namestr(:,:))
+        title('HbR in each repetition')
+
+        axx=[ax1 ax2 ax3];
+        grid (axx, 'on' )
+        ax1.YGrid='off';
+        xlabel(axx,'% of trial')
+        % ylabel(axx, '[M]')
+
+        ax1.GridAlpha=0.5;
+        xtickangle([ax1 ax2 ax3],0)
+
+        if delta1~=0 & delta2~=0
+            xticks([ax1 ax2 ax3],[0 x1(rs1-1) x1(rs1-1+rsdata) 100])
+            xticklabels([ax1 ax2 ax3],{'0','onset', 'end of the stimulus', '100'})
+
+
+        elseif delta1~=0 & delta2==0 %only baseline
+            xticks([ax1 ax2 ax3],[0 x1(rs1-1) 100])
+            xticklabels([ax1 ax2 ax3],{'0','onset', '100'})
+
+        elseif delta1==0 & delta2~=0 %no baseline but only resting period after the stim
+            xticks([ax1 ax2 ax3],[0  x1(end-(rs2-1)) 100])
+            xticklabels([ax1 ax2 ax3],{'0', 'end of the stimulus', '100'})
+
+        else %only exercises
+            xticks([ax1 ax2 ax3],[0  100])
+            xticklabels([ax1 ax2 ax3],{'0','100'})
+
+        end
+
+        hlink = linkprop([ax1 ax2 ax3],{'YGrid','GridAlpha'});
+
+
+        figs=[figs; f];
+    end
+    figs=(find(figs~=0));
+
+%saving figures
+    global maingui
+    folder=uigetdir();%decide where figures are stored.
+    namefig=maingui.dataTree.currElem.name; %can be modified according to how the .snirf file were named
+    namefig=namefig(1:end-6);
+    namefig=[namefig, '.fig'];
+    HRF=fullfile(folder,namefig);
+
+    savefig(figs,HRF,'compact');
+
+%add the data into the dcAvg and dcAvgStd Dataclasses
+
+    for nn=1:size(ml,1)
+        dcAvg.AddChannelHbO(ml(nn,1), ml(nn,2), ic);
+        dcAvg.AddChannelHbR(ml(nn,1), ml(nn,2), ic);
+        dcAvg.AddChannelHbT(ml(nn,1), ml(nn,2), ic);
+        dcAvgStd.AddChannelHbO(ml(nn,1), ml(nn,2), ic);
+        dcAvgStd.AddChannelHbR(ml(nn,1), ml(nn,2), ic);
+        dcAvgStd.AddChannelHbT(ml(nn,1), ml(nn,2), ic);
+
+    end
+
+    dcAvg.AppendDataTimeSeries(avghrf);
+    dcAvgStd.AppendDataTimeSeries(stdvalues);
+    dcAvg.SetTime(x1, true);
+    dcAvgStd.SetTime(x1, true);
+end
+
+
+
+