Home > bobstuff > pl64tap.m

pl64tap

PURPOSE ^

PL64TAP low-pass filters the hourly time series.

SYNOPSIS ^

function xf=pl64tap(x);

DESCRIPTION ^

 PL64TAP low-pass filters the hourly time series.

 xf=PL64TAP(x) low-pass filters the hourly time series x using 
 the PL64 filter described on p. 21, Rosenfeld (1983), WHOI 
 Technical Report 85-35. Filter half amplitude period = 33 hrs.,
 half power period = 38 hrs. The time series x is folded over 
 and cosine tapered at each end to return a filtered time series 
 xf of the same length.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ver. 1: 3/2/95 (SL)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001    function xf=pl64tap(x);
0002 
0003 % PL64TAP low-pass filters the hourly time series.
0004 %
0005 % xf=PL64TAP(x) low-pass filters the hourly time series x using
0006 % the PL64 filter described on p. 21, Rosenfeld (1983), WHOI
0007 % Technical Report 85-35. Filter half amplitude period = 33 hrs.,
0008 % half power period = 38 hrs. The time series x is folded over
0009 % and cosine tapered at each end to return a filtered time series
0010 % xf of the same length.
0011 %
0012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0013 % ver. 1: 3/2/95 (SL)
0014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0015 
0016 [npts,ncol]=size(x);
0017 
0018 if (npts>126) 
0019 
0020 %generate filter weights
0021 j=1:63;
0022 t=pi.*j;
0023 den=0.0009.*t.^3;
0024 wts=(2.*sin(0.06.*t)-sin(0.03.*t)-sin(0.09.*t))./den;
0025 % make symmetric filter weights
0026 % coefficient is to match fortran version of outputs
0027 wts=0.999531323.*[wts(63:-1:1),0.05997,wts];
0028 
0029 %fold tapered time series on each end
0030 cs=cos(t'./126);
0031 jm=[63:-1:1];
0032 
0033 % demean x
0034 xmean=x;
0035 for nc=1:ncol;
0036   x(:,nc)=x(:,nc)-xmean(nc);
0037 end
0038 % detrend
0039 xin=x;
0040 xt=zeros(npts,ncol);
0041 for nc=1:ncol;
0042   x(:,nc)=detrend(xin(:,nc));
0043   xt(:,nc)=xin(:,nc)-x(:,nc);
0044 end
0045 
0046 for ic=1:ncol
0047 ['column #',num2str(ic)]
0048 
0049  y=[cs(jm).*x(jm,ic);x(:,ic);cs(j).*x(npts+1-j,ic)];
0050 
0051 % filter
0052  yf=filter(wts,1.0,y);
0053 
0054 % strip off extra points
0055  xf(:,ic)=yf(127:npts+126);
0056  end
0057 
0058 % add mean and trend
0059 for nc=1:ncol;
0060 xf(:,nc)=xf(:,nc)+xt(:,nc)+xmean(nc);
0061 end
0062 
0063 else
0064 'warning less than 127 points'
0065 
0066 end
0067 end
0068 
0069

Generated on Thu 01-Dec-2005 07:46:03 by m2html © 2003