Skip to content

Commit

Permalink
new file
Browse files Browse the repository at this point in the history
  • Loading branch information
@annfal
annfal committed Dec 21, 2023
1 parent 592b270 commit 2dfa698
Showing 1 changed file with 273 additions and 0 deletions.
273 changes: 273 additions & 0 deletions FuncRegistry/UserFunctions/hmrR_TaskNorm.m
View file
Original file line number Diff line number Diff line change
@@ -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




0 comments on commit 2dfa698

Please sign in to comment.