pl66tn

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

 INPUT:  x=time series (must be column array)
         dt=sample interval time [hrs] (Default dt=1)
         T=filter half-amp period [hrs] (Default T=33)

 OUTPUT: xf=filtered series

0001 function xf=pl66tn(x,dt,T);
0002 % PL66TN: pl66t for variable dt and T
0003 % xf=PL66TN(x,dt,T) computes low-passed series xf from x
0004 % using pl66 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 (must be column array)
0009 %         dt=sample interval time [hrs] (Default dt=1)
0010 %         T=filter half-amp period [hrs] (Default T=33)
0011 %
0012 % OUTPUT: xf=filtered series
0013 
0014 % NOTE: both pl64 and pl66 have the same 33 hr filter
0015 % half-amplitude period. pl66 includes additional filter weights
0016 % upto and including the fourth zero crossing at 2*T hrs.
0017 
0018 % The PL64 filter is described on p. 21, Rosenfeld (1983), WHOI
0019 % Technical Report 85-35. Filter half amplitude period = 33 hrs.,
0020 % half power period = 38 hrs. The time series x is folded over
0021 % and cosine tapered at each end to return a filtered time series
0022 % xf of the same length. Assumes length of x greater than 132.
0023 
0024 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0025 % 10/30/00
0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0027 
0028 % default to pl64
0029 if (nargin==1); dt=1; T=33; end
0030 
0031 cutoff=T/dt;
0032 fq=1./cutoff;
0033 nw=2*T./dt;
0034 nw=round(nw);
0035 disp(['number of weights = ',int2str(nw)])
0036 nw2=2.*nw;
0037 
0038 [npts,ncol]=size(x);
0039 if (npts<ncol);x=x';[npts,ncol]=size(x);end
0040 xf=x;
0041 
0042 % generate filter weights
0043 j=1:nw;
0044 t=pi.*j;
0045 den=fq.*fq.*t.^3;
0046 wts=(2.*sin(2.*fq.*t)-sin(fq.*t)-sin(3.*fq.*t))./den;
0047 % make symmetric filter weights
0048 wts=[wts(nw:-1:1),2.*fq,wts];
0049 wts=wts./sum(wts);% normalize to exactly one
0050 % plot(wts);grid;
0051 % title(['pl64t filter weights for dt = ',num2str(dt),' and T = ',num2str(T)])
0052 % xlabel(['number of weights = ',int2str(nw)]);pause;
0053 
0054 % fold tapered time series on each end
0055 cs=cos(t'./nw2);
0056 jm=[nw:-1:1];
0057 
0058 for ic=1:ncol
0059 % ['column #',num2str(ic)]
0060  jgd=find(~isnan(x(:,ic)));
0061  npts=length(jgd);
0062  if (npts>nw2)
0063 %detrend time series, then add trend back after filtering
0064   xdt=detrend(x(jgd,ic));
0065   trnd=x(jgd,ic)-xdt;
0066   y=[cs(jm).*xdt(jm);xdt;cs(j).*xdt(npts+1-j)];
0067 % filter
0068   yf=filter(wts,1.0,y);
0069 % strip off extra points
0070   xf(jgd,ic)=yf(nw2+1:npts+nw2);
0071 % add back trend
0072   xf(jgd,ic)=xf(jgd,ic)+trnd;
0073  else
0074  'warning time series is too short'
0075  end
0076 end