function SmoothY=fastsmooth(Y,w,type,edge)
% fastbsmooth(Y,w,type,edge) smooths vector Y with smooth 
%  of width w. Version 2.0, April 2008.
% The argument "type" determines the smooth type:
%   If type=1, rectangular (sliding-average or boxcar) 
%   If type=2, triangular (2 passes of sliding-average)
%   If type=3, pseudo-Gaussian (3 passes of sliding-average)
% The argument "edge" controls how the "edges" of the signal 
% (the first w/2 points and the last w/2 points) are handled.
%   If edge=0, the edges are zero.  (In this mode the elapsed 
%     time is independent of the smooth width). The fastest.
%   If edge=1, the edges are smoothed with progressively 
%     smaller smooths the closer to the end. (In this mode the  
%     elapsed time increases with increasing smooth widths).
% fastsmooth(Y,w,type) smooths with edge=0.
% fastsmooth(Y,w) smooths with type=1 and edge=0.
%  T. C. O'Haver, 2008.
if nargin==2, edge=0; type=1; end
if nargin==3, edge=0; end
  switch type
    case 1
       SmoothY=sa(Y,w,edge);
    case 2   
       SmoothY=sa(sa(Y,w,edge),w,edge);
    case 3
       SmoothY=sa(sa(sa(Y,w,edge),w,edge),w,edge);
  end

function SmoothY=sa(Y,smoothwidth,edge)
w=round(smoothwidth);
SumPoints=sum(Y(1:w));
s=zeros(size(Y));
halfw=round(w/2);
for k=1:length(Y)-w,
   s(k+halfw)=SumPoints;
   SumPoints=SumPoints-Y(k);
   SumPoints=SumPoints+Y(k+w);
end
SmoothY=s./w;
% Taper the edges of the signal if edge=1.
  if edge==1,
  startpoint=(smoothwidth + 1)/2;
  L=length(Y);
  SmoothY(1)=(Y(1)+Y(2))./2;
  for k=2:startpoint,
     SmoothY(k)=mean(Y(1:(2*k-1)));
     SmoothY(L-k+1)=mean(Y(L-2*k+2:L));
  end
  SmoothY(L)=(Y(L)+Y(L-1))./2;
  end