% Glasgow Caledonian University: 2Nd November 2015

% Glasgow Caledonian University: 2nd November 2015

% Laura McKernan Ward

% Call in text file data

clc, clear all %cleans workspace to start afresh

FileName = uigetfile('*.txt','Select the M-File'); %opens up dialog box to select text file

FileID=fopen (FileName);%opens selected file - the data files MUST be in the same pathname

C = textscan(FileID,'%s %s %s %s %s %s %s %s %s %s %s'); %reads data from open text file into cell array C

fclose(FileID);

% Extracting correct data

D= [C{2} C{3} C{4} C{5} C{6} C{7} C{8} C{9} C{10} C{11}]; %creates a data file from text file with chromophores, time + marker data

D(1:4,:)=[]; %deleting emptpy first 5 rows so left with raw data

data=cellfun(@str2num,D); %data to be manipulated in matrix format

clear C D FileID ans

%% Arranging B0 data

ZD = data(1:70,10) < 1 & data(1:70,10);%raw values of baseline

B=[data(ZD,:)];B0=flipud(B);

B0=B0(1:20,:);B0=flipud(B0); % selecting last 20s of baseline stimulation

rawBmean=mean(B0); % mean of B0

normalised_baseline = B0 - rawBmean(ones(size(B0,1),1),:);

windowSize = 8;

nb = filter(ones(1,windowSize)/windowSize,1,normalised_baseline); % baseline post filtering of moving average filter

nb=[(1:20)',nb(:,2:10)];

clear off B PathName FileID rbm ZD sd

%% Deleting B0 from data and Normalise Data around B0

x=data(:,1)<70 & data(:,10)<1;

d = data(~x,:); %d=data(x,P) gives you only the baseline data

time =[1:length(d)]';

data=[time, d(:,2:10)];

nd = data(:,2:9) - repmat(rawBmean(:,2:9),size(data(:,1))); % normalised data around baseline values

data=[data(:,1),filter(ones(1,windowSize)/windowSize,1,nd),data(:,10)]; % filtering using moving average

x=detrend(data(:,2:9));data=[data(:,1),x,data(:,10)];

clear x time b Time B0_prep B rbm windowSize

%% Section up data into experimental epochs

x= data(1:35,:);y= data(20:65,:);

On1 = data(1:35,1)<35 & data(1:35,10)>1; %creating logical vector to identify the ON period

Off1= data(20:65,1)<65 & data(20:65,10)<1; %creating logical vector to identify the OFF period

HD1 = x(On1,:); ZD1 =y(Off1,:);

x = data(55:95,:); y = data(85:125,:);On2 = data(55:95,1)<100 & data(55:95,10)>1;

Off2= data(85:125,1)<130 & data(85:125,10)<1;HD2 =x(On2,:); ZD2 = y(Off2,:);

x = data(115:155,:); y=data(145:185,:);On3 = data(115:155,1)<155 & data(115:155,10)>1;

Off3= data(145:185,1)<185 & data(145:185,10)<1; HD3 = x(On3,:); ZD3 = y(Off3,:);

x= data(170:215,:); y = data(205:245,:);On4 = data(170:215,1)<215 & data(170:215,10)>1;

Off4= data(205:245,1)<245 & data(205:245,10)<1;HD4 = x(On4,:); ZD4 =y(Off4,:);

x=data(235:275,:);y=data(265:305,:);On5 = data(235:275,1)<275 & data(235:275,10)>1;

Off5= data(265:305,1)<305 & data(265:305,10)<1; HD5 = x(On5,:); ZD5 =y(Off5,:);

x=data(290:335,:);y=data(320:365,:);On6 = data(290:335,1)<335 & data(290:335,10)>1;

Off6= data(320:365,1)<365 & data(320:365,10)<1; HD6 = x(On6,:); ZD6 = y(Off6,:);

x=data(350:395,:);y=data(380:425,:);

On7 = data(350:395,1)<395 & data(350:395,10)>1;

Off7= data(380:425,1)<425 & data(380:425,10)<1;

HD7 = x(On7,:); ZD7 =y(Off7,:);

x=data(410:455,:);y=data(440:485,:);

On8= data(410:455,1)<455 & data(410:455,10)>1;

Off8= data(440:485,1)<485 & data(440:485,10)<1;

HD8 = x(On8,:); ZD8 = y(Off8,:);

x=data(470:515,:);y=data(500:545,:);

On9= data(470:515,1)<515 & data(470:515,10)>1;

Off9= data(500:545,1)<545 & data(500:545,10)<1;

HD9 = x(On9,:); ZD9 =y(Off9,:);

x=data(530:575,:);y=data(565:end,:);

On10= data(530:575,1)<575 & data(530:575,10)>1;

Off10= data(565:end,1)<610 & data(565:end,10)<1;

HD10 = x(On10,:); ZD10 = y(Off10,:);

clear d On1 On2 On3 On4 On5 On6 On7 On8 On9 On10 Off1 Off2 Off3 Off4 Off5 Off6 Off7 Off8 Off9 Off10 x y time Off11 On11

%% Organising all NORMALISED data

HD{1,1}=HD1;HD{1,2}=HD2;HD{1,3}=HD3;HD{1,4}=HD4;HD{1,5}=HD5;HD{1,6}=HD6;

ZD{1,1}=ZD1;ZD{1,2}=ZD2;ZD{1,3}=ZD3;ZD{1,4}=ZD4;ZD{1,5}=ZD5;ZD{1,6}=ZD6;

HD{1,7}=HD7;HD{1,8}=HD8;HD{1,9}=HD9;HD{1,10}=HD10;

ZD{1,7}=ZD7;ZD{1,8}=ZD8;ZD{1,9}=ZD9;ZD{1,10}=ZD10;

for i=1:length(HD) % arranging data back into on/off sections

ed{1,i}=[HD{1,i};ZD{1,i}];

end

ND=[nb;ed{1,1};ed{1,2};ed{1,3};ed{1,4};ed{1,5};ed{1,6};ed{1,7};ed{1,8};ed{1,9};ed{1,10}];

clear HD1 HD2 HD3 HD4 HD5 HD6 HD7 HD8 HD9 HD10 ZD1 ZD2 ZD3 ZD4 ZD5 ZD6 ZD7 ZD8 ZD9 ZD10

clear i off rawBmean ZD_means windowSize prep data a b c d x normalised_baseline first D

%% Calculating average HRF for stimulus ZD/HD (off/on)

D1_OSt=zeros(30,10);D2_OSt=zeros(30,10);D1_HbT=zeros(30,10);D2_HbT=zeros(30,10);

D1_HbO=zeros(30,10);D1_HbR=zeros(30,10);D2_HbO=zeros(30,10);D2_HbR=zeros(30,10);

for i=1:length(HD)

D1_OSt(:,i)= HD{1,i}(1:30,2); % oxygen saturation

D1_HbT(:,i)=HD{1,i}(1:30,3); % total haemoglobin

D1_HbO(:,i)= HD{1,i}(1:30,4); % HbO

D1_HbR(:,i)= HD{1,i}(1:30,5); % HbR

D2_OSt(:,i)=HD{1,i}(1:30,6);

D2_HbT(:,i)=HD{1,i}(1:30,7);

D2_HbO(:,i)= HD{1,i}(1:30,8);

D2_HbR(:,i)= HD{1,i}(1:30,9);

end

D1_OSt_off=zeros(29,10);D2_OSt_off =zeros(29,10);D1_HbT_off =zeros(29,10);D2_HbT_off =zeros(29,10);

D1_HbO_off =zeros(29,10);D1_HbR_off =zeros(29,10);D2_HbO_off =zeros(29,10);D2_HbR_off =zeros(29,10);

for i=1:length(ZD)

D1_OSt_off (:,i)= ZD{1,i}(1:29,2); % oxygen saturation

D1_HbT_off (:,i)=ZD{1,i}(1:29,3); % total haemoglobin

D1_HbO_off (:,i)= ZD{1,i}(1:29,4); % HbO

D1_HbR_off (:,i)= ZD{1,i}(1:29,5); % HbR

D2_OSt_off (:,i)=ZD{1,i}(1:29,6);

D2_HbT_off (:,i)=ZD{1,i}(1:29,7);

D2_HbO_off (:,i)= ZD{1,i}(1:29,8);

D2_HbR_off (:,i)= ZD{1,i}(1:29,9);

end

on_avg_hdr=[mean(D1_OSt,2),mean(D1_HbT,2),mean(D1_HbO,2),mean(D1_HbR,2),mean(D2_OSt,2),mean(D2_HbT,2),mean(D2_HbO,2),mean(D2_HbR,2)]; % OxySat/HbT/HbO/HbR for each detector 1+2

off_avg_hdr=[mean(D1_OSt_off,2),mean(D1_HbT_off,2),mean(D1_HbO_off,2),mean(D1_HbR_off,2),mean(D2_OSt_off,2),mean(D2_HbT_off,2),mean(D2_HbO_off,2),mean(D2_HbR_off,2)]; % OxySat/HbT/HbO/HbR for each detector 1+2

avg_hdr=[on_avg_hdr;off_avg_hdr];%average HDR for each detectors 1+2

V1=[(avg_hdr(:,1))+(avg_hdr(:,5)),(avg_hdr(:,2))+(avg_hdr(:,6)),(avg_hdr(:,3))+(avg_hdr(:,7)),(avg_hdr(:,4))+(avg_hdr(:,8))]/2;

OSt=[(avg_hdr(:,1)),(avg_hdr(:,5))];

HbT=[(avg_hdr(:,2)),(avg_hdr(:,6))];

HbO=[(avg_hdr(:,3)),(avg_hdr(:,7))];

HbR=[(avg_hdr(:,4)),(avg_hdr(:,8))];

s=[nanstd(OSt,0,2),nanstd(HbT,0,2),nanstd(HbO,0,2),nanstd(HbR,0,2)]; % calculate SD across rows

se=[s./sqrt(2)]./2; % standard error /2 for excel

V1=[V1,se];

clear i D1_HbO D1_HbR D2_HbR D2_HbO D1_HbO_off D1_HbR_off D2_HbR_off D2_HbO_off ans OSt HbT HbO HbR D1_OSt D1_HbT D2_OSt D2_HbT D1_HbT_off D1_OSt_off D2_HbT_off D2_OSt_off

%% Plotting the HDR

close all; figure('name',FileName);hold on

set(0,'defaultAxesFontName','Calibri');set(0,'defaultTextFontName','Calibri')

set(gcf,'Color','w','un','n','pos',[0.16,0.16,.55,.55]); %sets background as white, and positioning

set(gca,'FontSize',14,'LineWidth',2) %specifies look of graph

subplot(1,2,1); hold on;title('Average HDR - Detector A','FontSize',16,'HorizontalAlignment','center')

plot(avg_hdr(:,3),'-r','LineWidth',2.5);%plot HbO for O1

hold on; subplot(1,2,1);hold on

plot(avg_hdr(:,4),'Color',[0.2 0.4 1],'LineWidth',2.5); hold on; %plot Hb for O1

ylabel('[Chromophore] change (uM)','FontSize',12) % label for y axis

xlabel('Time (s)','FontSize',12) %label for x axis

subplot(1,2,2); hold on;title('Average HDR - Detector B','FontSize',16,'HorizontalAlignment','center')

plot(avg_hdr(:,7),'-r','LineWidth',2.5);%plot HbO for O1

hold on; subplot(1,2,2);hold on

plot(avg_hdr(:,8),'Color',[0.2 0.4 1],'LineWidth',2.5); hold on; %plot Hb for O1

ylabel('[Chromophore] change (uM)','FontSize',12) % label for y axis

xlabel('Time (s)','FontSize',12) %label for x axis

%% Taking differences of normalised data sections

% Need to take {1,1} HD - ZD

for i=1:length(HD)

a=HD{1,i}(1:29,:);

b=ZD{1,i}(1:29,:);

newd{1,i}=a-b;

end

for i=1:length(newd)

new{1,i}=newd{1,i}(:,2:9);

end

%new = differences between HD - ZD

for i=1:length(new)

new_means{1,i}=nanmean(new{1,i}(14:29,:)); % flipping data to capture last 15s

end

output=cell2mat(new_means');

B0=[mean(ND(1:20,:));((std(ND(1:20,:)))/sqrt(20))/2];

table= [nanmean(output);(nanstd(output)/sqrt(length(output)))/2;B0(:,2:9)];

clear yL nbnew d a b i normalised_baseline A B diff_data new_means newd diffdata B0 DA_HbO data off_avg_hdr on_avg_hdr ZD ZD_means HD_means HD i nd HbO HbR ed HDm ZDm rawBmean s se HDout HDout2 ZDout ZDout2 HDpeaks ZDpeaks Index_B Sorted_B prep a b c d n