% ex33a
% load Excel sheet as table, rename to "T"

% filter out all freq above 1hz
lpFilt = designfilt('lowpassfir','PassbandFrequency',0.002, ...
         'StopbandFrequency',0.004,'PassbandRipple',0.5, ...
         'StopbandAttenuation',65,'DesignMethod','kaiserwin');
% 1 hz = 0.002; 2= 0.004; 4=0.008; 6=0.012; 10=0.02

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

% this loop takes around 5+ mins to run; but since the table is in RAM, the
% graphs are then fast to plot and explore
for i=2:11
T.(i)=filtfilt(lpFilt,T.(i));
end

figure
plot(T.Time,T.s01,'DisplayName','s01 rST36');
hold on;
axis( [5 300 4.25 4.4])
yyaxis right
plot(T.Time,T.s02,'DisplayName','s02 rST36 +6mm');
axis( [5 300 4.77 4.92])
yyaxis left
xlabel('Time (seconds)')
ylabel('Impedance at 40KHz (K\Omega)')

legend('show', 'Location', 'northwest')
grid
title('Right ST-36 and +6mm')

figure
plot(T.Time,T.s08,'DisplayName','s08 rGB34');
hold on
axis( [5 300 2.87 3.05])
yyaxis right
plot(T.Time,T.s07,'DisplayName','s07 rGB34 +6mm');
axis( [5 300 3.11 3.29])

yyaxis left
xlabel('Time (seconds)')
ylabel('Impedance at 40KHz (K\Omega)')
legend('show', 'Location', 'southwest')
grid
title('Right GB-34 and +6mm')

% 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.s08(indexOfInterest)) % plot on new axes
axis([121 128 2.9095 2.935 ])
%axis tight
yyaxis right
plot(T.Time(indexOfInterest),T.s07(indexOfInterest))
axis([121 128 3.235 3.261 ])
title('Close-up of feature at 124 seconds')

% plot cross section of GB34
% import ex33aAPcoords, remane as "A"
figure
x=[A.x(1),A.x2(1)];
y=[A.y(1),A.y2(1)];
% 71 175 26 ; 0.773 0.684 0.102
plot(x,y,'-','Color',[0.773 0.684 0.102]);
hold on
x=[A.x(2),A.x2(2)];
y=[A.y(2),A.y2(2)];
plot(x,y,'-','Color',[0.773 0.684 0.102]);
x=[A.x(3),A.x2(3)];
y=[A.y(3),A.y2(3)];
plot(x,y,'-','Color',[0.277 0.473 0.098]);
x=[A.x(4),A.x2(4)];
y=[A.y(4),A.y2(4)];
plot(x,y,'-','Color',[0.277 0.473 0.098]);
axis ([-6.5 6.5 2.88  3.29])
xlabel('Distance from acupoint (mm)')
ylabel('Impedance at 40KHz (K\Omega)')
title ('Cross section of GB34')
grid


figure
plot(T.Time,T.s08,'DisplayName','s08 rGB34');
hold on
yyaxis right
plot(T.Time,T.s09,'DisplayName','s09 lGB34');
yyaxis left
xlabel('Time (seconds)')
ylabel('Impedance at 40KHz (K\Omega)')
legend('show', 'Location', 'southwest')
grid
title('Right and left GB-34')


figure
plot(T.Time,T.s09,'DisplayName','s09 lGB34');
hold on
yyaxis right
plot(T.Time,T.s10,'DisplayName','s10 lGB34 +6mm');
yyaxis left
xlabel('Time (seconds)')
ylabel('Impedance at 40KHz (K\Omega)')
legend('show', 'Location', 'southwest')
grid
title('Left GB-34 and +6mm')

%inset1
axes('position',[.65 .175 .25 .25])
box on % put box around new pair of axes
indexOfInterest = (T.Time > 158) & (T.Time < 187); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s09(indexOfInterest)) % plot on new axes
%axis([157 187 5.833 5.893])
%axis tight
yyaxis right
plot(T.Time(indexOfInterest),T.s10(indexOfInterest))
%axis([157 187 6.067 6.097 ])
title('Close-up of feature at 165 seconds')

%inset2
axes('position',[.65 .175 .25 .25])
box on % put box around new pair of axes
indexOfInterest = (T.Time > 196) & (T.Time < 209); % range of t near perturbation
plot(T.Time(indexOfInterest),T.s09(indexOfInterest)) % plot on new axes
%axis([157 187 5.833 5.893])
%axis tight
yyaxis right
plot(T.Time(indexOfInterest),T.s10(indexOfInterest))
%axis([157 187 6.067 6.097 ])
title('Close-up of feature at 200 seconds')


figure
title('Subject B, right and left ST36')
% plot(T.Time,T.s03,'DisplayName','s03 lST36');
% %yyaxis right
% plot(T.Time,T.s04,'DisplayName','s04 lST36 +6mm');
% %yyaxis left
% 
% figure
% plot(T.Time,T.s05,'DisplayName','s05 rNAP');
% hold on
% %yyaxis right
% plot(T.Time,T.s06,'DisplayName','s06 RNAP +6mm');
% 
% plot(T.Time,T.s09,'DisplayName','s09 lGB34');
% hold on
% %yyaxis right
% plot(T.Time,T.s10,'DisplayName','s10 lGB34 +6mm');
% 
% 
