statok2.m 1.08 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
function [xm,xs,xok]=statok2(x,varargin)

dimx=size(x);
if nargin < 3, dim=min(find(dimx ~= 1)); else dim=varargin{2}; end
if nargin < 2, weight=1; else weight=varargin{1}; end
if isempty(dim), dim = 1; end

if dim == 1
    if nargin < 4, thresh=NaN; else thresh=varargin{3}; end
    nc=prod(dimx(2:end));
    xp=reshape(x,dimx(1),nc);
    xokp=isfinite(xp);
    [xpm,xps]=deal(NaN*ones(1,nc));
    for j=1:nc
    	xq=xp(:,j);
	ok=xokp(:,j);
	nok=sum(ok);
	go=nok>0;
	if nok
	while go
	    xqm=sum(xq(ok),1)/nok;
	    xqs=sqrt(var(xq(ok),weight,1));
	    xq0=xq;
	    xq0(ok)=xqm;
	    out=abs(xq-xq0)>(thresh*xqs);
	    ok=ok & (~out);
	    nok=sum(ok);
	    go=(sum(out)>0) && (nok>0);
	end
	xpm(j)=xqm;
	xps(j)=xqs;
	xokp(:,j)=ok;
	else
	[xpm(j),xps(j)]=deal(NaN);
	xokp(:,j)=false;
	end
    end
    xm=reshape(xpm,[1 dimx(2:end)]);
    xs=reshape(xps,[1 dimx(2:end)]);
    xok=reshape(xokp,dimx);
else
   perm=1:ndims(x); perm(dim)=[]; perm=[dim perm];
   [xm1,xs1,xok1]= statok2(permute(x,perm),weight,1,varargin{3:end});
   xm=ipermute(xm1,perm);
   xs=ipermute(xs1,perm);
   xok=ipermute(xok1,perm);
end