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