function [bdata,nbad,nfixed,delmean,delvar,stdr]=...
    deglitch(bdata,nsd,ndt,iplot);
% Identify and remove suspicious points from ADV data
% [bdata,nbad,nfixed,delmean,delvar,stdr]=deglitch(bdata,nsd,ndt,iplot);
% 
% Flag and deglitch ADV data.
% Try using:
% nsd = 2.8; % num std. deviations that defines outliers
% ndt = 6; % cutoff frequency
% min_corr = 80;
% min_amp  = 140;
% Last revised 5/12/03 on desktop by Chris Sherwood, USGS

% number of points
n = bdata.sampleperburst;
u = bdata.Vx{:};
v = bdata.Vy{:};
w = bdata.Vz{:};
dt = 1 ./ (bdata.samprate);

% bad test - sticks out from smoothed time series
% low-pass filter (ends not treated well)
uf=lpfilt(u,dt,1/(ndt*dt));
vf=lpfilt(v,dt,1/(ndt*dt));
wf=lpfilt(w,dt,1/(ndt*dt));
X = (1:n)';

% other criteria
min_amp = 120;
min_corr = 65;

a = bdata.amp{1};
c = bdata.corr{1};

if(iplot),
  figure(1)
  clf
  vBins = [-100:1:100];
  [Nu, vBins]= hist(u(:),vBins);
  [Nv, vBins]= hist(v(:),vBins);
  [Nw, vBins]= hist(w(:),vBins);
  umean = mean(u);
  ustd  = std(u);
  vmean = mean(v);
  vstd  = std(v);
  wmean = mean(w);
  wstd  = std(w);
  
  subplot(331)
  bar(vBins,100*Nu./n);
  axis([-100.5 100.5 0 20])
  set(gca,'XTick',[-100:25:100])
  ylabel('Percent')
  hold on;
  plot([umean,umean],[0 20],'-r')
  plot([umean+nsd*ustd,umean+nsd*ustd],[0 20],'--r')
  plot([umean-nsd*ustd,umean-nsd*ustd],[0 20],'--r')
  subplot(334)
  bar(vBins,100*Nv./n);
  axis([-100.5 100.5 0 20])
  set(gca,'XTick',[-100:25:100])
  hold on
  plot([vmean,vmean],[0 20],'-r')
  plot([vmean+nsd*vstd,vmean+nsd*vstd],[0 20],'--r')
  plot([vmean-nsd*vstd,vmean-nsd*vstd],[0 20],'--r')
  ylabel('Percent')

  subplot(337)
  bar(vBins,100*Nw./n);
  axis([-50.5 50.5 0 20])
  set(gca,'XTick',[-50:10:50])
  hold on
  plot([wmean,wmean],[0 20],'-r')
  plot([wmean+nsd*wstd,wmean+nsd*wstd],[0 20],'--r')
  plot([wmean-nsd*wstd,wmean-nsd*wstd],[0 20],'--r')
  ylabel('Percent')
  xlabel('Velocity')
  
  % Correlation histograms
  cBins = [0:1:100];
  worst_corr = min(min(c))
  [Nc1, cBins]= hist(c(:,1),cBins);
  [Nc2, cBins]= hist(c(:,2),cBins);
  [Nc3, cBins]= hist(c(:,3),cBins);
  
  subplot(332)
  bar(cBins,100*Nc1./n);
  axis([.5 100.5 0 20])
  set(gca,'XTick',[0:20:100])
  ylabel('Percent')
  hold on
  plot([min_corr,min_corr],[0 20],'-r')
  subplot(335)
  bar(cBins,100*Nc2./n);
  axis([.5 100.5 0 20])
  set(gca,'XTick',[0:20:100])
  ylabel('Percent')
  hold on
  plot([min_corr,min_corr],[0 20],'-r')
  subplot(338)
  bar(cBins,100*Nc3./n);
  axis([.5 100.5 0 20])
  set(gca,'XTick',[0:20:100])
  ylabel('Percent')
  hold on
  plot([min_corr,min_corr],[0 20],'-r')
  xlabel('Correlation')
  
  % Amplitude histograms
  aBins = [0:5:250];
  worst_amp = min(min(a))
  [Na1, aBins]= hist(a(:,1),aBins);
  [Na2, aBins]= hist(a(:,2),aBins);
  [Na3, aBins]= hist(a(:,3),aBins);
  subplot(333)
  bar(aBins,100*Na1./n);
  axis([.5 250.5 0 20])
  set(gca,'XTick',[0:50:250])
  hold on
  plot([min_amp,min_amp],[0 20],'-r')
  subplot(336)
  bar(aBins,100*Na2./n);
  axis([.5 250.5 0 20])
  set(gca,'XTick',[0:50:250])
  hold on
  plot([min_amp,min_amp],[0 20],'-r')
  subplot(339)
  bar(aBins,100*Na3./n);
  axis([.5 250.5 0 20])
  set(gca,'XTick',[0:50:250])
  hold on
  plot([min_amp,min_amp],[0 20],'-r')
  xlabel('Amplitude')
  drawnow;
end
  
nfixed = 0;
ur = u;
vr = v;
wr = w;
warning off MATLAB:polyfit:RepeatedPointsOrRescale;
uresid = u-uf;
vresid = v-vf;
wresid = w-wf;
stdur=std(uresid);
stdvr=std(vresid);
stdwr=std(wresid);
stdr = [stdur,stdvr,stdwr];
fprintf(1,'std devs of residuals %g %g %g\n',stdr);
% Could use other flags to find replacement points,
% but I think the low-pass filter works pretty well
bad = ( (abs(uresid) > (nsd*stdur)) | ( abs(vresid) > (nsd*stdvr) ));
bad_corr = ( (c(:,1)< min_corr)|(c(:,2)<min_corr)|(c(:,3)<min_corr) );
bad_amp  = ( (a(:,1)< min_amp )|(a(:,2)<min_amp )|(a(:,3)<min_amp ) );

rep = find(bad);
nrep = find(~bad);
if(any(rep)),
  for kk=1:length(rep),
    if(rep(kk)>(7) & rep(kk)<n-7), % skip endpoints
      % surrounding values
      lopts =  ((rep(kk)-7):(rep(kk)-1));
      hipts =  ((rep(kk)+1):(rep(kk)+7));
      blopts = bad(lopts,1);
      bhipts = bad(hipts,1);
      lopts = lopts(find(blopts~=1));
      hipts = hipts(find(bhipts~=1));
      if( (sum(blopts)<6) & (sum(bhipts)<6) ),
	pts = [lopts hipts];
	% at least 2 adjacent points on each side are ok...
	% replace with quadratic fit
	%[pts' X(pts',1) u(pts',1) v(pts',1) w(pts',1)];
	P = polyfit(X(pts,1),u(pts,1),2);
	ur(rep(kk),1)= polyval(P,X(rep(kk),1));
	P = polyfit(X(pts,1),v(pts,1),2);
	vr(rep(kk),1)=polyval(P,X(rep(kk),1));
	P = polyfit(X(pts,1),w(pts,1),2);
	wr(rep(kk),1)=polyval(P,X(rep(kk),1));
	nfixed=nfixed+1;
      else
	fprintf(1,'Warning: bad pts surrounting kk=%d\n',rep(kk))
      end
    end
  end
end
nbad = sum(bad);
fprintf(1,'Number of points flagged: %d\n',nbad );
fprintf(1,'Number. of endpoints fixed: %d\n',nfixed );
fprintf(1,'Statistics: before (after replacement) [after removal].\n');
m = [ mean([ u, ur]), mean( ur(nrep) ) ];
fprintf(1,'ur     Mean: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   m, 100*[m(2)/m(1), m(3)/m(1)] )
vv = [ std([ u, ur]), std( ur(nrep) ) ].^2;
fprintf(1,'ur Variance: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   vv, 100*[vv(2)/vv(1), vv(3)/vv(1)] );
mx = [max([ u, ur ]), max(u(nrep))];
mn = [min([ u, ur ]), min(u(nrep))];
fprintf(1,'ur      Max: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   mx, 100*[mx(2)/mx(1), mx(3)/mx(1)] );
fprintf(1,'ur      Min: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   mn, 100*[mn(2)/mn(1), mn(3)/mn(1)] );
fprintf('\n')
delmean = zeros(3,1);
delvar  = zeros(3,1);
delmean(1)=m(3)/m(1);
delvar(1) = vv(3)/vv(1);

m = [ mean([ v, vr]), mean( v(nrep) ) ];
fprintf(1,'vr     Mean: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   m, 100*[m(2)/m(1), m(3)/m(1)] )
vv = [ std([ v, vr]), std( v(nrep) ) ].^2;
fprintf(1,'vr Variance: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   vv, 100*[vv(2)/vv(1), vv(3)/vv(1)] );
mx = [max([ v, vr ]), max(v(nrep))];
mn = [min([ v, vr ]), min(v(nrep))];
fprintf(1,'vr      Max: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   mx, 100*[mx(2)/mx(1), mx(3)/mx(1)] );
fprintf(1,'vr      Min: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   mn, 100*[mn(2)/mn(1), mn(3)/mn(1)] );
fprintf('\n')
delmean(2)=m(3)/m(1);
delvar(2) = vv(3)/vv(1);

m = [ mean([ w, wr]), mean( wr(nrep) ) ];
fprintf(1,'wr     Mean: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   m, 100*[m(2)/m(1), m(3)/m(1)] )
vv = [ std([ w, wr]), std( w(nrep) ) ].^2;
fprintf(1,'wr Variance: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   vv, 100*[vv(2)/vv(1), vv(3)/vv(1)] );
mx = [max([ w, wr ]), max(w(nrep))];
mn = [min([ w, wr ]), min(w(nrep))];
fprintf(1,'wr      Max: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   mx, 100*[mx(2)/mx(1), mx(3)/mx(1)] );
fprintf(1,'wr      Min: %9f (%9f) [%9f]; Pct. orig: (%f) [%f]\n',...
   mn, 100*[mn(2)/mn(1), mn(3)/mn(1)] );
fprintf('\n')
delmean(3)=m(3)/m(1);
delvar(3) = vv(3)/vv(1);

if(iplot),
% Plots u,v, smoothed u,v, and flagged points
dt = 0.5;
bt = (0:dt:n*dt-dt)';
nd(:,1) = (u-mean(u))./std(u);
nd(:,2) = (v-mean(v))./std(v);
nd(:,3) = (w-mean(w))./std(w);
ndr(:,1) = (ur-mean(u))./std(u);
ndr(:,2) = (vr-mean(v))./std(v);
ndr(:,3) = (wr-mean(w))./std(w);
ndf(:,1) = lpfilt(nd(:,1),dt, 1/(nsd*dt));
ndf(:,2) = lpfilt(nd(:,2),dt ,1/(nsd*dt));
ndf(:,3) = lpfilt(nd(:,3),dt ,1/(nsd*dt));
sym = ['.g';'.m';'oc';'+g';'xr'];
nsda = 4.5;
npts=200;
for ich=1:n/npts;
  ie = npts*ich; is = ie-npts+1;

  figure(2); clf
  subplot(211)
  plot(bt(is:ie),u(is:ie,1),'-k')
  hold on
  plot(bt(is:ie),v(is:ie,1),'-b')
  plot(bt(is:ie),w(is:ie,1),'-c')
  
  subplot(212)
  plot(bt(is:ie),nd(is:ie,1),'.k')
  axis([bt(is) bt(ie) -nsda nsda])
  xlabel('Time (s)')
  ylabel('Std. Devs. about Detrended ur, vr')
  hold on
  plot(bt(is:ie),nd(is:ie,2),'.b')
  plot(bt(is:ie),ndr(is:ie,1),'-k')
  plot(bt(is:ie),ndr(is:ie,2),'-b')
  plot(bt(is:ie),ndf(is:ie,1),'--k')
  plot(bt(is:ie),ndf(is:ie,2),'--b')
  blist = find( bad(:,1) );
  if(any(blist)),
    plot(bt(blist),nd(blist,1),'xr');
    plot(bt(blist),nd(blist,2),'xr');
  end
  drawnow
  shg
  pause
end
end

bdata.Vx{:}=ur;
bdata.Vy{:}=vr;
bdata.Vz{:}=wr;