% ex34 impedance
% import excel sheet as table; rename to "T0"
% import ex34 Ren12ache; rename to "R12"

% the following filtering may take some time, but since the results are
% then held in memory, the plots can be produce and changed quickly

% filter out the equipment artifact at about 0.7 to1.2 Hz,
lpFilt = designfilt('lowpassfir','PassbandFrequency',0.001, ...
         'StopbandFrequency',0.0014,'PassbandRipple',0.5, ...
         'StopbandAttenuation',65,'DesignMethod','kaiserwin');
   
T1=T0;
for i=2:13
T1.(i)=filtfilt(lpFilt,T1.(i));
end  

%LineWidth
set(groot,'defaultLineLineWidth',1)

% ===== left ST36 and +6mm ============
T=T1;
figure
plot(T.Time,T.s04,'DisplayName','lST36');
hold on;
yyaxis right
plot(T.Time,T.s03,'DisplayName','lST36 +6mm');

yyaxis left
y=1.58;
y1=y+.25;
axis( [110 150 y y1])

% xlabel('Time (seconds)')
ylabel('Impedance at40kHz (K\Omega)')
legend('show', 'Location', 'northwest')
grid
title('Left ST-36 and +6mm; lowpass filter at 0.7Hz')
% marker

yyaxis right
marker=121.821;
x=[marker,marker];
y=[1,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[1,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');

y=2.83;
y1=y+.25;
axis( [110 150 y y1])



% ======== next plot; closeup and inset ============

T=T1;
% left ST36 
figure
plot(T.Time,T.s04,'DisplayName','lST36');
hold on;
yyaxis right
plot(T.Time,T.s03,'DisplayName','lST36 +6Med');
yyaxis left

y=1.58;
y1=y+0.06;
axis( [110 480 y y1])

% xlabel('Time (seconds)')
ylabel('Impedance @40KHz (KOhms)')
xlabel('Time (seconds)')
legend('show', 'Location', 'northwest')
grid
title('Left ST-36 and +6mm medial; lowpass filter 0.7Hz; signals superimposed')
% marker

yyaxis right
marker=121.821;
x=[marker,marker];
y=[1,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[1,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');

y=2.85;
y1=y+0.06;
axis( [110 150 y y1])
% left ST36 and Ren12 ache

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

marker=121.821;
x=[marker,marker];
y=[1,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[1,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');


% ===ST19 left=================
% left ST19 = s011 12; right st19 = s09 10

T=T1;
figure
plot(T.Time,T.s11,'DisplayName','lST19');
hold on;
yyaxis right
plot(T.Time,T.s12,'DisplayName','lST19 +6med');
hold on;

yyaxis left
ylabel('Impedance @40KHz (KOhms)')
xlabel('Time (seconds)')
legend('show', 'Location', 'southwest')
grid
title('Left ST-19 and +6mm medial; lowpass filter 0.7Hz')

y1=0.931;
y2=y1+0.35;
axis ([110 480 y1 y2])

yyaxis right
y1=5.973;
y2=y1+0.35;
axis ([110 480 y1 y2])

% marker
marker=121.821;
x=[marker,marker];
y=[5,6.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[5,6.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');


% closeup box
axes('position',[.65 .175 .25 .25])
box on % put box around new pair of axes
indexOfInterest = (T.Time > 121) & (T.Time < 128); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s11(indexOfInterest)) % plot on new axes
hold on
y1=1.11;
y2=y1+0.04;
axis([121 128 y1 y2])
yyaxis right
plot(T.Time(indexOfInterest),T.s12(indexOfInterest))
hold on
y1=6.283;
y2=y1+0.04;
axis([121 128 y1 y2])
title('Close-up of features between marks')
grid

marker=121.821;
x=[marker,marker];
y=[6,7];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[6,7];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');


%====== right st19
T=T1;
figure
plot(T.Time,T.s10,'DisplayName','rST19');
hold on;
yyaxis right
plot(T.Time,T.s09,'DisplayName','rST19 +6med');
hold on;

yyaxis left
ylabel('Impedance @40KHz (KOhms)')
xlabel('Time (seconds)')
legend('show', 'Location', 'northwest')
grid
title('Right ST-19 and +6mm medial; lowpass filter 0.7Hz')

y1=4.95;
y2=y1+0.35;
axis ([110 480 y1 y2])

yyaxis right
y1=5.85;
y2=y1+0.35;
axis ([110 480 y1 y2])

% marker
marker=121.821;
x=[marker,marker];
y=[1,7];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[1,7];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');


% closeup box
T=T1;
axes('position',[.65 .175 .25 .25])
box on % put box around new pair of axes
indexOfInterest = (T.Time > 121) & (T.Time < 139); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s10(indexOfInterest)) % plot on new axes
hold on
y1=5.03;
y2=y1+0.04;
axis([121 139 y1 y2])
yyaxis right
plot(T.Time(indexOfInterest),T.s09(indexOfInterest))
hold on
y1=6.1;
y2=y1+0.04;
axis([121 139 y1 y2])
title('Close-up of features between marks')
grid

marker=121.821;
x=[marker,marker];
y=[6,7];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[6,7];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');


%======= L st21
T=T1;

figure
plot(T.Time,T.s07,'DisplayName','lST21');
hold on;
yyaxis right
plot(T.Time,T.s08,'DisplayName','lST21 +6med');

yyaxis left
xlabel('Time (seconds)')
ylabel('Impedance @40KHz (KOhms)')
legend('show', 'Location', 'northwest')
grid
title('Left ST-21 and +6mm medial; lowpass filter at 4Hz')

y=2.634;
y1=y+0.3;
axis( [110 480 y y1])

yyaxis right
y=3.68;
y1=y+0.3;
axis( [110 480 y y1])

%marker
marker=121.821;
x=[marker,marker];
y=[3.5,4.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[3.5,4.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');


% closeup box
T=T1;
axes('position',[.65 .175 .25 .25])
box on % put box around new pair of axes
indexOfInterest = (T.Time > 121) & (T.Time < 130); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s07(indexOfInterest)) % plot on new axes
hold on

y1=2.75;
y2=y1+0.03;
axis([121 130 y1 y2])
yyaxis right
plot(T.Time(indexOfInterest),T.s08(indexOfInterest))
hold on

y1=3.86;
y2=y1+0.03;
axis([121 130 y1 y2])
title('Close-up of features between marks')
grid

marker=121.821;
x=[marker,marker];
y=[3,4];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[3,4];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');



%====== rST21,
T=T1;
figure
plot(T.Time,T.s05,'DisplayName','rST21');
hold on;
yyaxis right
plot(T.Time,T.s06,'DisplayName','rST21 +6mm');

yyaxis left
xlabel('Time (seconds)')
ylabel('Impedance at 40kHz (K\Omega)')
legend('show', 'Location', 'northwest')
grid
title('Right ST-21 and +6mm; lowpass filter at 0.7Hz')

yyaxis right
% marker
marker=121.821;
x=[marker,marker];
y=[5.8,6.2];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark 1');

marker=127.148;
x=[marker,marker];
y=[5.8,6.2];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark 2');

y=3.74;
y1=y+0.18;
axis( [110 480 y y1])

yyaxis right
y=5.88;
y1=y+0.18;
axis( [110 480 y y1])

%====== rst21 close up; Fig04
T=T1;
figure
plot(T.Time,T.s05,'DisplayName','rST21');
hold on;
yyaxis right
plot(T.Time,T.s06,'DisplayName','rST21 +6mm');

yyaxis left
xlabel('Time (seconds)')
ylabel('Impedance at 40KHz (k\Omega)')
legend('show', 'Location', 'west')
grid
title('Right ST-21 and +6mm; lowpass filter at 0.7Hz; signals superimposed')

y1=3.749;
y2=y1+0.148;
axis([98 215 y1 y2])

% marker
yyaxis right
y1=5.929;
y2=y1+0.148;
axis([98 215 y1 y2])

marker=121.821;
x=[marker,marker];
y=[5.5,6.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[5.5,6.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');

% closeup box
axes('position',[.65 .175 .25 .25])
box on % put box around new pair of axes
indexOfInterest = (T.Time > 117) & (T.Time < 131); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s05(indexOfInterest)) % plot on new axes
hold on
y1=3.84;
y2=y1+0.05;
axis([117 131 y1 y2])

yyaxis right
plot(T.Time(indexOfInterest),T.s06(indexOfInterest))
hold on
y1=6.01;
y2=y1+0.05;
axis([117 131 y1 y2])
title('Signals superimposed at marks')

marker=121.821;
x=[marker,marker];
y=[6,6.1];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[6,6.1];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');


%======= rst36 filt 1hz =======================

T=T1;

figure
plot(T.Time,T.s01,'DisplayName','rST36');
hold on;
yyaxis right
plot(T.Time,T.s02,'DisplayName','rST36 +6mm');

yyaxis left
xlabel('Time (seconds)')
ylabel('Impedance at 40KHz (k\Omega)')
legend('show', 'Location', 'northwest')
grid
title('Right ST-36 and +6mm; lowpass filter at 0.7Hz')

% marker
yyaxis right
marker=121.821;
x=[marker,marker];
y=[4,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[4,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');
axis([110 480 4.747 5.247])
yyaxis left
axis([110 480 1.388 1.478])

% ren12
subplot(2,1,2)
stairs(R12.Time,R12.Ren12, 'Linewidth',2,'DisplayName', 'Ren12 ache')
ylabel('Intensity')
xlabel('Time (seconds)')
legend('show')
grid
axis([110 480 1 10])

% closeup box
axes('position',[.65 .175 .25 .25])
box on % put box around new pair of axes
indexOfInterest = (T.Time > 110) & (T.Time < 150); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s01(indexOfInterest)) % plot on new axes
hold on
y1=1.396;
y2=y1+0.01;
axis([110 150 y1 y2])

yyaxis right
plot(T.Time(indexOfInterest),T.s02(indexOfInterest))
hold on
y1=5.19;
y2=y1+0.01;
axis([110 150 y1 y2])
title('Signals superimposed at marks')

marker=121.821;
x=[marker,marker];
y=[5,5.3];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[5,5.3];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');



% close up
% ALTERNATIVE FILTER, to show fine detail when zoomed in     
% filter out all freq above 6-10hz
lpFilt = designfilt('lowpassfir','PassbandFrequency',0.012, ...
         'StopbandFrequency',0.02,'PassbandRipple',0.5, ...
         'StopbandAttenuation',65,'DesignMethod','kaiserwin');
% 1 hz = 0.002; 2= 0.004; 4=0.008; 6=0.012; 10=0.02; 12=0.024;
% 16=0.032; 28=0.056; 34=0.068; 46=0.092 
T=T0;
for i=2:3
T.(i)=filtfilt(lpFilt,T.(i));
end

figure
indexOfInterest = (T.Time > 121) & (T.Time < 130); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s01(indexOfInterest),'DisplayName','rST36') % plot on new axes
hold on
y1=1.398;
y2=y1+0.012;
axis([120 128 y1 y2])

yyaxis right
plot(T.Time(indexOfInterest),T.s02(indexOfInterest),'DisplayName','rST36 +6mm')
hold on
y1=5.19;
y2=y1+0.012;
axis([120 128 y1 y2])
yyaxis left
xlabel('Time (seconds)')
ylabel('Impedance @40KHz (KOhms)')
legend('show', 'Location', 'northwest')
title('Right ST-36; features between marks; signals superimposed')
grid

yyaxis right
marker=121.821;
x=[marker,marker];
y=[5,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark1');

marker=127.148;
x=[marker,marker];
y=[5,5.5];
plot(x,y,'--','Color',[0.277 0.473 0.098],'DisplayName','mark2');

% ============= peaks ================
%
% plot(T.Time,T.s01,'DisplayName','rST36');
% plot(T.Time,T.s02,'DisplayName','rST36 +6mm');
% plot(T.Time,T.s04,'DisplayName','lST36');
% plot(T.Time,T.s03,'DisplayName','lST36 +6Med');
% 
% plot(T.Time,T.s05,'DisplayName','rST21');
% plot(T.Time,T.s06,'DisplayName','rST21 +6mm med');
% plot(T.Time,T.s07,'DisplayName','lST21');
% plot(T.Time,T.s08,'DisplayName','lST21 +6med');
% 
% plot(T.Time,T.s11,'DisplayName','lST19');
% plot(T.Time,T.s12,'DisplayName','lST19 +6med');
% plot(T.Time,T.s10,'DisplayName','rST19');
% plot(T.Time,T.s09,'DisplayName','rST19 +6med');

%========= peaks
T=T1;

[R21pks,R21locs] = findpeaks(T.s05);
R21locs=R21locs/1000;

[L21pks,L21locs] = findpeaks(T.s07);
L21locs=L21locs/1000;

[R19pks,R19locs] = findpeaks(T.s10);
R19locs=R19locs/1000;

[L19pks,L19locs] = findpeaks(T.s11);
L19locs=L19locs/1000;

[R36pks,R36locs] = findpeaks(T.s01);
R36locs=R36locs/1000;

[L36pks,L36locs] = findpeaks(T.s04);
L36locs=L36locs/1000;
