[6563] in s-news-athena
median filtering of time series data
daemon@ATHENA.MIT.EDU (Steve Wofsy)
Fri Mar 3 22:19:28 1995
Date: Fri, 3 Mar 95 22:04:18 -0500
From: scw@ariel.harvard.edu (Steve Wofsy)
To: s-news@utstat.toronto.edu
Median filtering of time-series data.
We often desire to apply a median filter to a long time series
of data, presented in order. This can be achieved using apply()
or tapply(), but the procedure is extremely slow.
The attached function does the job efficiently, using a Fortran
function from Numerical Recipes (Press et al.) wrapped with a
suitable Splus function.
Feedback appreciated.
scw
#Splus function starts here -------------------------------------
dyn.load2("medfil.o")
medfil<-function(co2x,utx,tint){
if(tint<=0)stop("tint > 0 required")
#function to median filter time series co2x, data at times utx,
#into bins tint long (units of utx)
tint<-as.double(tint)
if(!is.double(co2x))co2x<-as.double(co2x)
if(!is.double(utx))utx<-as.double(utx)
maxints<-as.integer(diff(range(utx,na.rm=T))/tint +1)
x<-as.double(rep(0,maxints))
tmean<-x
nout<-maxints
nx<-maxints
nco2x<-length(co2x)
if(nco2x!=length(utx))stop("lengths of data and time series must match")
vv<-.Fortran("medfil",co2x,utx,x,tmean,nout,tint,nco2x,nx)
nout<-vv[[5]]
list(y=vv[[3]][1:nout],t=vv[[4]][1:nout])
}
#------------------------------------------------------------------
#usage: result<-medfil(YY,TT,tint) ; returns a list.
#to create medfil.o: f77 -c -G0 medfil.f
#Splus function ends here -------------------------------------
C **** medfil.f **** fortran starts here
subroutine medfil(co2x,utx,x,tmean,nout,tint,nco2x,nx)
c median filter in interval bins size tint
c (or smaller in case there are fewer data).
c co2x, data to be filtered; utx, time line for co2x
c x=filtered results, tmean=times for filtered results
c nco2x=length of co2x, nx=length of x ; nout=number filtered data pts
real*8 co2x(nco2x),utx(nco2x)
real*8 x(nx),tmean(nx),median,tint,time
integer nx,nco2x,interval,ndrop,nstart
c number of complete calibrations
i=1
j=0
j0=1
time=utx(1)
do 1001 k=1,nco2x
c--------------
25 if(utx(k).le.time+tint.and.k.lt.nco2x)then
j=j+1
c this one in the interval
else
c index of first value to be median/d
nmedian=j
x(i)=median(nmedian,co2x(j0))
c median fcn must not re-order
tmean(i)=median(nmedian,utx(j0))
c print99,j0,j,i,co2x(j0),utx(j0),x(i),tmean(i)
99 format(3i5,8f8.2)
time=utx(k-1)
j0=k
j=1
i=i+1
endif
c--------------
1001 continue
nout=i-1
return
end
function median(n,art)
c median of n points; declare median==real.
code from p227 of Press et al.,, 1st ed.
real*8 arr(1000),art(1000),median,amedian,a
integer n
do 1999 k=1,n
if(n.gt.1000)then
print*,'Error: # pts in median (',n,') exceeds 1000'
stop
endif
arr(k)=art(k)
1999 CONTINUE
do 12 j=2,n
a=arr(j)
do 11 i=j-1,1,-1
if(arr(i).le.a) goto 10
arr(i+1)=arr(i)
11 continue
i=0
10 arr(i+1)=a
12 continue
n2=n/2
if (n2.eq.float(n)/2) amedian = (arr(n2)+arr(n2+1))/2.
if (n2.ne.float(n)/2) amedian = arr(n2+1)
median =amedian
2 return
end
C **** medfil.f **** fortran ends here
