function [y,ok] = std2(varargin)
%STD2 Standard deviation, ignoring NaN values.
% For vectors, Y = STD2(X) returns the standard deviation. For matrices,
% Y is a row vector containing the standard deviation of each column. For
% N-D arrays, STD operates along the first non-singleton dimension of X.
%
% STD2 normalizes Y by N, where N is the sample size. It produces the square
% root of the second moment of the sample about its mean.
% STD2(X,1) is the same as STD(X).
%
% Y = STD2(X,0) normalizes by N-1. This is the sqrt of an unbiased estimator
% of the variance of the population from which X is drawn, as long as X
% consists of independent, identically distributed samples.
%
% Y = STD2(X,FLAG,DIM) takes the standard deviation along the dimension
% DIM of X. Pass in FLAG==1 to use the default normalization by N, or
% 0 to use N-1.
%
% Example: If X = [4 -2 1
% 9 5 7]
% then std2(X,0,1) is [3.5355 4.9497 4.2426] and std2(X,0,2) is [3.0
% 2.0]
% See also MEAN2.
if nargin > 3
x=varargin{1};
w=varargin{2};
dim=varargin{3};
ok=isfinite(x);
tile = ones(1,max(ndims(x),dim)); tile(dim) = size(x,dim);
threshold=varargin{4};
while true
y=sqrt(var2(x,w,dim));
overs=abs(x-repmat(mean2(x,dim),tile)) > repmat(threshold*y,tile);
if ~any(overs(:)), break; end
ok(overs)=false;
x(overs)=NaN;
end
else
y=sqrt(var2(varargin{:}));
end