c . 1 . 2 . 3 . 4 . 5 . 6 . 7.. c c program SINCAL1 c c Determination of the three parameters of an electromagnetic c seismometer: free period, damping and generator constant (gain) c from a record of 4 to 9 test frequencies. c c The exact frequencies of the test signals must be specified c in the parameter file sincal.par. If they are not known, use c the program sincal2. c c ******************************************************************* c c This program assumes that you have recorded the response of a c simple seismometer (with a geophone-type transfer function) to c at least four sinusoidal test signals of equal amplitude. When c f0 is the eigenfrequency of the sensor, then the signal frequencies c 0.2*f0, 0.5*f0, f0, 2*f0, and 5*f0 are a good choice. The amplitudes c must be equal in terms of calibration coil current or ground c acceleration. If your calibration-signal generator puts out c signals corresponding ground velocity, the amplitudes must be c chosen inversely proportional to the frequency. Exact values c of the test frequencies must be communicated to the program as the c first parameter t0 of the 'ban' commands in the 'sincal.par' file. c c The records of the test signals must be contained in a single c data file. Each test signal should contain about 25 full periods, c depending on the half-width of the gaussian bandpass filters used c in the analysis. This half-width, measured in periods, is entered c as the second parameter, m, of the 'ban' command. m=3 is generally c sufficient; higher values result in more selective filters, thus c a better immunity to noise. Each test frequency must be present c for at least 5*m periods at constant amplitude (more if tapered). c amplitude (more if tapered). Between the test signals, quiet c intervals corresponding to m periods of the slower one of the c adjacent signals should be present. c c The program is controlled by a parameter file 'sincal,par' c unless a different file name is specified in the runstring. c This file contains file names and processing commands. (see below). c Several data files can be processed in one run if the test c frequencies are the same. c c c ******************************************************************* c c c Derived from program SEIFE written by e. wielandt, eth zurich, 1984 c sinewave amplitudes inverted for seismometer parameters Oct. 2002 c Updated July 2004. c c this program has limited time-domain signal processing capabilities c like SEIFE, but is not a general data-processing tool, and produces c no output signal. c c CAUTION! Since this program is a modification of SEIFE, some comments c may apply to the original program rather than the modified one. I have c however tried to update the comments as well. c c convention for the 'seife' type file format: c - one header line, up to 72 characters c - up to 48 lines of signal description and processing log c these lines have '%' as the first character c (the processing log is updated by the present program) c - one line specifying the number of samples, the format, the sampling c interval, and the start time in min and sec after midnight. this c line may be repeated after the data when its contents is not known c at the time when the file header is written; however the format c must always be specified. c - data in any FORTRAN format or in free ASCII format c Other supported formats are: c - single-column ASCII with header lines ignored (read only) c - MARS88 - ASCII with header line (read only) c c the desired time-domain processes are defined in a command file whose c name is entered either interactively or as a runstring. this file also c specifies the names of the input files (up to 24; all c processes apply to all input files). each line of the command file c consists of a thee-letter code followed by one or more blanks and some c numerical parameters (such as corner period, damping, window limits). c parameters must be separated by blanks or commas and fit into columns c 5 to 50. c c the following codes are accepted: c c blank first character, or rem: line is ignored c c commands for input and output: c c fil nam1: specify input files c fil will read 'seife' files. c mar nam1 like fil, but for MARS88-ASCII files from MARSDUMP c asc nam1 dt lin where dt is the sampling interval and lin the c number of header lines to be skipped. Will read any ascii data. c The last character of the last line must be a , or the last c sample will be skipped. c c lim n: limit the length of the input series to n samples c c end: end of input file, start the inversion for system parameters c c Commands for signal processing: c c ban t0 m: gaussian FIR bandpass-filter with center period t0 c and a total width of 2*m periods. The maximum amplitude of c the output signal is stored but the output signal itself c is discarded. c avg n: remove average (determined from first n samples) c dec n: decimate to every n'th sample c dif differentiate the signal c exp n: remove exponential trend determined from first n samples c ext determine extremal samples (signal is unchanged) c fac f: multiply each sample by f c int integrate the signal c lin remove linear trend through first and last sample c nor a: normalize amplitude to a c nop a: normalize peak amplitude to a c nox a: normalize extrema to +-a c nul n1 n2: nulls the signal between samples n1 and n1 (inclusive) c pad n; pad with zeroes to length n c pol n: remove polynomial trend of degree n c rev reverse signal in time (useful for backward filtering) c rms determine rms amplitude (signal is unchanged) c skp n: skip n samples c spl clp,n1,n2,apex: interpolate clipped samples (value >= clp) with c cubic spline from sample n1 to n2. clp=0 is replaced by a default c value of slightly less than 2**18. if apex is specified other than c zero, interplation is done with a parabola with specified height of c apex. c tap n1 n2: signal is cosine- tapered between samples n1 and n2 c (signal is unchanged between samples 1 and n1) c tre n: remove linear trend (determined from first n samples) c twi tmin1 tsec1 tmin2 tsec2: window defined by time after midnight c swi sec1 sec2: window defined by seconds after first sample c win n1 n2: time window from sample n1 to n2 (n1 becomes first sample) c c recursive (IIR) filters: c c filters are specified by a three-character code for the type c and 1 to 4 parameters. c type: one of lp1, hp1, le1, he1, lp2, hp2, bp2, le2, he2, lpb, hpb c lp = low pass, hp= high pass, bp = band pass c parameters: corner period [,damping if second-order] c le = inverse low pass, he = inverse high pass (equalizers) c parameters: old corner period, [old damping], new corner c period, [new damping] c 1 or 2: number of poles, i.e. order of denominator polynomial c b: Butterworth, specify period and order c Example: lpb 30. 6 is a 30 sec, sixth-order Butterworth low-pass filter. c Additional type: sez, inverse second-order high-pass filter (restitution) c sez t0 h: remove seismo response with free period t0 and damping h c c end of sincal1 info. program sincal1 c the following statement is required with NAS-FortranPlus use nas_system parameter(ndim=120000,ndf=ndim/12) dimension x(ndim),f1(ndf),f2(ndf) dimension ix(ndim),fdt(32),lin(32),t0(9),resp(9) equivalence (x,ix) character typ*3,text*72,ttext*72,zeile*60,par*56,info*78 character*24 nam1(32),nam3 logical nil data fdt/32*0./, lin/32*0/ write(*,*) write(*,*) 'HELP: >> sincal1 info << or >> sincal1 -help <<' write(*,*) if(iargc().eq.0) then nam3='sincal.par' else call getarg(1,nam3) endif if(nam3(1:4).eq.'info'.or.nam3(1:5).eq.'-help') then open(7,file='sincal1.f',status='old') do i=1,199 read(7,'(a)',end=178) info write(6,'(a)') info if (info(1:8).eq.'c end of') stop enddo 178 stop endif write(6,'("reading commands from file ",a)') nam3 open(7,file=nam3,status='old') nmax=ndim nf=ndf call files(nam1,fdt,lin,nseis,zeile,nmax) close(7) do 1 ns=1,nseis n1=1 dt=fdt(ns) n=ndim+1-n1 n=min(n,nmax) call inpasc(nam1(ns),lin(ns),ttext,n,dt,ttmin,ttsec,x(n1)) if(n1.eq.1) then tmin=ttmin tsec=ttsec text=ttext endif if(n.eq.0.or.dt.eq.0.) then write(6,'("n or dt unspecified - check data header")') stop endif open(7,file=nam3,status='old') 2 continue read(7,'(a)',end=4) zeile if(zeile(1:1).eq.' ') goto 2 typ=zeile(1:3) par=zeile(5:60) if(typ.eq.'fil'.or.typ.eq.'asc'.or.typ.eq.'lim' & .or.typ.eq.'mar') goto 2 if(typ.eq.'prt') then write(6,*) par(1:56) goto 2 endif if(typ.eq.'rem') goto 2 nil=.true. n=min(n,nmax) if(typ.eq.'sei') then ic=0 read(par,*) ts,hs nil=.false. endif if(typ.eq.'fac') call factor(nil,par,x,n) if(typ.eq.'nor'.or.typ.eq.'nop'.or.typ.eq.'nox') & call norm(nil,typ,par,x,n) if(typ.eq.'ban') call band(nil,par,x,n,dt,f1,f2,nf,t0,resp,ic) if(typ.eq.'dif') call deriv(nil,par,x,n,dt) if(typ.eq.'spl') call interp(nil,par,x,n) if(typ.eq.'win'.or.typ.eq.'skp'.or.typ.eq.'tap'.or.typ.eq.'swi') & call window(nil,typ,par,x,n,dt,tmin,tsec) if(typ.eq.'twi') call timwin(nil,par,x,n,dt,tmin,tsec) if(typ.eq.'nul') call zero(nil,par,x,n) if(typ.eq.'dec') call decim(nil,par,x,n,dt) if(typ.eq.'avg') call mean(nil,par,x,n) if(typ.eq.'tre'.or.typ.eq.'lin') call trend(nil,typ,par,x,n) if(typ.eq.'pol') call polytrend(nil,par,x,n) if(typ.eq.'exp') call expon(nil,par,x,n) if(typ.eq.'rev') call reverse(nil,par,x,n) if(typ.eq.'ext') call extrem(nil,x,n) if(typ.eq.'rms') call rms(nil,x,n) if(typ.eq.'pad') call pad(nil,par,x,n) if(nil) call filter(nil,zeile,x,n,dt,typ,par) c all done - store result if(typ.eq.'end') then call calib(par,ts,hs,t0,resp,ic) goto 1 endif c produce error message if code unrecognized if(nil) then write(6,3) typ 3 format('prozedur ',a3,' ist nicht definiert --- stop') stop else goto 2 endif 1 continue stop 4 write(6,'("end-zeile fehlt im befehlsfile --- stop")') stop end c subroutine calib(par,ts,hs,t00,resp0,ic0) dimension t00(9),t0(9),resp0(9),resp(9),tresp(9) dimension rest0(9),rest1(9),rest2(9),tb(9),hb(9),gb(9) character par*56 logical v c invert data (ib=0) and then perform bootstrap test 1 format("Checkpoint ",i2) do ib=0,ic0 if(ib.eq.0) then v=.true. ic=ic0 do i=1,ic t0(i)=t00(i) resp(i)=resp0(i) enddo else v=.false. if(ic0.lt.4) then write(*,*) write(*,*) '*** insufficient number of test frequencies ***' stop endif if(ib.eq.1) then write(6,*) write(6,'(a)') ' Bootstrap test:' endif ic=ic0-1 iadd=0 do i=1,ic if(i.eq.ib) iadd=1 t0(i)=t00(i+iadd) resp(i)=resp0(i+iadd) enddo endif if(v) write(6,*) if(v) write(6,'(a)') ' Fitting a resonance curve to the data:' do iter=1,3 gain=0. do i=1,ic tau=t0(i)/ts tresp(i)=sqrt((tau-1./tau)**2+4.*hs**2) gain=gain+resp(i)*tresp(i) enddo gain=gain/ic err=0. do i=1,ic tresp(i)=gain/tresp(i) rest0(i)=resp(i)-tresp(i) err=err+rest0(i)**2 enddo err=sqrt(err/ic) if(v.and.iter.eq.1) write(6,111) " ts=",ts," hs=",hs, & " gain=",gain," rms=",err 111 format(a,f8.3,a,f7.4,a,f10.3,a,f8.3) gain=0. do i=1,ic tau=t0(i)/ts/1.001 tresp(i)=sqrt((tau-1./tau)**2+4.*hs**2) gain=gain+resp(i)*tresp(i) enddo gain=gain/ic do i=1,ic tresp(i)=gain/tresp(i) rest1(i)=resp(i)-tresp(i)-rest0(i) enddo gain=0. do i=1,ic tau=t0(i)/ts tresp(i)=sqrt((tau-1./tau)**2+4.008*hs**2) gain=gain+resp(i)*tresp(i) enddo gain=gain/ic err=0. do i=1,ic tresp(i)=gain/tresp(i) rest2(i)=resp(i)-tresp(i)-rest0(i) enddo c11=0. c12=0. c22=0. b1=0. b2=0. do i=1,ic c11=c11+rest1(i)*rest1(i) c12=c12+rest1(i)*rest2(i) c22=c22+rest2(i)*rest2(i) b1=b1+rest1(i)*rest0(i) b2=b2+rest2(i)*rest0(i) enddo c21=c12 det=c11*c22-c21*c12 x1=(b1*c22-b2*c12)/det x2=(c11*b2-c21*b1)/det ts=ts*(1-0.001*x1) hs=hs*(1-0.001*x2) gain=0. do i=1,ic tau=t0(i)/ts tresp(i)=sqrt((tau-1./tau)**2+4.*hs**2) gain=gain+resp(i)*tresp(i) enddo gain=gain/ic err=0. do i=1,ic tresp(i)=gain/tresp(i) rest0(i)=resp(i)-tresp(i) err=err+rest0(i)**2 enddo err=sqrt(err/ic) if(v) write(6,111) " ts=",ts," hs=",hs," gain=",gain, & " rms=",err enddo if(v) then ts0=ts hs0=hs gs0=gain write(6,*) write(6,'(a)') ' Final fitting errors:' do i=1,ic write(6,112) " t0=",t0(i)," resp=",resp(i)," theor=", & tresp(i)," err=",resp(i)-tresp(i) 112 format(a,f10.3,a,f10.3,a,f10.3,a,f8.3) enddo endif if(.not.v) then write(6,113) " w/o t0=",t00(ib),": ts=",ts, & " hs=",hs," gain=",gain," rms=",err 113 format(a,f8.3,a,f9.3,a,f7.4,a,f10.3,a,f8.3) tb(ib)=ts hb(ib)=hs gb(ib)=gain endif enddo c evaluate bootstrap test tba=0. hba=0. gba=0. tbaq=0. hbaq=0. gbaq=0. do i=1,ic0 tba=tba+tb(i) hba=hba+hb(i) gba=gba+gb(i) enddo ts1=tba/ic0 hs1=hba/ic0 gs1=gba/ic0 trms=0. hrms=0. grms=0. do i=1,ic0 trms=trms+(tb(i)-ts1)**2 hrms=hrms+(hb(i)-hs1)**2 grms=grms+(gb(i)-gs1)**2 enddo ic1=ic0-1 trms=sqrt(trms/ic1) hrms=sqrt(hrms/ic1) grms=sqrt(grms/ic1) write(*,*) write(*,*) 'RESULT: best fit boots.avg. std.dev.' write(*,114) ' Period ',ts0,ts1,' +-',trms write(*,115) ' Damping ',hs0,hs1,' +-',hrms write(*,116) ' Gain ',gs0,gs1,' +-',grms 114 format(a,f10.3,f12.3,a,f8.3) 115 format(a,f10.4,f12.4,a,f8.4) 116 format(a,f10.2,f12.2,a,f8.2,' normalized with nox?') return end c subroutine inpasc(name,lin,text,n,dt,tmin,tsec,x) dimension x(n) character iform*20,name*24,text*72,pro*80 write(6,*) '================================================' write(6,*) 'Opening file ',name open(7,file=name,status='old') c ASCII files with header namlen=index(name,' ') if(dt.gt.0.) then write(6,22) name(1:namlen),0,n,dt do i=1,lin read(7,'(a)') text enddo text=name c NAS FortranPlus read(7,*,err=30,end=30) (x(j),j=1,n) c Other compilers c read(7,*,err=31,end=31) (x(j),j=j,n) c FortranPlus does not recognize the end of data and fills c the series with zeros. Therefore, we remove trailing zeros. 30 continue do j=n,1,-1 if(x(j).ne.0.) exit enddo if(j.eq.0) then write(6,*) 'All samples are zero - stop' stop endif n=j write(6,32) n close(7) return 31 continue n=j write(6,32) n 32 format(' end of data after',i8,' samples') close(7) return c SEIFE format else if(dt.eq.0.) then read(7,'(a)') text write(6,'(1x,a)') text 21 read(7,'(a)') pro if(pro(1:1).ne.'%'.and.pro(1:1).ne.'}') goto 20 goto 21 20 read(pro,1) nn,iform,dt,tmin,tsec 1 format(i10,a20,3f10.3) write(6,22) name(1:namlen),nn,n,dt 22 format(' reading file ',a,' ',i7,'/',i7,' samples',' dt=', & f10.6) if(n.lt.nn.or.nn.eq.0) then c write(6,'("number of samples limited to",i7)') n nn=n endif n=nn if(iform(1:3).eq.'fre') then read(7,*,err=25,end=23) (x(j),j=1,n) else read(7,iform,end=23) (x(j),j=1,n) endif c MARS88-ASCII else if(dt.lt.0.) then read(7,'(a)') text text(61:72)=' ' read(text(15:17),*) idt dt=idt/1000. read(text(53:54),*) thr read(text(56:57),*) tmin read(text(59:60),*) tsec tmin=tmin+60.*thr write(6,22) name,n,n,dt do j=1,n read(7,*,err=23,end=23) x(j) enddo endif close(7) return 23 n=j-1 write(6,24) n 24 format(' end of file after ',i8,' samples') close(7) return 25 write(6,26) name,j 26 format(' Input error reading file ',a,' at sample # ',i7) stop end c subroutine gauss(aik,m,n,rs,f) c solve linear equations implicit real*8 (a-h,o-z) dimension aik(n,n),rs(n),f(n),h(14),imax(13) do 1401 j=1,m aikmax=0.d0 do 1402 k=1,m h(k)=aik(j,k) if(abs(h(k)).le.aikmax) go to 1402 aikmax=abs(h(k)) index=k 1402 continue h(m+1)=rs(j) do 1403 k=1,m q=aik(k,index)/h(index) do 1404 l=1,m 1404 aik(k,l)=aik(k,l)-q*h(l) 1403 rs(k)=rs(k)-q*h(m+1) do 1405 k=1,m 1405 aik(j,k)=h(k) rs(j)=h(m+1) 1401 imax(j)=index do 1406 j=1,m index=imax(j) 1406 f(index)=rs(j)/aik(j,index) return end subroutine rekfl(x,y,n,f0,f1,f2,g1,g2) c perform recursive filtering implicit real*8 (a-h,o-z) real*4 x(n),y(n) xa=0.d0 xaa=0.d0 ya=0.d0 yaa=0.d0 do 1 j=1,n xn=dble(x(j)) yn=f0*xn+f1*xa+f2*xaa+g1*ya+g2*yaa y(j)=sngl(yn) xaa=xa xa=xn yaa=ya 1 ya=yn return end subroutine rfk(it,t0,h,t0s,hs,f0,f1,f2,g1,g2) c determine coefficients for recursive filter implicit real*8 (a-h,o-z) data zero,one,two,four,eight/0.d0,1.d0,2.d0,4.d0,8.d0/ if(it.ne.1) goto 10 f0=one/two/t0 f1=f0 f2=zero g1=one g2=zero return 10 zpi=eight*datan(one) eps=zpi/t0 f2=zero g2=zero if(it.gt.20) goto 20 g1=(two-eps)/(two+eps) if(it.gt.11) goto 12 f0=eps/(two+eps) f1=f0 goto 14 12 if(it.gt.12) goto 13 f0=two/(two+eps) f1=-f0 goto 14 13 if(it.gt.13) return epss=zpi/t0s f0=(epss+two)/(eps+two) f1=(epss-two)/(eps+two) 14 return 20 epsq=eps*eps a=one-eps*h+epsq/four b=-two+epsq/two c=one+eps*h+epsq/four g1=-b/c g2=-a/c if(it.gt.21) goto 22 f0=epsq/four/c f1=f0+f0 f2=f0 goto 25 22 if(it.gt.22) goto 23 f0=one/c f1=-f0-f0 f2=f0 goto 25 23 if(it.gt.23) goto 24 epss=zpi/t0s epssq=epss*epss as=one-epss*hs+epssq/four bs=-two+epssq/two cs=one+epss*hs+epssq/four f0=cs/c f1=bs/c f2=as/c goto 25 24 if(it.gt.24) return f0=eps/two/c f1=zero f2=-f0 25 return end subroutine files(nam1,fdt,lin,nseis,zeile,npts) character nam1(32)*24,zeile*60,z3*3 dimension fdt(32),lin(32) nseis=0 101 read(7,'(a)',end=2) zeile write(6,'(a)') zeile(1:60) z3=zeile(1:3) if(z3.eq.'fil'.or.z3.eq.'asc'.or.z3.eq.'mar') then nseis=nseis+1 if(z3.eq.'mar') then fdt(nseis)=-1. lin(nseis)=1 endif do 30 i=5,60 30 if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 31 31 ii=i do 11 i=ii,60 11 if (zeile(i:i).eq.' '.or.zeile(i:i).eq.',') goto 34 34 if(i.gt.ii+24) goto 24 nam1(nseis)=zeile(ii:i-1) if(z3.eq.'asc') then ii=i do 14 i=ii,60 14 if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 4 4 if(i.gt.58) goto 26 if(z3.eq.'asc') then read(zeile(i:60),*,err=26) fdt(nseis),lin(nseis) else read(zeile(i:60),*,err=26) lin(nseis) lin(nseis)=-lin(nseis) endif endif else if(z3.eq.'lim') then read(zeile,'(4x,i10)') npts else if(z3.eq.'end') then if(nseis.gt.0) then return else write(6,3) 3 format(' Kein Filename angegeben - stop') stop endif endif goto 101 2 write(6,'("end-zeile fehlt im befehlsfile --- stop")') stop 24 write(6,25) zeile 25 format(' Filename zu lang: ',a) stop 26 write(6,27) 27 format(' Anweisung asc nicht lesbar - stop') stop end subroutine factor(nil,par,x,n) c multiply by a constant factor dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) f do 1 j=1,n 1 x(j)=f*x(j) write(6,'(" fac ",f10.3)') f nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for fac' stop end subroutine norm(nil,typ,par,x,n) c normalize rms, peak, or peak-to-peak amplitude dimension x(n) character typ*3,par*56 logical nil read(par,*,err=99,end=99) fn xrms=0. xmin=1e24 xmax=-xmin do 1 j=1,n xrms=xrms+x(j)**2 xmin=min(xmin,x(j)) 1 xmax=max(xmax,x(j)) xrms=sqrt(xrms/n) xp=max(abs(xmin),abs(xmax)) if(typ.eq.'nox') then f=2.*fn/(xmax-xmin) off=(xmax+xmin)/2. do 2 j=1,n 2 x(j)=f*(x(j)-off) else if(typ.eq.'nor') then f=fn/xrms else f=fn/xp endif do 3 j=1,n 3 x(j)=f*x(j) endif write(6,'(1x,a,a,e12.6)') typ,': normalization factor ',f nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for ',typ stop end subroutine pad(nil,par,x,n) c Nullen ergaenzen dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) nul do 1 j=n+1,nul 1 x(j)=0. n=nul write(6,2) nul 2 format(" pad: Es wurden Nullen aufgefuellt bis Punkt",i8) nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for pad' stop end subroutine zero(nil,par,x,n) c null time series from n1 to n2 dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) n1,n2 do 1 j=n1,n2 1 x(j)=0. write(6,'(a,2i10)') ' null ',n1,n2 nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for null' stop end subroutine deriv(nil,par,x,n,dt) c derivative (non-recursive, non-causal, symmetric-difference) dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) tau if(tau.eq.0.) tau=1. twodt=2.*dt/tau do 1 j=1,n-2 1 x(j)=(x(j+2)-x(j))/twodt do 2 j=n-1,2,-1 2 x(j)=x(j-1) x(n)=x(n-1) write(6,*) 'dif - derivative by symmetric differences' nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for dif' stop end subroutine interp(nil,par,x,n) c interpolate clipped samples real*8 aik(4,4),xx(4),f(4),t(4) dimension x(n) character*56 par logical nil data zero,one,two/0.d0,1.d0,2.d0/ c define illegal samples read(par,*,err=99,end=99) clip,n1,n2,apex if(clip.eq.0.) clip=262000. n1=max(n1,1) if(n2.eq.0.or.n2.gt.n) n2=n c search for illegal samples naus=0 nabt=0 10 do 1 j=n1,n2 if(abs(x(j)).gt.clip) goto 2 1 continue c if no more, done nil=.false. write(6,7) naus,nabt,clip,apex 7 format(" spl - an",i4," stellen wurden",i5," werte >",f9.0, & " interpoliert. apex =",f9.0) return 2 j1=j c determine the last bad sample (two good ones must follow) do 3 j2=j1,n-2 if(abs(x(j2+1)).lt.clip.and.abs(x(j2+2)).lt.clip) goto 4 3 continue if(j1.lt.n1+2.or.j2.gt.n2-2) then write(6,'(" *** spl *** zu wenig stuetzwerte - stop")') stop endif 4 if(apex.eq.zero) then c compute spline xx(1)=x(j1-2) t(1)=-two xx(2)=x(j1-1) t(2)=-one xx(3)=x(j2+1) t(3)=j2-j1+1 xx(4)=x(j2+2) t(4)=t(3)+one do 6 j=1,4 aik(j,1)=one do 6 k=2,4 6 aik(j,k)=t(j)*aik(j,k-1) else c compute parabola xx(1)=x(j1-1) t(1)=-one xx(3)=x(j2+1) t(3)=j2-j1+1 xx(2)=apex t(2)=(t(1)+t(3))/two do 16 j=1,3 aik(j,1)=one aik(j,4)=zero aik(4,j)=0.d0 do 16 k=2,3 16 aik(j,k)=t(j)*aik(j,k-1) aik(4,4)=one xx(4)=zero endif call gauss(aik,4,4,xx,f) c interpolate do 5 j=j1,j2 tt=j-j1 5 x(j)=f(1)+tt*(f(2)+tt*(f(3)+tt*f(4))) n1=j2+1 naus=naus+1 nabt=nabt+j2-j1+1 goto 10 99 write(6,*) ' insufficient or wrong parameters for spl' stop end subroutine window(nil,typ,par,x,n,dt,tmin,tsec) c time window by sample number dimension x(n) character*56 par character*3 typ logical nil if(typ.eq.'skp') then read(par,*,err=99,end=99) np1 write(6,*) ' ',typ,np1 else if(typ.eq.'swi') then read(par,*,err=99,end=99) sec1,sec2 np1=max(nint(sec1/dt)+1,1) np2=min(nint(sec2/dt),n) write(6,*) ' ',typ,(np1-0.5)*dt,(np2-0.5)*dt np1=np1-1 n=np2-np1 else read(par,*,err=99,end=99) np1,np2 np1=max(np1,1) np2=min(np2,n) write(6,*) ' ',typ,np1,np2 np1=np1-1 n=np2-np1 endif if(typ.ne.'tap') goto 16 n=np2 pih=2.*atan(1.) fak=pih/(np2-np1) do 14 j=np1,np2 14 x(j)=x(j)*cos(fak*(j-np1))**2 goto 15 16 if(typ.ne.'cos') goto 13 n1=n+1 pih=2.*atan(1.) fak=pih/np1 do 17 j=1,np1 taper=sin(fak*j)**2 x(j)=x(j)*taper 17 x(n1-j)=x(n1-j)*taper goto 15 13 if(typ.eq.'win'.or.typ.eq.'swi') then do 11 j=1,n 11 x(j)=x(np1+j) endif if(typ.eq.'skp') then n=n-np1 do j=1,n x(j)=x(np1+j) enddo endif tsec=60.*tmin+tsec+np1*dt tmin=int(tsec/60.) tsec=tsec-60.*tmin 15 nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for ',typ stop end subroutine timwin(nil,par,x,n,dt,tmin,tsec) c time window by time after midnight dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) tmin1,tsec1,tmin2,tsec2 t0=60.*tmin+tsec t1=60.*tmin1+tsec1 t2=60.*tmin2+tsec2 np1=1.5+(t1-t0)/dt np2=0.5+(t2-t0)/dt if(np1.lt.1.or.np2.gt.n.or.np2.lt.np1) then write(6,1) tmin1,tsec1,tmin2,tsec2 1 format(" twi ",4f10.1/" *** impossible window ***") if(np1.lt.1) then tsec1=tsec1+(1-np1)*dt np1=1 endif if(np2.gt.n) then tsec2=tsec2-(np2-n)*dt np2=n endif if(np2.lt.np1) stop write(6,'(" new window from sample ",i7," to",i7)') np1,np2 endif write(6,2) tmin1,tsec1,tmin2,tsec2 2 format(' twi ',4f10.1) np1=np1-1 n=np2-np1 do 3 j=1,n 3 x(j)=x(np1+j) tsec=t0+np1*dt tmin=int(tsec/60.) tsec=tsec-60.*tmin nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for twi' stop end subroutine decim(nil,par,x,n,dt) c decimate time series by a factor up to 100 dimension x(n),g(100) character*56 par logical nil read(par,*,err=99,end=99) idez if(idez.gt.100) stop 100 pih=2.*atan(1.) fidez=idez do 3 j=1,idez 3 g(j)=(cos((j-idez)/fidez*pih))**2/fidez n=n/idez dt=dt*idez do 5 k=1,n-1 jk=k*idez+1 x(k)=g(idez)*x(jk) do 6 l=1,idez-1 6 x(k)=x(k)+g(l)*(x(jk-idez+l)+x(jk+idez-l)) 5 continue do 7 j=n,2,-1 7 x(j)=x(j-1) write(6,'(" dec ",i10)') idez nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for dec' stop end subroutine mean(nil,par,x,n) c remove average character*56 par dimension x(n) logical nil read(par,*,err=99,end=99) n2 if(n2.lt.1.or.n2.gt.n) n2=n sum=0. do 1 j=1,n2 1 sum=sum+x(j) sum=sum/n2 do 2 j=1,n 2 x(j)=x(j)-sum write(6,'(" avg ",e11.3)') sum nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for avg' stop end subroutine trend(nil,typ,par,x,n) c remove trend implicit real*8 (a-h,o-z) real*4 x(n) logical nil character typ*3,par*56 common /tre/ a,b if(typ.eq.'tre') then read(par,*,err=99,end=99) n2 if(n2.le.1.or.n2.gt.n) n2=n gn = n2 alpha = 0.5d0*gn*(gn+1.d0) beta = (2.d0*gn+1.d0)*(gn+1.d0)*gn/6.d0 det = gn*beta-alpha*alpha sx = 0.d0 sjx = 0.d0 do 1001 j=1,n2 sx = sx+x(j) 1001 sjx = sjx+x(j)*j a = (sx*beta-sjx*alpha)/det b = (sjx*gn-sx*alpha)/det else if(typ.eq.'zsp') then n2=n gn = n alpha = 0.5d0*gn*(gn+1.d0) beta = (2.d0*gn+1.d0)*(gn+1.d0)*gn/6.d0 det = gn*beta-alpha*alpha sx = 0.d0 sjx = 0.d0 do 1003 j=1,n2 sx = sx+x(j) 1003 sjx = sjx+x(j)*j a = (sx*beta-sjx*alpha)/det b = (sjx*gn-sx*alpha)/det return else b=(x(n)-x(1))/float(n-1) a=x(1)-b endif do 1002 j=1,n 1002 x(j) = x(j)-sngl(a+b*j) write(6,'(" tre ",2f10.3)') a,b nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for ',typ stop end subroutine polytrend(nil,par,x,n) c remove polynomial trend parameter(ndi=6) implicit real*8 (a-h,o-z) real*4 x(n) dimension b(ndi),c(ndi,ndi),a(ndi) logical nil character*56 par read(par,*,err=99,end=99) m m=min(m,ndi-1) fnh=n/2. one=1.d0 do j=1, m+1 do k=1, m+1 c(j,k)=0 do i=1, n c(j,k)=c(j,k)+(dble(i)/fnh-one)**(j+k-2) enddo enddo b(j)=0 do i=1,n b(j)=b(j)+(dble(i)/fnh-one)**(j-1)*x(i) enddo enddo call gauss(c,m+1,ndi,b,a) write(6,100) (j-1,a(j),j=1,m+1) 100 format(i5,e15.6) do i=1,n xpol=a(m+1) do j=m,1,-1 xpol=xpol*(dble(i)/fnh-one)+a(j) enddo x(i)=x(i)-sngl(xpol) enddo write(6,'(" pol ",i5)') m nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for pol' stop end subroutine expon(nil,par,x,n) c remove exponential trend implicit real*8 (a-h,o-z) real*4 x(n) logical nil character*56 par dimension sum(3) read(par,*,err=99,end=99) nex if(nex.lt.3.or.nex.gt.n) nex=n n3=nex/3 do 1 i=1,3 sum(i)=0.d0 ii=(i-1)*n3 do 1 j=1,n3 1 sum(i)=sum(i)+x(j+ii) q=(sum(2)-sum(3))/(sum(1)-sum(2)) c=-dlog(q)/dble(n3) b=c*(sum(1)-sum(2))/(1.d0-q)**2 a=(sum(1)+sum(2)+sum(3)-b/c*(1.d0-q**3))/dble(n3)/3.d0 do 2 j=1,n 2 x(j) = x(j)-sngl(a+b*dexp(-c*dble(j))) write(6,'(" exp ",2f10.3,f10.6)') a,b,c nil=.false. return 99 write(6,*) ' insufficient or wrong parameters for exp' stop end subroutine reverse(nil,par,x,n) c reverse signal in time logical nil dimension x(n) character*56 par do j=1,n/2 jj=n+1-j xx=x(j) x(j)=x(jj) x(jj)=xx enddo write(6,'(" rev - signal was reversed")') nil=.false. return end c subroutine filter(nil,zeile,x,n,dts,typ,par) implicit real*8 (a-h,o-z) character typ*3,typ1*3,typ2*3,par*56,zeile*60 real*4 x(n),dts logical nil dt=dts typ=zeile(1:3) par=zeile(5:60) c Subtract value of first sample before high-pass filtration if(typ(1:1).eq.'h') then x0=x(1) do i=1,n x(i)=x(i)-x0 enddo endif c decode filter type and read parameters if(typ.eq.'lp1') then read(par,*,err=99,end=99) t0 it=11 else if(typ.eq.'hp1') then read(par,*,err=99,end=99) t0 it=12 else if(typ.eq.'lp2') then read(par,*,err=99,end=99) t0,h it=21 else if(typ.eq.'hp2') then read(par,*,err=99,end=99) t0,h it=22 else if(typ.eq.'bp2') then read(par,*,err=99,end=99) t0,h it=24 else if(typ.eq.'int') then t0=1.d0 it= 1 else if(typ.eq.'he1'.or.typ.eq.'le1') then read(par,*,err=99,end=99) t0s,t0 it=13 else if(typ.eq.'he2'.or.typ.eq.'le2') then read(par,*,err=99,end=99) t0s,hs,t0,h it=23 else if(typ.eq.'sez') then read(par,*,err=99,end=99) t0s,hs t0=1.e12 h=1. it=23 else goto 2 endif c calculate filter weights call rfk(it,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2) if(typ.eq.'le1') then fac=t0s/t0 f0=f0*fac f1=f1*fac endif if(typ.eq.'le2') then fac=(t0s/t0)**2 f0=f0*fac f1=f1*fac f2=f2*fac endif c perform recursive filtration call rekfl(x,x,n,f0,f1,f2,g1,g2) 8 format(1x,a3,2x,5f10.3) nil=.false. c confirm execution if(it.eq.1) write(6,8) typ,t0 if(it.eq.11.or.it.eq.12) write(6,8) typ,t0 if(it.eq.13) write(6,9) typ,t0s,t0 9 format(1x,a3,2x,f10.3,f16.3) if(it.eq.21.or.it.eq.22.or.it.eq.24) write(6,8) typ,t0,h if(it.eq.23) then if(typ.eq.'sez') then write(6,8) typ,t0s,hs else write(6,7) typ,t0s,hs,t0,h 7 format(1x,a3,2x,2f10.3,f16.3,f10.3) endif endif return c Butterworth filters 2 if(typ.eq.'lpb'.or.typ.eq.'hpb') then read(par,*,err=99,end=99) t0,m mm=m/2 if(typ.eq.'lpb') then it1=11 it2=21 typ1='lp1' typ2='lp2' else it1=12 it2=22 typ1='hp1' typ2='hp2' endif if(m.gt.2*mm) then call rfk(it1,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2) call rekfl(x,x,n,f0,f1,f2,g1,g2) write(6,8) typ1,t0 endif pih=2.d0*datan(1.d0) do 3 j=1,mm h=dsin(pih*(2*j-1)/m) call rfk(it2,t0/dt,h,t0s/dt,hs,f0,f1,f2,g1,g2) call rekfl(x,x,n,f0,f1,f2,g1,g2) write(6,8) typ2,t0,h 3 continue nil=.false. endif return 99 write(6,*) ' insufficient or wrong parameters for ',typ stop end subroutine extrem(nil,x,n) dimension x(n) logical nil nil=.false. xmin=1e24 xmax=-xmin do 1 j=1,n if(x(j).lt.xmax) goto 3 xmax=x(j) nmax=j 3 if(x(j).gt.xmin) goto 1 xmin=x(j) nmin=j 1 continue write(6,2) xmin,nmin,xmax,nmax 2 format(" ext xmin=",e11.3," (j=",i6,") xmax=",e11.3," (j=",i6, &")") return end subroutine rms(nil,x,n) dimension x(n) logical nil nil=.false. xrms=0. do 1 j=1,n 1 xrms=xrms+x(j)**2 write(6,2) sqrt(xrms/n) 2 format(" rms =",e11.3) return end subroutine band(nil,par,x,n,dt,f1,f2,nf,t0,resp,ic) character*56 par logical nil dimension x(n),f1(nf),f2(nf),t0(9),resp(9) ic=ic+1 if(ic.eq.1) write(6,'(/,a)') & ' Bandpass filtration and amplitude measurement:' if(ic.gt.9) then write(*,*) "This program accepts only 9 test frequencies" stop endif read(par,*,err=99,end=99) t0(ic),m lenh=int(m*t0(ic)/dt) if(lenh.gt.nf) then write(6,*) 'ban: FIR filter is too long - stop' stop endif twopi=8.*atan(1.) om=twopi/t0(ic) fn1=1. fn2=0. do k=1,lenh co=cos(om*k*dt) si=sin(om*k*dt) f1(k)=co*exp(-(2.5*k/lenh)**2) f2(k)=si*exp(-(2.5*k/lenh)**2) fn1=fn1+2.*f1(k)*co fn2=fn2+2.*f2(k)*si enddo do k=1,lenh f1(k)=f1(k)/fn1 f2(k)=f2(k)/fn2 enddo xmax=0. do i=1,n yy=x(i)/fn1 zz=0. do k=1,lenh if(i.gt.k) then yy=yy+f1(k)*x(i-k) zz=zz-f2(k)*x(i-k) endif if(i+k.le.n) then yy=yy+f1(k)*x(i+k) zz=zz+f2(k)*x(i+k) endif enddo xmax=max(xmax,yy**2+zz**2) enddo xmax=sqrt(xmax) resp(ic)=xmax nil=.false. write(6,90) " filter t0=",t0(ic)," m=",m," resp=",resp(ic) 90 format(a,f8.3,a,i2,a,f10.3) return 99 write(6,*) ' insufficient or wrong parameters for ban' stop end