function s4std(channel)

if (~exist('z_data'))
    load('z_data');
end
load('s1');

fig = channel+400;
z = z_data(:,channel);

figure(fig); clf;
plot(0:.05:9999.95,z);

threshold = 4*std(z);
s4std_indices = find( abs(z) > threshold );  % finds all +ve and -ve peaks
s4std = 0.05 * s4std_indices;                % save these peak times in s4std
hold on;
plot(s4std,zeros(length(s4std),1),'k.','markersize',10);

ss1=conv(z,s1);
ss1=ss1(11:end-10);
plot(0:.05:9999.95, ss1,'r');
AX = axis;
plot( [AX(1) AX(2)], [threshold threshold], 'g' );
plot( [AX(1) AX(2)], [-threshold -threshold], 'g' );

figure(fig+1);clf;
plot(0:.05:9999.95,z);

threshold = -4*std(z);
s4std_indices = find( z < threshold );        % finds -ve peak only
s4std = 0.05 * s4std_indices - 0.05;          % save these peak times in s4std
hold on;
plot(s4std,zeros(length(s4std),1),'k.');
AX = axis;
plot( [AX(1) AX(2)], [threshold threshold], 'g' );

k=0; m=0; mi=0;
sp=[];
for i=1:length(s4std)
   if (s4std(i)-k <= 0.4)           % if we haven't jumped past a 0.4ms vicinity in this step
     if (z(s4std_indices) < m)
       m=z(s4std_indices);          % accumulates local minimum z(t)
       mi=s4std(i);                 % gets actual time of minimum
     end
   else
     sp(end+1)=mi;                  % saves the min peak time in sp
     m=z(s4std_indices);            % load in a new minimum amplitude
     mi=s4std(i);                   % and a new minimum time
   end
   k=s4std(i);                      % store the previous time in k
end
sp(end+1)=mi;
sp=sp(2:end);                       % ignore the first point

plot(sp,zeros(length(sp),1),'rx');