From 2dfa69802696f261787b0f03f66fccca0cefce15 Mon Sep 17 00:00:00 2001 From: annfal Date: Thu, 21 Dec 2023 14:20:48 +0100 Subject: [PATCH] new file --- FuncRegistry/UserFunctions/hmrR_TaskNorm.m | 273 +++++++++++++++++++++ 1 file changed, 273 insertions(+) create mode 100644 FuncRegistry/UserFunctions/hmrR_TaskNorm.m diff --git a/FuncRegistry/UserFunctions/hmrR_TaskNorm.m b/FuncRegistry/UserFunctions/hmrR_TaskNorm.m new file mode 100644 index 00000000..81fb8b2b --- /dev/null +++ b/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 + + + +