function [f,ibig] = medfilt(v, fl, tol);
% [f,ibig] = medfilt(v, fl, tol);
% f is filtered (replace outliers with median)
% ibig is indices of those replaced
% 
% v is array of values to filter
% fl is filter length   (no half-width)
% tol is tolerance (defines outliers)


if rem(fl,2) ~= 1
  error('Filter length must be odd-valued')
end

[mrow ncol] = size(v);
if mrow ~= 1 & ncol ~= 1
  error('Input must be a vector')
end

if mrow > 1     % want to work with row vector
  f = v';
  ncol = mrow;
else
  f = v;
end

half_fl = fix(fl/2);

% create array of staggered vectors:
w = NaN * ones(fl,ncol+fl-1);
for ii = 1:fl
  w(ii,ii:ii+ncol-1) = f;
end

% use median of middle section:
m = median_nan(w(:,half_fl+1:half_fl+ncol));

% check against tolerance; replace with median if needed
ibig = find(abs(f-m) > tol);
if ~isempty(ibig)
  f(ibig) = m(ibig);
end

% retranspose if needed
if mrow > 1 
  f = f';
end

