function [x_hub, varargout] = huber(x)
% weighted arithmetic mean with Huber weights
% mean = huber(x); 
% [mean, sigma] = huber(x); 
% [mean, sigma, sum_w] = huber(x);
% [mean, sigma, sum_w, w] = huber(x); 
% [mean, sigma, sum_w, w, rms] = huber(x); 
%
% July 2002, Nils Olsen, DSRI

x = x(~isnan(x)); % remove NaN
epsilon = 0.001;  % interrupt iteration if relative change < epsilon
Breakpoint = 1.5;
% Breakpoint = 9e9;

x_old = mean(x);
c = Breakpoint * std(x);
w = min(c./abs(x - x_old), 1);
sum_w = sum(w);
x_hub = sum(w.*x)/sum_w;     
sigma = sqrt(sum(w.*(x-x_hub).^2)/sum_w);
rms = sqrt(sum(w.*x.^2)/sum_w);
c = Breakpoint * sigma;
while abs((x_old - x_hub)/x_hub) > epsilon
   x_old = x_hub;
   w = min(c./abs(x - x_old), 1);
   sum_w = sum(w);
   x_hub = sum(w.*x)/sum_w;
   sigma = sqrt(sum(w.*(x-x_hub).^2)/sum_w);
   rms = sqrt(sum(w.*x.^2)/sum_w);
   c = Breakpoint * sigma;
end

if nargout > 1; varargout{1} = sigma; end
if nargout > 2; varargout{2} = sum_w; end
if nargout > 3; varargout{3} = w; end
if nargout > 4; varargout{4} = rms; end
