function [ no_spike y_smooth ] = nospike( x, y, n )

%function [ no_spike y_smooth ] = nospike( x, y, n )
%
%   This function finds spikes in y (>5 sigma excursions) and then smooths 
%   the remaining data across n points.  It returns no_spike, which is the 
%   index to the y elements that are not spikes.  y_smooth is the smoothed
%   version of y, incorporating NaNs for data points that are spikes.
%   
%   J. G. Bellingham, November 25, 2013

%%
% Do median filter then figure out distribution of differences from the
% median value

y_med = medfilt1(y, 5);
y_dif = (y - y_med);

% Now construct distribution of deviations from median
comparison = std(y_dif(find(~isnan(y_dif))));

i_spike = find(abs(y_dif) > 5*comparison);
no_spike = find(abs(y_dif) <= 5*comparison);

%%
% Do moving line fit...this seems to be a very interesting way to smooth
% data.

if n > 1
    [y_slope, y_intercept, y_fit_error] = slope(x, y, n);
    yy = x .* y_slope + y_intercept;
    
    % version wo spikes...
    
    [y_slope, y_intercept, y_fit_error] = slope(x(no_spike), y(no_spike), 30);
    y_smooth = NaN * yy;
    y_smooth(no_spike) = x(no_spike) .* y_slope + y_intercept;
else
    y_smooth = y(no_spike);
end;

