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) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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