pl33tn

 PL33TN: pl33t with variable dt and T
 xf=PL33TN(x,dt,T) computes low-passed series xf from x
 using pl33 filter, with optional sample interval dt (hrs)
 and filter half-amplitude period T (hrs) as input for
 non-hourly series.

 INPUT:  x=time series
         dt=sample interval [hrs] (Default: dt=1)
         T=filter half-amplitude period [hrs] (Default: T=33)

 OUTPUT: xf=filtered series

0001 function xf=pl33tn(x,dt,T);
0002 % PL33TN: pl33t with variable dt and T
0003 % xf=PL33TN(x,dt,T) computes low-passed series xf from x
0004 % using pl33 filter, with optional sample interval dt (hrs)
0005 % and filter half-amplitude period T (hrs) as input for
0006 % non-hourly series.
0007 %
0008 % INPUT:  x=time series
0009 %         dt=sample interval [hrs] (Default: dt=1)
0010 %         T=filter half-amplitude period [hrs] (Default: T=33)
0011 %
0012 % OUTPUT: xf=filtered series
0013 
0014 % The PL33 filter is described on p. 21, Rosenfeld (1983), WHOI
0015 % Technical Report 85-35. Filter half amplitude period = 33 hrs.,
0016 % half power period = 38 hrs. The time series x is folded over
0017 % and cosine tapered at each end to return a filtered time series
0018 % xf of the same length. Assumes length of x greater than 67.
0019 
0020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 % 10/30/00
0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0023 
0024 % default to pl33
0025 if (nargin==1); dt=1;T=33; end
0026 
0027 cutoff=T/dt;
0028 fq=1./cutoff;
0029 nw=T./dt;
0030 nw=floor(cutoff);
0031 disp(['number of weights = ',int2str(nw)])
0032 nw2=2.*nw;
0033 
0034 [npts,ncol]=size(x);
0035 if (npts<ncol);x=x';[npts,ncol]=size(x);end
0036 xf=x;
0037 
0038 % generate filter weights
0039 j=1:nw;
0040 t=pi.*j;
0041 den=fq.*fq.*t.^3; % assumes f1=f2/2
0042 wts=(2.*sin(2.*fq.*t)-sin(fq.*t)-sin(3.*fq.*t))./den;
0043 % make symmetric filter weights
0044 wts=[wts(nw:-1:1),2.*fq,wts];
0045 wts=wts./sum(wts);% normalize to exactly one
0046 % plot(wts)
0047 
0048 % fold tapered time series on each end
0049 cs=cos(t'./nw2);
0050 jm=[nw:-1:1];
0051 
0052 for ic=1:ncol
0053 % ['column #',num2str(ic)]
0054  jgd=find(~isnan(x(:,ic)));
0055  npts=length(jgd);
0056 
0057  if (npts>nw2)
0058 % detrend time series, then add trend back after filtering
0059   xdt=detrend(x(jgd,ic));
0060   trnd=x(jgd,ic)-xdt;
0061   y=[cs(jm).*xdt(jm);xdt;cs(j).*xdt(npts+1-j)];
0062 
0063 % filter
0064   yf=filter(wts,1.0,y);
0065 
0066 % strip off extra points
0067   xf(jgd,ic)=yf(nw2+1:npts+nw2);
0068 
0069 % add back trend
0070   xf(jgd,ic)=xf(jgd,ic)+trnd;
0071 
0072  else
0073  'warning time series is too short'
0074  end
0075 
0076 end