function y = nanstd(x)
%NANSTD Standard deviation ignoring NaNs.
%   NANSTD(X) returns the same standard deviation treating NaNs
%   as missing values.
%
%   For vectors, NANSTD(X) is the standard deviation of the
%   non-NaN elements in X.  For matrices, NANSTD(X) is a row
%   vector containing the standard deviation of each column,
%   ignoring NaNs.
%
%   Y=nanstd(X,DIM) is also supported
%
%   See also NANMEAN, NANMEDIAN, NANMIN, NANMAX, NANSUM.
 
%   Copyright 1993-2002 The MathWorks, Inc.
%   $Revision: 2.10 $  $Date: 2002/01/17 21:31:14 $
%
%   Revised by Simon to include averaging over a specified dimension.
%   2004/5/19
  
if isempty(x) % Check for empty input.
    y = NaN;
    %y = []; % Simon thinks empty in -> empty out
    return
end

% put specified dimension first (SPdeS)
if exist('dim','var')
  if dim > ndims(x)
    error(['DIM exceeds number of dimensions in X. '...
	   num2str(dim) '>' num2str(ndims(x))]);
  elseif dim <= 0
    error(['DIM = ' num2str(dim) ' <= 0.']);
  elseif dim > 1
    x = shiftdim(x,dim-1); % put dim first in x
  end
end

% column-vectorize other dimensions (SPdeS)
m = size(x);
if length(m)>2
  x = reshape(x,m(1),prod(m(2:end)));
end

nans = isnan(x);
i = find(nans);
 
% Find mean
avg = nanmean(x);
 
if min(size(x))==1,
   count = length(x)-sum(nans);
   x = x - avg;
else
   count = size(x,1)-sum(nans);
   x = x - avg(ones(size(x,1),1),:);
end
 
% Replace NaNs with zeros.
x(i) = zeros(size(i));
 
% Protect against a column of all NaNs
i = find(count==0);
count(i) = ones(size(i));
y = sqrt(sum(x.*x)./max(count-1,1));
y(i) = i + NaN;

% put the dimensions of y back the way they were (SPdeS)
if length(m)>2
  y = reshape(y,[1,m(2:end)]);
end

if exist('dim','var')
  y = shiftdim(y,length(m)-(dim-1));
end