% 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;