function x = dSpikeEliminate( x, MaxPeakWidth, tol )
% dSpikeEliminate - eliminates spike noise in 1D arrays
%	y = dSpikeEliminate( x, MaxPeakWidth )
%	x - 1D column vector of signal data
%	MaxPeakWidth=4 (optional parameter) peaks up to this number of points wide 
%			will be found and eliminated
%	y - output vector of the same dimension as x
%	Peaks are substituted by lineary interpolated points. Algorithm eliminates peaks
%	in several runs, detecting them relatively to the standard deviation of the 
%	first derivative of the signal.
%	If more peaks are found than 10% of the length of the vector, no action is taken.

if nargin <= 1, MaxPeakWidth = 4; end;
if nargin <= 2, tol=3; end;
n = size(x,1); %length of x
flag = 1;
while (flag),
flag=0;
	dx = diff(x);
	s  = std(dx);
	for (tau=1:MaxPeakWidth),
		I = find( abs(dx) > tol*s );
		J = find( ( diff(I) == tau ));
		K = find( sign(dx(I(J)).*dx(I(J)+tau)) < 0 );
		t = I(J(K))+1;
		if (length(t) > n/20),
			flag=0; break;
		end; 
		if (~isempty(t))
			flag = 1;
			k = (x(t+tau)-x(t-1))/(tau+1);
			t2 = t*ones(1,tau) + ones(size(t,1),1)*(0:tau-1);
			x(t2) = x(t-1)*ones(1,tau) + k*(1:tau);		
		end; %if
	end; %for
end;
return;