00001
00002 c
00003 c Copyright (c) 1986,1987,1988,1989,1990,1991,1992,1993,
00004 c by Steve McMillan, Drexel University, Philadelphia, PA.
00005 c
00006 c All rights reserved.
00007 c
00008 c Redistribution and use in source and binary forms are permitted
00009 c provided that the above copyright notice and this paragraph are
00010 c duplicated in all such forms and that any documentation,
00011 c advertising materials, and other materials related to such
00012 c distribution and use acknowledge that the software was developed
00013 c by the author named above.
00014 c
00015 c THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
00016 c IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
00017 c WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00018 c
00019
00020 subroutine smooth(t,x,xsmooth,n,dts,iopt,iw)
00021 save
00022 c
00023 c Smooth the input data, using one of a variety of methods.
00024 c
00025 c Options:
00026 c
00027 c iopt = 0: median
00028 c 1: arithmetic mean
00029 c 2: harmonic mean
00030 c 3: geometric mean
00031 c 11: clipped arithmetic mean
00032 c 12: clipped harmonic mean
00033 c 13: clipped geometric mean
00034 c
00035 c iwindow = 0: flat window
00036 c 1: triangular window
00037 c
00038 parameter (NMAX = 10000, FCLIP = 0.1)
00039 c
00040 dimension t(1),x(1),xsmooth(1),work(NMAX),tt(NMAX)
00041 c
00042 external comp
00043 c
00044 c Smooth the input data over a window of width dts.
00045 c Brute-force approach recalculates the sum at each data point!
00046 c
00047 c Window function:
00048 c
00049 window(xx) = max(0., 1. - iwindow*abs(xx/dt2))
00050 c
00051 if (dts.le.0.) return
00052 iwindow = iw
00053 if (iwindow.ne.0) iwindow = 1
00054 dt2 = .5*dts
00055 c
00056 if (iopt.gt.0.and.iopt.le.10.and.iwindow.eq.0) then
00057 c
00058 c Flat window, straight average only:
00059 c
00060 sum = 0.
00061 ns = 0
00062 j0 = 1
00063 j1 = 0
00064 c
00065 do 190 i=1,n
00066 c
00067 120 continue
00068 c
00069 c if (j1.lt.n
00070 c & .and.t(j1+1).le.t(i)+min(dt2,t(i)-t(1))) then
00071 c
00072 if (j1.lt.n.and.t(j1+1).le.t(i)+dt2) then
00073 j1 = j1 + 1
00074 sum = sum + f(x(j1),iopt)
00075 ns = ns + 1
00076 go to 120
00077 end if
00078 c
00079 130 continue
00080 c if (t(j0).lt.t(i)-min(dt2,t(n)-t(i))) then
00081 if (t(j0).lt.t(i)-dt2) then
00082 sum = sum - f(x(j0),iopt)
00083 ns = ns - 1
00084 j0 = j0 + 1
00085 go to 130
00086 end if
00087 c
00088 xsmooth(i) = finv(sum/max(ns,1),iopt)
00089 c
00090 190 continue
00091 else
00092 c
00093 c All other options:
00094 c
00095 j0 = 1
00096 j1 = 0
00097 c
00098 c Locate the range containing the desired interval.
00099 c
00100 do 290 i=1,n
00101 c
00102 220 if (j1.lt.n
00103 & .and.t(j1+1).le.t(i)+min(dt2,t(i)-t(1))) then
00104 j1 = j1 + 1
00105 go to 220
00106 end if
00107 c
00108 230 if (t(j0).lt.t(i)-min(dt2,t(n)-t(i))) then
00109 j0 = j0 + 1
00110 go to 230
00111 end if
00112 c
00113 c Range is j0 to j1.
00114 c
00115 do 240 j=j0,j1
00116 work(j-j0+1) = x(j)
00117 tt(j-j0+1) = t(j)
00118 240 continue
00119 nw = j1-j0+1
00120 c
00121 if (iopt.eq.0.or.iopt.gt.10) then
00122 c
00123 c Sort the data.
00124 c
00125 if (nw.gt.1) call sort2(nw,work,tt)
00126 c
00127 if (iopt.eq.0) then
00128 c
00129 c Interpolate to the median.
00130 c
00131 xj = .5*(nw+1)
00132 jj = xj
00133 xsmooth(i) = work(jj)
00134 & + (xj-jj)*(work(jj+1) - work(jj))
00135 else
00136 c
00137 c Clip the data.
00138 c
00139 xj1 = FCLIP*(nw+1)
00140 jj1 = max(1,nint(xj1))
00141 xj2 = (1.-FCLIP)*(nw+1)
00142 jj2 = min(nw,nint(xj2))
00143 c
00144 c New range is jj1 to jj2.
00145 c
00146 do 250 j=jj1,jj2
00147 work(j-jj1+1) = work(j)
00148 tt(j-jj1+1) = tt(j)
00149 250 continue
00150 nw = jj2-jj1+1
00151 end if
00152 end if
00153 c
00154 if (iopt.gt.0) then
00155 c
00156 c Perform averaging on the work data.
00157 c
00158 sum = 0.
00159 wsum = 0.
00160 c
00161 do 260 j=1,nw
00162 weight = window(tt(j)-t(i))
00163 wsum = wsum + weight
00164 sum = sum + weight*f(work(j),iopt)
00165 260 continue
00166 c
00167 if (wsum.gt.0.) then
00168 xsmooth(i) = finv(sum/wsum,iopt)
00169 else
00170 xsmooth(i) = x(i)
00171 end if
00172 else
00173 end if
00174 c
00175 290 continue
00176 end if
00177 c
00178 end
00179
00180
00181
00182 function f(x,iopt)
00183 save
00184 c
00185 if (iopt.eq.1.or.iopt.eq.11) then
00186 f = x
00187 else if (iopt.eq.2.or.iopt.eq.12) then
00188 if (x.ne.0.) then
00189 f = 1./x
00190 else
00191 f = 0.
00192 end if
00193 else if (iopt.eq.3.or.iopt.eq.13) then
00194 if (x.gt.1.e-20) then
00195 f = log(x)
00196 else
00197 f = -46.05
00198 end if
00199 end if
00200 c
00201 end
00202
00203
00204
00205 function finv(x,iopt)
00206 save
00207 c
00208 c Compute the inverse function to f above.
00209 c
00210 if (iopt.eq.1.or.iopt.eq.11) then
00211 finv = x
00212 else if (iopt.eq.2.or.iopt.eq.12) then
00213 if (x.ne.0.) then
00214 finv = 1./x
00215 else
00216 finv = 0.
00217 end if
00218 else if (iopt.eq.3.or.iopt.eq.13) then
00219 if (x.gt.-46.05) then
00220 finv = exp(x)
00221 else
00222 finv = 1.e-20
00223 end if
00224 end if
00225 c
00226 end
00227
00228
00229
00230 integer*2 function comp(x,y)
00231 save
00232 c
00233 c Comparison function for use with qsort.
00234 c
00235 if (x.gt.y) then
00236 comp = 1
00237 else if (x.lt.y) then
00238 comp = -1
00239 else
00240 comp = 0
00241 end if
00242 c
00243 end
