function y = rem_outliers(x, std_lim)
%  function y = rem_outliers(x, std_lim)
%  The returned matrix Y is the same as the input matrix X but with
%  the outliers removed.  Outliers are defined as elements which are
%  more than STD_LIM standard deviations away from the mean.  Both the
%  mean and the standard deviation are calculated by ignoring any NaNs
%  in the matrix.  Note that that since the value of the standard
%  deviation is strongly influenced by outliers then this technique
%  will only be useful in those cases where there are very few
%  outliers.

%     Copyright J. V. Mansbridge, CSIRO, august 30 1993

[m, n] = size(x);
y = reshape(x, m*n, 1);
x_mean = nanmean(y);
x_std = nanstd(y);
x_norm = abs(x - x_mean)/x_std;
index = find(x_norm > std_lim);
y = x;
y(index) = nan*ones(size(index));
disp([ 'mean = ', num2str(x_mean) ' std = ', num2str(x_std) ])
disp([ num2str(length(index)) ' outliers removed'])
