function uu=despike(uz)
%disp('despik');
[np,a]=size(uz);
if a==1;
    ux=1:np;
    uz=[uz ux'];
else
    np=a;
    ux=1:a;
    uz=[uz' ux'];
end;
uu2=sortrows(uz,1);
uu=uu2(1:np,1);
mp=floor(np/2);
mu=uu(mp);
sp=floor(.84*np);
sm=floor(.16*np);
sig=(uu(sp)-uu(sm))/2;
dsig=max(4*sig,.5);
im=1;
while abs(mu-uu(im))>dsig;
   im=im+1;
end;
ip=np;
while abs(uu(ip)-mu)>dsig;
   ip=ip-1;
end;
pct=(im+np-ip)/np*100;
%disp(['total spikes= ' num2str(pct) ' %    median=  ' num2str(mu) '  sigma= ' num2str(sig)] );
uu2(1:im,1)=NaN;
uu2(ip:np,1)=NaN;
uy=sortrows(uu2,2);
%find NaN and substitute a nearest neighbor
ii=find(isnan(uy));ab=length(ii);
if ~isempty(ii)
    uy(ii(2:ab))=uy(ii(2:ab)-1);
    if ii(1)>1
        uy(ii(1))=uy(ii(1)-1);
    else
        uy(ii(1))=uy(ii(1)+1);
    end;
end;
%Do it again to eliminate double spikes
ii=find(isnan(uy));ab=length(ii);
if ~isempty(ii)
    uy(ii(2:ab))=uy(ii(2:ab)-1);
    if ii(1)>1
        uy(ii(1))=uy(ii(1)-1);
    else
        uy(ii(1))=uy(ii(1)+1);
    end;
end;

%Do it again to eliminate triple spikes
ii=find(isnan(uy));ab=length(ii);
if ~isempty(ii)
    uy(ii(2:ab))=uy(ii(2:ab)-1);
    if ii(1)>1
        uy(ii(1))=uy(ii(1)-1);
    else
        uy(ii(1))=uy(ii(1)+1);
    end;
end;

if a==1
    uu=uy(1:np,1);
else
    uu=uy(1:a,1)';
end