function compute_neff,x,nlag=nlag,tfilter=tfilter,tscale=tscale ; ; ***Although there are no programming errors, as far as I know, the ; ***method would seem to be in error, since neff(raw) is always greater ; ***than neff(hi) plus neff(low) - which shouldn't be true, otherwise ; ***some information has somehow been lost. For now, therefore, run ; ***compute_neff for unfiltered series, then for low-pass series, and ; ***subtract the results to obtain the neff for the high-pass series! ; ; The aim of this function is to estimate the effective independent sample ; size of a filtered time series. The result of this function is [neff,rfac], ; where neff=effective independent sample size, which equals the sample size ; divided by the reduction factor, rfac. The original timeseries is x(t). ; The filter to be applied is given as tfilter=[t1,t2], where t1=0 implies ; a high pass filter with cutoff t2, t2=999 implies a low pass filter with ; cutoff t1, and t1<>0 and t2<>999 implies a band pass filter from t1 to t2. ; If tfilter is not specified or if tfilter=[0,999] then no filtering is done. ; If tscale=[s1,s2] is specified, then a fraction (s1) of the high frequency ; variability is added back on to the series, and a fraction (s2) of the low ; frequency variability is added back on to the series. ; ; The sample size itself is computed by computing the lag-1 autocorrelation ; of x(t), then generating a synthetic 10000-term AR1 series, filtering this ; as described above, and then computing the lag-1 to lag-nlag (default 30) ; autocorrelations and evaluating 1+sum(abs(rk)). ; ; Check inputs and set defaults ; if n_elements(nlag) eq 0 then nlag=30 if (nlag lt 1) or (nlag gt 200) then message,'Wrong number of lags' ; if n_elements(tfilter) eq 0 then tfilter=[0,999] if n_elements(tfilter) ne 2 then message,'Specify two filter cutoffs' if tfilter(0) ge tfilter(1) then message,'Filter cutoffs in wrong order' ; if (tfilter(0) eq 0) and (tfilter(1) ne 999) then $ print,'Not sure that it is working quite right for high-pass filters!' ; if n_elements(tscale) eq 0 then tscale=[0,0] if n_elements(tscale) ne 2 then message,'Specify two filter scaling factors' if (tscale(0) lt 0) or (tscale(0) gt 1) then message,'Scaling factor wrong' if (tscale(1) lt 0) or (tscale(1) gt 1) then message,'Scaling factor wrong' ; nx=n_elements(x) seed=98541L nlen=20000L ; ; Compute lag1r of the input series ; lagr=a_correlate(x,[1]) ar1=lagr(0) ; ; Generate the AR1 series ; ; Need a slightly longer series to avoid end-effects with the filtering if tfilter(1) eq 999 then nplus=tfilter(0) else nplus=tfilter(1) if (nplus mod 2) eq 1 then nplus=nplus+1 nextra=nlen+nplus wnoi=randomn(seed,nextra) y=fltarr(nextra) y(0)=wnoi(0) for i = 1L , nextra-1 do begin y(i)=ar1*y(i-1)+wnoi(i) endfor xval=findgen(nlag) ; ; Now filter it appropriately ; t1=tfilter(0) if t1 gt 0 then begin filter_cru,t1,tsin=y,tshigh=yh1,tslow=yl1 endif else begin yh1=y*0. yl1=y endelse ; t2=tfilter(1) if t2 ne 999 then begin filter_cru,t2,tsin=y,tshigh=yh2,tslow=yl2 endif else begin yh2=y yl2=y*0. endelse ; yfilt=yl1-yl2+tscale(0)*yh1+tscale(1)*yl2 ; ; Remove end segments to avoid end-effects ; nph=nplus/2 y=yfilt(nph:nph+nlen-1) ; ; Now evaluate the sample size according to the autocorrelation function of ; the filtered series. Because the autocorrelation doesn't converge to zero ; at increasing lags (due to limited sample length and resultant sampling ; noise), we need to remove the effect of any small but non-zero ; autocorrelation coefficient - e.g., anything less than 0.02 in magnitude. ; lagnr=a_correlate(y,indgen(nlag)+1) zlist=where(abs(lagnr) lt 0.1,nzero) if nzero gt 0 then lagnr(zlist)=0. ;print,lagnr ;print,total(lagnr),total(lagnr > 0),total(abs(lagnr)) rrr=abs(lagnr) rfac=1.+2.*total(rrr) return,[nx/rfac,rfac] ; end