% ex40
% import temp readings as table, rename to TT0
% import imp readings as table, rename to T0


% The thermistor's thermal time constant is given as 5 s, measured in air
% But when comparing the readings to the markers, and also the artefact on
% left Lu1, there apears to be no significant delay in the temperature
% response, and certainly nothing = 5 s, or even near to 1 s. I
% corresponded with the manufacturer, and they said the time constant
% depends on the application, and a time constant of around 100 ms, would
% be normal in some circumstances. Therefore, no allowance has been made in
% the readings for any delay in the temperature readings.


% filter out any equipment artifact at about 0.7 to1.2 Hz,
lpFilt = designfilt('lowpassfir','PassbandFrequency',0.001, ...
         'StopbandFrequency',0.0015,'PassbandRipple',0.5, ...
         'StopbandAttenuation',65,'DesignMethod','kaiserwin');

     % 0.4 hz= 0.0008; 0.5hz =0.001; 1 hz = 0.002; 2= 0.004; 3=0.006; 4=0.008; 6=0.012; 10=0.02
TT=TT0;
for i=2:3
TT.(i)=filtfilt(lpFilt,TT.(i));
end


% low pass filter from 0.5 Hz, to filter out pulse artefact, and all above
lpFilt = designfilt('lowpassfir','PassbandFrequency',0.0008, ...
         'StopbandFrequency',0.001,'PassbandRipple',0.5, ...
         'StopbandAttenuation',65,'DesignMethod','kaiserwin');
% 0.4 hz= 0.0008; 0.5hz =0.001; 1 hz = 0.002; 2= 0.004; 3=0.006; 4=0.008; 6=0.012; 10=0.02
     
T1=T0;
for i=2:13
T1.(i)=filtfilt(lpFilt,T1.(i));
end


% low pass filter from 2 Hz, to show any pulse artefact
lpFilt = designfilt('lowpassfir','PassbandFrequency',0.004, ...
         'StopbandFrequency',0.006,'PassbandRipple',0.5, ...
         'StopbandAttenuation',65,'DesignMethod','kaiserwin');
% 0.4 hz= 0.0008; 0.5hz =0.001; 1 hz = 0.002; 2= 0.004; 3=0.006; 4=0.008; 6=0.012; 10=0.02

TP=T0;
for i=12:13
TP.(i)=filtfilt(lpFilt,TP.(i));
end


% ==== right lu9 and pulse artefact
T=TP;
figure
plot(T.Time,T.s12,'DisplayName','r LU9, lowpass filter 2 Hz');
hold on;
yyaxis right;
T=T1;
plot(T.Time,T.s12,'DisplayName','r LU9, lowpass filter 0.5 Hz');

yyaxis left;
xlabel('Time (seconds)') 
ylabel('Impedance at 40 kHz (K\Omega)')
legend('show', 'Location', 'northwest')
title('Right Lu-1, showing pulse artifact')
grid
xlim([75 105])


axes('position',[0.13 0.11 0.775 0.815]);
box on % put box around new pair of axes
plot (TT.Time,TT.s02,'Color',[0.4660    0.6740    0.1880],'DisplayName','breath at nostril');
hold on
axis([75 105 30 40 ])
%axis tight
legend('show', 'Location', 'north');
grid

%====== peaks
[pks,locs] = findpeaks(TT.s02);
locs=locs/1000;


% ============== LU9 ===========
T=T1;
figure
set(groot,'defaultLineLineWidth',1)
set(gca,'FontSize',14)
hold all;

plot(T.Time,T.s12,'DisplayName','r LU9');
hold on;
yyaxis right;
plot(T.Time,T.s11,'DisplayName','r LU9 +6mm');

yyaxis left
xlabel('Time (seconds)') 
ylabel('Impedance at 40 kHz (K\Omega)')
legend('show', 'Location', 'northwest')
title('Impedance at right LU-9 and at +6 mm, and breathing pattern; lowpass filter 0.5 Hz')
grid

yMax=45;
yMin=0;
yyaxis right
marker=168.107;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','hold1');

marker=180.165;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','end1');

marker=225.156;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.91 0.41 0.17],'DisplayName','hold2');

marker=250.662;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.91 0.41 0.17],'DisplayName','end2');

% closeup box
axes('position',[.65 .175 .25 .25])
box on % put box around new pair of axes
indexOfInterest = (T.Time > 166) & (T.Time < 182); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s12(indexOfInterest)) % plot on new axes
hold on
%axis([120 128 1.62 1.635 ])
%axis tight
yyaxis right
plot(T.Time(indexOfInterest),T.s11(indexOfInterest))
hold on
%axis([120 128 2.875 2.89 ])
title('Close-up between first marks')
grid

yMax=3;
yMin=2;
yyaxis right
marker=168.107;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','hold1');

marker=180.165;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','end1');

% breath at nostril
axes('position',[0.13 0.11 0.775 0.815]);
box on % put box around new pair of axes
plot (TT.Time,TT.s02,'Color',[0.4660    0.6740    0.1880],'DisplayName','breath at nostril');
hold on
axis([0 300 17.9 44.9 ])
%axis tight
legend('show', 'Location', 'north');


% LU1, ========================

T=T1;
figure
set(groot,'defaultLineLineWidth',1)
set(gca,'FontSize',14)
hold all;

plot(T.Time,T.s01,'DisplayName','r LU1');
hold on;
plot(T.Time,T.s02,'DisplayName','r LU1 +6mm');

plot(T.Time,T.s03,'DisplayName','l LU1');
plot(T.Time,T.s04,'DisplayName','l LU1 +6mm');

% plot(T.Time,T.s08,'DisplayName','r LU6');
% hold on;
% yyaxis right;
% plot(T.Time,T.s07,'DisplayName','r LU6 +6mm');
% 
% plot(T.Time,T.s10,'DisplayName','l LU6');
% hold on;
% yyaxis right;
% plot(T.Time,T.s09,'DisplayName','l LU6 +6mm');
% 
% plot(T.Time,T.s12,'DisplayName','r LU9');
% hold on;
% yyaxis right;
% plot(T.Time,T.s11,'DisplayName','r LU9 +6mm');


yyaxis right;
plot(TT.Time,TT.s02,'DisplayName','Breath at nostril');
hold on;

ylabel('Breath at nostril (Celsius)')

yyaxis left
xlabel('Time (seconds)') 
ylabel('Impedance at 40 kHz (K\Omega)')
legend('show', 'Location', 'northwest')
title('Impedance at left and right Lu-1, compared to breathing pattern; lowpass filter 0.5 Hz')
grid


yMax=42;
yMin=20;
yyaxis right
marker=168.107;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','hold1');

marker=180.165;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','end1');

marker=225.156;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.91 0.41 0.17],'DisplayName','hold2');

marker=250.662;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.91 0.41 0.17],'DisplayName','end2');


% breath at nostril
axes('position',[0.13 0.11 0.775 0.815]);
box on % put box around new pair of axes
plot (TT.Time,TT.s02,'Color',[0.4660    0.6740    0.1880],'DisplayName','breath at nostril');
hold on
axis([60 90 30 40 ])
%axis tight
legend('show', 'Location', 'north');
grid

% closeup box
axes('position',[.65 .175 .25 .25])
box on % put box around new pair of axes
indexOfInterest = (T.Time > 60) & (T.Time < 90); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s01(indexOfInterest)) % plot on new axes
hold on
%axis([120 128 1.62 1.635 ])
%axis tight
yyaxis right
plot(T.Time(indexOfInterest),T.s02(indexOfInterest))
hold on
%axis([120 128 2.875 2.89 ])
title('Close-up of right LU-1')
grid

% ==== LU6,left ====
T=T1;
figure
set(groot,'defaultLineLineWidth',1)
set(gca,'FontSize',14)
hold all;

plot(T.Time,T.s10,'DisplayName','l LU6');
hold on;
yyaxis right;
plot(T.Time,T.s09,'DisplayName','l LU6 +6mm');

yyaxis left
xlabel('Time (seconds)') 
ylabel('Impedance at 40 kHz (K\Omega)')
legend('show', 'Location', 'northwest')
title('Impedance at left LU-6 and at 6 mm radius, compared to breathing pattern')
grid

yMax=5;
yMin=4;
yyaxis right
marker=168.107;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','hold1');

marker=180.165;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','end1');

marker=225.156;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.91 0.41 0.17],'DisplayName','hold2');

marker=250.662;
x=[marker,marker];
y=[yMin,yMax];
plot(x,y,'--','Color',[0.91 0.41 0.17],'DisplayName','end2');

% breath at nostril
axes('position',[0.13 0.11 0.775 0.815]);
box on % put box around new pair of axes
plot (TT.Time,TT.s02,'Color',[0.5 0.5 0.5],'DisplayName','breath at nostril');
hold on
axis([0 300 17.9 44.9 ])
%axis tight
legend('show', 'Location', 'north');


%====== peaks
[pks,locs] = findpeaks(T1.s02); %
locs=locs/1000;

[L6pks,L6locs] = findpeaks(T.s07);
L6locs=L6locs/1000;


%============= r LU1 =6mm invertd
rLU1invert= -T.s02;

[pks,locs] = findpeaks(rLU1invert); %
locs=locs/1000;





