function [ufilt,jdf]=plfilt(u,jd,nsamp)
% PLFILT  low-pass filters hourly data using the pl33 filter
%         and cuts off 33 hours at beginning and end
%
%  USAGE: [ulp]=plfilt(u)
%     or
%         [ulp,jdlp]=plfilt(u,jd)
%     or
%         [ulp,jdlp]=plfilt(u,jd,nsub)
%
%    where u is a matrix with time series in each column and
%    jd is the optional julian day vector, and nsub subsamples
%    the output each nsub points.

load pl33
[m,n]=size(u);
if(min(m,n)==1),
  m=length(u);
  n=1;
  u=u(:);
end
if (m<length(pl33)),
  disp('data too short for this filter!!');
  return
end
ufilt=zeros(m-66,n);
for i=1:n,
  uf=conv(u(:,i),pl33);
  ufilt(:,i)=uf(67:length(uf)-66);
end
[m,n]=size(ufilt);
if (nargin>1),
  jdf=jd(34:34+m-1);
end
if (nargin==3)
 ufilt=ufilt([1:nsamp:m],:);
 jdf=jdf([1:nsamp:m]);
end
