text: function [v,y,w]=nearnonz(x,d)
  %NEARNONZ replace each zero element with the nearest non-zero element [V,Y,W]=nearnonz(X,D)
  %
  %  Inputs:  x         input vector, matrix or larger array
  %           d         dimension to apply filter along [default 1st non-singleton]
  %
  % Outputs:  v         v is the same size as x but with each zero entry replaced by
  %                     the nearest non-zero value along dimension d
  %                     elements equidistant from two non-zero entries will be taken
  %                     from the higher index
  %           y         y is the same size as x and gives the index along dimension d
  %                     from which the corresponding entry in v was taken
  %                     If there are no non-zero entries, then the corresponding
  %                     elements of y will be zero.
  %           w         w is the same size as x and gives the distance (+ or -) to the
  %                     nearest non-zero entry in x
  
  %      Copyright (C) Mike Brookes 1997
  %      Version: $Id: nearnonz.m 713 2011-10-16 14:45:43Z dmb $
  %
  %   VOICEBOX is a MATLAB toolbox for speech processing.
  %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %   This program is free software; you can redistribute it and/or modify
  %   it under the terms of the GNU General Public License as published by
  %   the Free Software Foundation; either version 2 of the License, or
  %   (at your option) any later version.
  %
  %   This program is distributed in the hope that it will be useful,
  %   but WITHOUT ANY WARRANTY; without even the implied warranty of
  %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  %   GNU General Public License for more details.
  %
  %   You can obtain a copy of the GNU General Public License from
  %   http://www.gnu.org/copyleft/gpl.html or by writing to
  %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  e=size(x);
  p=prod(e);
  if nargin<2             % if no dimension given, find the first non-singleton
      d=find(e>1,1);
      if ~numel(d)
          d=1;
      end
  end
  k=e(d);                 % size of active dimension
  q=p/k;                  % size of remainder
  if d==1
      z=reshape(x,k,q);
  else
      z=shiftdim(x,d-1);
      r=size(z);
      z=reshape(z,k,q);
  end
  xx=z~=0;
  cx=cumsum(xx);
  [i,j]=find(z);
  qq=cx(xx);
  pos=full(sparse(qq,j,i,k,q)); % list the positions of non-zero elements in each column
  mp=ceil((pos(1:end-1,:)+pos(2:end,:))*0.5); % find the mid point between consecutive non-zero elements
  [i2,j2]=find(pos(2:end,:)>0);
  zz=1+cumsum(full(sparse(mp(pos(2:end,:)>0),j2,1,k,q)));
  y=pos(zz+repmat((0:q-1)*k,k,1));
  v=z(max(y,1)+repmat((0:q-1)*k,k,1));
  w=y-repmat((1:k)',1,q);
  w(y==0)=0;
  if d==1
      y=reshape(y,e);
      v=reshape(v,e);
      w=reshape(w,e);
  else
      y=shiftdim(reshape(y,r),length(e)+1-d);
      v=shiftdim(reshape(v,r),length(e)+1-d);
      w=shiftdim(reshape(w,r),length(e)+1-d);
  end
  