%
% function to clean bad pts and fill gaps in data
% uses a combination of despiking and spline fit routines
%
%
% INPUT:
%    where yd = yearday
%           x = raw data
%        hrs1  = despike window width desired
%        hrs2  = spline fit interval desired (or time vector)
%        nstd = number of standard deviations for acceptable
%               data (default=2)
%        miss  = missing value
%
% OUTPUT:
%         newt = new equally spaced yearday vector
%           y = cleaned & equally spaced data
%        npts = number of points in each despike window
% 
function [newt,y,npts]=clean(yd,x,hrs1,hrs2,nstd,miss);
days=hrs1/24;
%
if nargin< 5; nstd=2; end;
if nargin< 6; miss=999; end;
%
% First, replace bad pts > nstd*std of a "hrs1" window with mean of that
% window
%
for i=1:length(yd);
   id=find(yd>yd(i)-days/2&yd<yd(i)+days/2);
   id2=find(x(id)~=miss);
   npts(i)=length(id(id2));
   m=mean(x(id(id2)));
   s=std(x(id(id2)));
   if abs(x(i)-m)>nstd*s;
     y1(i)=m;
   elseif x(i)==miss;
     y1(i)=NaN;
   else;
     y1(i)=x(i);
   end;
end;
%
%
% now, we want to fill the gaps
% try spline fit to regular spaced pts
%
if length(hrs2)>1;
 newt=hrs2;
else;
 newt=min(yd):hrs2/24:max(yd);
end;
%
y2=spline(yd,y1,newt);
%
% plot the clean data
%
subplot(311);
plot(yd,x,'r',yd,y1,'g');
title('raw (red) & despiked (green)');
subplot(312);
plot(yd,x,'r');
hold on;
plot(newt,y2,'g');
title('spline fit');
hold off;
%
% Now, redo the despiking to rid of artificial spline peaks
%
count=0;
%
for i=1:length(y2);
 id2=find(yd>newt(i)-days/2&yd<newt(i)+days/2);
 id3=find(y1(id2)~=NaN);
%
% Keep track of missing values count
%
   if id2==[]|id3==[];
     count=count+1;
   elseif length(id2(id3))==0;
     count=count+1;
   else
     m=mean(x(id2(id3)));
     s=std(x(id2(id3)));
     if count>0;
        if length(hrs2)>1;
          misstim=mean(diff(hrs2))*24*count+(days/2)*24.0;
        else;  
          misstim=hrs2*count+(days/2)*24.0;
        end;
       disp(['Note: ',num2str(misstim/24),' day hole in the time series']);
%       k=menu(['Warning: Missing ',num2str(misstim),' hours prior to yearday ',num2str(newt(i))] ...
%              ,'Fill with overall mean' ...
%              ,'Fill with zeroes' ...
%              ,'Fill with "NaN" (not a number)' ...
%              ,'trust the despiked spline fit');
       k=3;
       if k==1;
          for jj=i-1:-1:i-count;
          y(jj)=mean(x);
          end;
       elseif k==2;
          for jj=i-1:-1:i-count;
          y(jj)=0;
          end;
       elseif k==3;
          for jj=i-1:-1:i-count;
          y(jj)=NaN;
          end;
        end; 
        count=0;
     end;
    end;
%
   if abs(y2(i)-m)>nstd*s;
     y(i)=m;
   else
     y(i)=y2(i);
   end;
end;
%
% Now, plot the despiked, gap-filled, spline fitted data
%
subplot(313);
plot(yd,x,'r');
title('despiked spline fit');
hold on;
plot(newt,y,'g');
%disp('press return to continue');
pause(2);
hold off;

