c program SEIFE c c REVISIONS (there were many; this protocol starts November 2012) c Nov. 2012: the global 'frq' command permits specification of corner c frequencies in place of corner periods c prewarping activated for butterworth filters c July 2013: "nul" parameters limited to 1..n c Sept. 2013: intermediate results can be stored with 'sav' c Oct. 2014: output format 'asc' implemented (MATLAB compatible) c c time-domain signal processing (filtering, windowing, decimation, etc.) c Author: E. Wielandt c The original version was written in 1984 at ETH Zurich 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 The following INPUT formats are also supported: c - ASCII with a variable number of samples per line, header lines ignored (asc) c - single-column ASCII, header lines ignored (sca) c - MARS88 - ASCII with header line c - ASL (Peterson) ASCII format such as produced by Quanterra's CIMARRON c - BDF ASCII as used in the Albuquerque Seismological Lab c - PSN textfiles from WINQUAKE (Larry Cochrane) c c OUTPUT can be either in SEIFE format or in one-column or two-column headerless c ASCII format compatible with GNUPLOT, MATLAB and LABVIEW 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 and output files (up to 24 pairs; all c processes apply to all input files). each line of the command file c defines one processing step (a time window, filter, etc.) and 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 c columns 5 to 60. The following codes are accepted: c c blank first character, or rem: line is ignored c c commands for input and output: c c Reading data files: c fil nam1 nam2: read files in SEIFE format. nam1=input file, nam2=output file c nam2='_join_' will append the next input file to the present one c mar nam1 nam2: read MARS88-ASCII files from MARSDUMP c asl nam1 nam2: read ASL ascii format (from Cimarron) c bdf nam1 nam2: read BDF ascii (headers interspersed with data - caution, c anything that is not a number will be ignored) c psn nam1 nam2: read PSN text files from WINQUAKE (Larry Cochrane) c asc nam1 nam2 dt lin: read ASCII files where dt is the sampling interval and c lin the number of header lines to be skipped. (Header lines are echoed c but not interpreted.) Will read almost any ASCII data. c The last character of the last line must be a , or the last c sample will be skipped. asc is slow; use sca if possible. c sca (single-column ascii) like asc but reads only one sample per line (much c faster than asc) c c global commands: c c frq specify filter corners as frequencies (default is periods) c c chp nch: chop signal into sections of nch samples for output. c Append .00, .01 etc. to name of output file (maximum 24 files). c c com insert comment into file header c rem insert comment into file header c c lim n: limit the length of the input series to n samples c c ncl no comment lines c c npl don't write a plot-parameter file c c pfi set a path for the input files. Don't use quotation marks. c The path may also be entered as part of the file name. c c pfo set a path for the output files. Don't use quotation marks. c The path may also be entered as part of the file name. c c pli set up winplot.par to plot input and output files. Works for c SEIFE format only, and for up to six input signals. Otherwise ignored. c Do not apply a time window together with pli (traces will not be c synchronous if you do). c c plt w h: set nonstandard plot width w and height h in winplot.par c c plw n1 n2: set window for plot, rms, and ext. Set 'remove average' for plot. c c cfi specify name for the Labview cfg file. Default is labview.cfg c c sav save intermediate result; attach -1, -2 etc. to name of output file c c end [format]: write output file. c The output file may be identical to the input file. c Legal format specifications are: c (any fortran format) - specified format, in parentheses but without c quotation marks, must fit into 20 characters. Example:(5x,10f8.3) c fre or = 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 sqr square each sample (used for the detection of quadratic distortions) c srt heapsort c swi sec1 sec2: window defined by seconds after first sample c zsp n1 n2 n3 n4 detrend with a z-spline with corners n1 .. n4 c tap n1: signal is tapered from the first sample to sample n1 and c from sample (last-n1) to last. The weight is a cosine half-cycle. c The amplitude is normalized so that the energy of a stationary c signal is unchanged. c tid [n]: remove tides. tides are interpolated over n samples. c default (n=0 or blank) is interpolation over 5 minutes. c tr3 n1 t0 h: fit and remove remove transient with eigenperiod t0 c and damping h, after sample n1 c tr6 n1 t0 h f1 f2 f3: remove transient using amplitudes determined c by tr3 (possibly from a different signal) c tre n: remove linear trend (determined from first n samples) c twi tmin1 tsec1 tmin2 tsec2: window defined by time after midnight c tim tmin tsec set time stamp 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, sez 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]. These are APPROXIMATE inverses. 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 sez t0 h is the exact inverse to hp2 t0 h c end of seife info. program seife c the following statement is compiler-specific, it makes runstring c parameters available use nas_system parameter(ndim=200000000) dimension x(ndim) dimension fdt(32),lin(32) real*8 ddt,tcom1,tcom2,tcom1m,tcom1s,tcom2m,tcom2s character typ*3,text*72,ttext*72,zeile*99,par*56,info*78 character*60 prolog(48) character*48 nam1(32),nam2(32),nam3,nam4,cfgname,inpath,outpath character*256 cfgline(32) character timestring*12,einsbisneun*2 character*3 datform(32) logical nil,join,bin,bina,ncl,npl,npw,pli,frq data fdt/32*0./, lin/32*0/, breit/24./, hoch/16./ data cfgname,inpath,outpath /3*'#'/ data ncl,npl,npw,pli,frq /5*.false./ write(*,*) write(*,*) 'Latest revision: Feb. 2013. For a list of options,' write(*,*) 'copy the source code into the working directory and' write(*,*) 'type "seife info". See also the program description', & ' seife.doc' write(*,*) if(iargc().eq.0) then nam3='seife.par' else call getarg(1,nam3) endif if(nam3(1:4).eq.'info') then open(7,file='seife.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)') trim(nam3) open(7,file=nam3,status='old') nchin=ndim nmax=ndim ntrend=0 call files(nam1,nam2,fdt,lin,nseis,zeile,nmax,nchin,ncl,npl,npw, & pli,frq,cfgname,inpath,outpath,datform) close(7) tcom1=-1.d12 tcom2=1.d12 lines=0 join=.false. npl1=1 npl2=nmax do 1 ns=1,nseis nsav=0 nch=nchin npro=0 if(join) then n1=n+1 else n1=1 endif bina=bin bin=.false. dt=fdt(ns) if(join.and.(bina.ne.bin)) then stop ' Cannot join binary and ascii files' endif n=ndim+1-n1 if(bin) then stop ' binary input not implemented' else if (lin(ns).ge.0) then n=min(n,nmax) call inpasc(nam1(ns),lin(ns),ttext,prolog,npro,n,dt,ttmin,ttsec, & x(n1),ndim,datform(ns)) else stop ' unknown format specification' endif npl2=min(npl2,n) tcom1=max(tcom1,60.d0*ttmin+ttsec) tcom2=min(tcom2,60.d0*ttmin+ttsec+(n-1)*dt) if(n1.eq.1) then tmin=ttmin tsec=ttsec text=ttext endif if(n.eq.0.or.dt.eq.0.) then stop ' n or dt unspecified - check data header' endif if(join) then n=n1+n-1 write(6,11) n 11 format(' data appended to previous data, total ',i8,' samples') endif join=nam2(ns).eq.'_join_' if(join) goto 1 npro=npro+1 write(prolog(npro),5) trim(nam1(ns)),trim(nam2(ns)) 5 format('from ',a,' to ',a) 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) if(zeile(4:4).ne.' ') write(6,*) 'WARNING: parameter file has non- &blank entry in column 4' par=zeile(5:60) if(typ.eq.'fil'.or.typ.eq.'asc'.or.typ.eq.'lim'.or.typ.eq.'chp' & .or.typ.eq.'asl'.or.typ.eq.'bdf'.or.typ.eq.'ncl'.or.typ.eq.'mar' & .or.typ.eq.'pfi'.or.typ.eq.'pfo'.or.typ.eq.'cfi'.or.typ.eq.'npl' & .or.typ.eq.'npw'.or.typ.eq.'pli'.or.typ.eq.'sca'.or.typ.eq.'psn' & .or.typ.eq.'frq') & goto 2 npro=npro+1 if(npro.gt.48) then stop ' prolog has more than 48 lines' endif write(prolog(npro),'(a)') zeile if(typ.eq.'ext'.or.typ.eq.'rms') npro=npro-1 if(typ.eq.'rem'.or.typ.eq.'com') goto 2 nil=.true. n=min(n,nmax) npl2=min(npl2,n) ddt=dt if(typ.eq.'sqr') call square (nil,x,n) 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.'add') call add(nil,par,x,n) if(typ.eq.'dif') call deriv(nil,par,x,n,dt) if(typ.eq.'dis') call distort(nil,par,x,n) if(typ.eq.'spl') call interp(nil,par,x,n) if(typ.eq.'shl'.or.typ.eq.'shr') call shift(nil,typ,par,x,n) if(typ.eq.'csi') call intpo (nil,par,dt,x,n) if(typ.eq.'win'.or.typ.eq.'sin'.or.typ.eq.'skp'.or.typ.eq.'tap' & .or.typ.eq.'sis'.or.typ.eq.'cos'.or.typ.eq.'swi') & call window(nil,pli,typ,par,x,n,dt,tmin,tsec) if(typ.eq.'twi') call timwin(nil,par,x,n,dt,tmin,tsec) if(typ.eq.'tim') call timeset(nil,par,tmin,tsec) if(typ.eq.'nul') call zero(nil,par,x,n) if(typ.eq.'set') call set(nil,par,x,n) if(typ.eq.'clp') call clip(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.'zsp') call zspline(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.'tid') call tides(nil,par,x,n,ddt) if(typ.eq.'tr3'.or.typ.eq.'tr6') call trans(nil,typ,par,x,n,ddt) if(typ.eq.'ext') call extrem(nil,x,npl1,npl2,n) if(typ.eq.'rms') call rms(nil,x,npl1,npl2,n) if(typ.eq.'pad') call pad(nil,par,x,n,ndim) if(typ.eq.'del') call delay(nil,par,x,n) if(typ.eq.'srt') call heapsort(nil,x,n) if(typ.eq.'sig'.or.typ.eq.'abs') call sig(nil,zeile,x,n) if(typ.eq.'plt') call plotpar(nil,par,breit,hoch) if(typ.eq.'qnt') call quant(nil,par,x,n) if(typ.eq.'plw') call plotwindow(nil,par,npl1,npl2,n,ntrend) if(nil) call filter(nil,zeile,x,n,ddt,npw,frq) c all done - store result if(typ.eq.'sav'.or.typ.eq.'end') then nch=min(max(nch,int((n+9)/24)),n) nfile=int(n-1)/nch+1 if(typ.eq.'sav'.and.nfile.gt.1) then write(6,*) 'sav ignored, it is incompatible with chp' goto 2 endif if(typ.eq.'sav'.and.nseis.gt.1) then write(6,*) 'sav ignored, it accepts only one input file' goto 2 endif do nfi=1,nfile n1=1+(nfi-1)*nch n2=min(nch,n-n1+1) nam4=nam2(ns) if(nfile.gt.1) then i=index(nam4,' ') write(nam4(i:i+2),6) nfi 6 format('.',i2.2) tsec=60.*tmin+tsec+(n1-1)*dt tmin=int(tsec/60.) tsec=tsec-60.*tmin endif if(typ.eq.'sav') then nsav=nsav+1 if(nsav.gt.9) goto 2 write(einsbisneun,'("-",i1)') nsav nam4=trim(nam2(ns))//einsbisneun endif call outform(nam4,par,text,prolog,npro,ncl,n2,dt,tmin,tsec, & x(n1),cfgline,lines) enddo if(typ.eq.'sav') goto 2 goto 1 endif c produce error message if code unrecognized if(nil) then write(*,*) ' undefined procedure ',typ stop else goto 2 endif 1 continue c print common start and end time tcom1m=int(tcom1/60.d0) tcom1s=tcom1-60.d0*tcom1m tcom2m=int(tcom2/60.d0) tcom2s=tcom2-60.d0*tcom2m write(6,20) tcom1m,tcom1s,tcom2m,tcom2s 20 format(" common time interval: from ",f6.0," min ",f7.3," sec"/ & " to ",f6.0," min ",f7.3," sec") if(lines.eq.0) then c produce plot-parameter file open(9,file='winplot.par') if(.not.npl.and.nsav.eq.0) then jcount=0 do ns=1,ns if(trim(nam2(ns)).eq.'_join_') jcount=jcount+1 enddo if(pli.and.jcount.eq.0) then write(9,7) ntrend,2*nseis,breit,hoch,npl1,npl2 do ns=1,nseis write(9,'(a)') trim(nam1(ns)) write(9,'(a)') trim(nam2(ns)) enddo else write(9,7) ntrend,nseis-jcount,breit,hoch,npl1,npl2 7 format(i3,',',i3,",",f5.1,",",f5.1,",",i8,",",i8,", 0.8") do ns=1,nseis if(trim(nam2(ns)).ne.'_join_') write(9,'(a)') trim(nam2(ns)) enddo endif else if(.not.npl.and.nsav.gt.0) then nsav=min(nsav,9) if(pli) then write(9,7) ntrend,nsav+2,breit,hoch,npl1,npl2 write(9,'(a)') trim(nam1(1)) do i=1,nsav write(einsbisneun,'("-",i1)') i write(9,'(a)') trim(nam2(1))//einsbisneun enddo write(9,'(a)') trim(nam2(1)) else write(9,7) ntrend,nsav+1,breit,hoch,npl1,npl2 do i=1,nsav write(einsbisneun,'("-",i1)') i write(9,'(a)') trim(nam2(1))//einsbisneun enddo write(9,'(a)') trim(nam2(1)) endif close(9) write(6,*) ' Plot-parameter file winplot.par was generated' else write(6,*) ' this should not happen - check program' stop endif else c produce cfg file for Labview if(index(cfgname,'#').ne.0) then cfgname='labview.cfg' endif if(index(outpath,'#').eq.0) then cfgname=trim(outpath)//'\'//cfgname endif cfgline(lines+1)='Type: Single' ithr1=int((tcom1m*60.d0+tcom1s)/3600.d0) itmin1=int((tcom1m*60.d0+tcom1s-3600.d0*ithr1)/60.d0) tcom1s=int((tcom1m*60.d0+tcom1s-3600.d0*ithr1-60.d0*itmin1)) write(timestring,77)ithr1,itmin1,tcom1s 77 format(i2,':',i2,':',f6.3) cfgline(lines+2)='Startzeit: 00/00/00 '//timestring ithr2=int((tcom2m*60.d0+tcom2s)/3600.d0) itmin2=int((tcom2m*60.d0+tcom2s-3600.d0*ithr2)/60.d0) itsec2=int((tcom2m*60.d0+tcom2s-3600.d0*ithr2-60.d0*itmin2)) write(timestring,77)ithr2,itmin2,tcom2s cfgline(lines+3)='Endzeit: 00/00/00 '//timestring open(9,file=cfgname) do i=1,lines+3 write(9,'(a)') trim(cfgline(i)) enddo close(9) write(6,*) ' configuration file ',trim(cfgname),' was generated' endif stop 4 write(*,*) ' end line missing in file ',trim(nam3) stop end c subroutine inpasc(name,lin,text,prolog,npro,n,dt,tmin,tsec,x,ndim, &datform) dimension x(ndim),y(60),ix(60) character iform*20,name*48,text*72,pro*80,code*4,rest*40 character prolog(48)*60,code1*6,code2*6,datform*3 character zeile*120,b,zarr(120),zeil*40 common /asl/ code1,code2,year,day,thr,tm,ts,t1000,srate,cpv real*8 srate8 logical first c b=' ' cpv=1.000 do j=1,n x(j)=0. enddo write(6,*) '===' write(6,*) ' Opening file ',name open(7,file=name,status='old') c read ASCII files if(datform.eq.'asc'.or.datform.eq.'sca') then write(6,22) trim(name),0,ndim,dt do i=1,lin read(7,'(a)') text write(6,*) ' header line: ',trim(text) prolog(i)=text enddo npro=lin text=name if(datform.eq.'sca') then c single-column ascii do j=1,n read(7,*,err=232,end=231) x(j) enddo 231 n=j-1 write(6,*) n,' samples read from file ',trim(name) return 232 write(6,*) 'error (maybe a non-numeric character) in data line ' & ,n,' of file ',trim(name) stop else c multi- and variable-column ascii c note: zeile is padded with blanks, therefore there is always a blank after the last digit nz=lin nval=0 call cpu_time(cputime1) do read(7,'(a)',end=131) zeile c read(zeile,'(120a)') zarr nn=0 do i=1,119 c if(zarr(i).ne.b.and.zarr(i+1).eq.b) nn=nn+1 if(zeile(i:i).ne.b.and.zeile(i+1:i+1).eq.b) nn=nn+1 enddo nz=nz+1 read(zeile,*,err=132) (y(j),j=1,nn) do j=1,nn x(nval+j)=y(j) enddo nval=nval+nn enddo 131 write(6,*) nval,' samples read from file ',trim(name) n=min(n,nval) c call cpu_time(cputime2) c write(*,'(a,f7.3,a)') ' in ',cputime2-cputime1,' seconds' close(7) return 132 write(6,*) 'error (maybe a non-numeric character) in data line ' & ,nz,' of file ',trim(name) stop endif c read data in SEIFE format else if(datform.eq.'fil') then read(7,'(a)') text write(6,*) ' header line: ',trim(text) npro=0 21 read(7,'(a)') pro if(pro(1:1).ne.'%'.and.pro(1:1).ne.'}') goto 20 npro=npro+1 if(npro.gt.48) stop ' cannot handle more than 48 header lines' prolog(npro)=pro(3:62) goto 21 20 read(pro,1) nn,iform,dt,tmin,tsec 1 format(i10,a20,3f10.3) write(6,22) trim(name),nn,ndim,dt 22 format(' reading file ',a,' ',i9,'/',i9,' 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'.or.index(iform,'i').gt.0.or. & index(iform,'*').gt.0) then write(6,*) ' data will be read in free-format' read(7,*,err=25,end=23) (x(j),j=1,n) else read(7,iform,end=23) (x(j),j=1,n) endif c read data in MARS88 format else if(datform.eq.'mar') 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) trim(name),n,ndim,dt do j=1,n read(7,*,err=23,end=23) x(j) enddo else if(datform.eq.'asl') then c read data in ASL format (such as written by Quanterra's Cimarron) read(7,'(a)') zeile text=zeile(1:72) write(6,*) 'header: ',trim(text) read(zeile,'(120a)') zarr nn=0 do i=1,119 if(zarr(i).ne.b.and.zarr(i+1).eq.b) then nn=nn+1 ix(nn)=i endif enddo read(zeile(1:ix(1)),'(a)') code1 read(zeile(ix(1)+2:ix(2)),'(a)') code2 read(zeile(ix(2)+2:ix(7)),*) iyear,iday,ithr,itmin,itsec read(zeile(ix(7)+2:ix(11)),*) t1000,srate8,n,cpv srate=srate8 thr=ithr tmin=itmin tsec=itsec dt=nint(1.d6/srate8)/1.d6 tmin=tmin+60.*ithr tsec=tsec+t1000/1000. nz=1 nval=0 do read(7,'(a)',end=31) zeile c read(zeile,'(120a)') zarr nn=0 call cpu_time(cputime1) do i=1,119 c if(zarr(i).ne.b.and.zarr(i+1).eq.b) nn=nn+1 if(zeile(i:i).ne.b.and.zeile(i+1:i+1).eq.b) nn=nn+1 enddo read(zeile,*,err=32) (ix(j),j=1,nn) do j=1,nn x(nval+j)=ix(j) enddo nval=nval+nn enddo 31 write(6,*) nval,' samples read from file ',trim(name) n=min(n,nval) c call cpu_time(cputime2) c write(*,'(a,f7.3,a)') ' in ',cputime2-cputime1,' seconds' close(7) return 32 write(6,*) 'error (maybe a non-numeric character) in data line ' & ,nz,' of file ',trim(name) stop else if(datform.eq.'psn') then c read PSN single-column ascii with multi-line header read(7,'(a)') zeil if(zeil(1:11).ne.'! PSN ASCII') stop 'not a PSN ASCII file' do read(7,'(a)') zeil if(zeil(1:11).eq.'Start Time:') then read(zeil(24:25),*) thr read(zeil(27:28),*) tmin read(zeil(30:35),*) tsec tmin=tmin+60.*thr else if(zeil(1:4).eq.'SPS:') then read(zeil(6:9),*) srate8 srate=srate8 dt=nint(1.d6/srate8)/1.d6 else if(zeil(1:18).eq.'Number of Samples:') then read(zeil(20:30),*) n else if(zeil(1:5).eq.'Data:') then exit endif enddo do j=1,n read(7,*,err=233,end=234) x(j) enddo write(6,*) n,' samples read from file ',trim(name) return 234 n=j-1 write(6,*) n,' samples read from file ',trim(name) return 233 write(6,*) 'error (maybe a non-numeric character) in data line ' & ,n,' of file ',trim(name) stop else if(datform.eq.'bdf') then c read data in BDF format first=.true. nnn=0 40 read(7,'(a)',end=23) text code=text(1:4) rest=text(5:44) c write(*,*) code,rest if(code.eq.'DATA') then goto 41 else if(code.eq.'RATE') then read(rest,*) srate8 srate=srate8 dt=nint(1.d6/srate8)/1.d6 else if(first.and.code.eq.'TIME') then read(rest(1:4),*) year read(rest(6:8),*) day read(rest(10:11),*) thr read(rest(13:14),*) tmin read(rest(16:22),*) tsec tmin=tmin+60.*thr first=.false. else if(code.eq.'NSAM') then read(rest,*) nsam endif goto 40 41 read(7,*,err=26,end=23) (x(j),j=nnn+1,nnn+nsam) nnn=nnn+nsam goto 40 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) trim(name),j 26 format(' Input error reading file ',a,' at sample # ',i7) stop end subroutine outform(name,par,text,prol,npro,ncl,n,dt,tmin,tsec,x, & cfgline,lines) c formatted output dimension x(n) character name*48,text*72,par*56,iform*20,fform*22 character*6 code1*6,code2*6,prol(48)*60,cfgline(32)*256 common /asl/ code1,code2,year,day,thr,tm,ts,t100,srate,cpv logical ncl if(npro.gt.48) stop ' cannot handle more than 48 header lines' iform=' ' xmax=0. do 2 j=1,n 2 xmax=max(xmax,abs(x(j))) if(index(par,'lab').gt.0) then nvor=log10(max(xmax/cpv,2.))+1. ndec=max(0,10-nvor) else nvor=log10(max(xmax,2.))+1. ndec=max(0,10-nvor) endif i1=index(par,'(') i2=index(par,')') if(trim(par).eq.'') par(1:3)='fre' if(index(par,'fre').gt.0.) then write(iform,3) ndec 3 format('(5f13.',i1,')') else if(index(par,'asc').gt.0.) then write(iform,7) ndec 7 format('(f13.',i1,')') else if(index(par,'lab').gt.0) then write(iform,4) ndec 4 format('(1x,f10.3,1x,f13.',i1,') ') else if(index(par,'nil').gt.0) then return else if(i1.gt.0.and.i2.gt.i1+4) then read(par(i1:i2),'(a)') iform else write(*,*) ' illegal format ',par stop endif write(6,5) trim(name),n,iform 5 format(' writing file ',a/i9,' samples in format ',a) if(index(par,'lab').gt.0) then open(8,file=name,status='new') write(8,iform) ((j-1)*dt,x(j)/cpv,j=1,n) call cfgfile(name,n,cfgline,lines,xmax/cpv,dt) else if(index(par,'asc').gt.0.) then open(8,file=name) write(8,iform) (x(j),j=1,n) else ndf=min(8,8-floor(alog10(dt))) write(fform,'(a13,i1,a8)') '(i10,a20,f10.',ndf,',2f10.3)' open(8,file=name) write(8,'(a)') text if(.not.ncl) write(8,'("% ",a)') (prol(j),j=1,npro) write(8,fform) n,iform,dt,tmin,tsec 1 format(i10,a20,f10.5,2f10.3) write(8,iform) (x(j),j=1,n) endif close(8) write(6,'(" finished!")') return 99 stop ' insufficient or wrong parameters for end command' end subroutine plotwindow(nil,par,npl1,npl2,n,ntrend) dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) npl1,npl2 npl1=max(npl1,1) npl2=min(npl2,n) if(npl1.gt.npl2) then stop 'plw: impossible plot window' endif write(6,*) 'plotwindow set: ',npl1,npl2 nil=.false. ntrend=1 return 99 stop 'plw parameters missing or unreadable' end subroutine cfgfile(name,n,cfgline,lines,xmax,dt) c Erzeugt eine Zeile des cfg-Files character name*48,code1*6,code2*6 character*256 cfgline(32) common /asl/ code1,code2,year,day,thr,tm,ts,t100,srate,cpv lines=lines+1 g=1. ifilt=0 tlp=2.5*dt thp=0. write(cfgline(lines),10) trim(name),2.*xmax,g,srate,ifilt, & tlp,thp,1./dt,n 10 format(a,'; Range:',e10.3,'; Gain:',f6.1,'; A/DFreq:',f8.3, & '; A/DFiltCode:',i2,'; BWTPCornerPer:',f8.3,'; BWHPCornerPer:', & f8.3,'; PrintFreq:',f8.3,'; Samples:',i7) return end 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,nam2,fdt,lin,nseis,zeile,npts,nch,ncl,npl, & npw,pli,frq,cfgname,inpath,outpath,datform) character nam1(32)*48,nam2(32)*48,zeile*99,z3*3 character*48 inpath,outpath,cfgname character*3 datform(32) dimension fdt(32),lin(32) logical ncl,npl,npw,pli,nopli,frq data nopli /.false./ nseis=0 101 read(7,'(a)',end=2) zeile if(zeile(1:1).ne.' ') write(6,'(1x,a)') trim(zeile) z3=zeile(1:3) if(z3.eq.'fil'.or.z3.eq.'asc'.or.z3.eq.'asl'.or.z3.eq.'mar' & .or.z3.eq.'bdf'.or.z3.eq.'sca'.or.z3.eq.'psn') then nseis=nseis+1 datform(nseis)=z3 if(z3.ne.'fil') nopli=.true. if(z3.eq.'mar'.or.z3.eq.'asl'.or.z3.eq.'bdf'.or.z3.eq.'psn') then lin(nseis)=1 endif do 30 i=5,99 30 if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 31 31 ii=i do 11 i=ii,99 11 if (zeile(i:i).eq.' '.or.zeile(i:i).eq.',') goto 34 34 if(i.gt.ii+48) goto 24 nam1(nseis)=zeile(ii:i-1) ii=i do 12 i=ii,99 12 if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 32 32 ii=i do 13 i=ii,99 13 if (zeile(i:i).eq.' '.or.zeile(i:i).eq.',') goto 33 33 if(i.gt.ii+48) goto 24 nam2(nseis)=zeile(ii:i-1) if(z3.eq.'asc'.or.z3.eq.'sca') then ii=i do 14 i=ii,99 14 if (zeile(i:i).ne.' '.and.zeile(i:i).ne.',') goto 4 4 if(i.gt.97) goto 26 c if(z3.eq.'asc') then read(zeile(i:99),*,err=26) fdt(nseis),lin(nseis) c else c read(zeile(i:99),*,err=26) lin(nseis) c lin(nseis)=-lin(nseis) c endif endif else if(z3.eq.'lim') then read(zeile,'(4x,i10)') npts else if(z3.eq.'chp') then read(zeile,'(4x,i10)') nch else if(z3.eq.'ncl') then ncl=.true. else if(z3.eq.'npl') then npl=.true. else if(z3.eq.'npw') then npw=.true. else if(z3.eq.'pli') then pli=.true. else if(z3.eq.'pfi') then inpath=zeile(5:28) else if(z3.eq.'pfo') then outpath=zeile(5:28) else if(z3.eq.'cfi') then cfgname=zeile(5:28) else if(z3.eq.'frq') then frq=.true. else if(z3.eq.'end') then if(pli.and.nopli) write(6,*) ' pli command will be ignored' pli=pli.and..not.nopli if(nseis.gt.0) then c Pfade vor die Filenamen setzen if(inpath(1:1).ne.'#') then do i=1,nseis nam1(i)=trim(inpath)//'\'//nam1(i) enddo endif if(outpath(1:1).ne.'#') then do i=1,nseis nam2(i)=trim(outpath)//'\'//nam2(i) enddo endif return else stop ' Filename missing' endif endif goto 101 2 stop 'missing end line or missing after end line' 24 write(*,*) ' Filename too long: ',zeile stop 26 stop ' asc or sca statement incorrect or incomplete' 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 ",e13.6)') f nil=.false. return 99 stop ' parameter for fac missing or unreadable' 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,a3,a9,e13.6)') typ,' Faktor ',f nil=.false. return 99 write(*,*) ' illegal parameter for ',typ stop end subroutine pad(nil,par,x,n,nd) c Nullen ergaenzen dimension x(nd) character*56 par logical nil read(par,*,err=99,end=99) nul if(nul.le.n.or.nul.gt.nd) goto 99 do 1 j=n+1,nul 1 x(j)=0. n=nul write(6,2) nul 2 format(" pad: signal was padded with zeros to length",i8) nil=.false. return 99 stop ' parameter for pad missing, illegal, or unreadable' end subroutine delay(nil,par,x,n) c vorne Nullen einfuegen dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) nul if(nul.le.0) goto 99 do 1 j=n,nul+1,-1 1 x(j)=x(j-nul) do 3 j=nul,1,-1 3 x(j)=0. write(6,2) nul 2 format(" del: the signal was delayed by",i8," samples") nil=.false. return 99 stop ' parameter for del missing, wrong, or unreadable' end subroutine plotpar(nil,par,breit,hoch) c Plotgroesse aendern character*56 par logical nil read(par,*,err=99,end=99) breit,hoch write(6,2) breit,hoch 2 format(" plt: plotsize ",f5.1," * ",f5.1) nil=.false. return 99 stop ' parameter for plt missing, wrong, or unreadable' end subroutine add(nil,par,x,n) c add a constant to each sample 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,'(" add ",e13.6)') f nil=.false. return 99 stop ' insufficient or wrong parameters for add' end subroutine distort(nil,par,x,n) c add nonlinear distortion dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) k,a xmax=0. do 2 j=1,n 2 xmax=max(xmax,abs(x(j))) do 1 j=1,n 1 x(j)=x(j)+a*xmax*(x(j)/xmax)**k write(6,'(" distortion added",f12.6," * x**",i1)') a,k nil=.false. return 99 stop ' insufficient or wrong parameters for dis' 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 n1=max(n1,1) n2=min(n2,n) if(n1.gt.n2) goto 99 do 1 j=n1,n2 1 x(j)=0. write(6,'(" null ",2i10)') n1,n2 nil=.false. return 99 stop ' insufficient or wrong parameters for null' end subroutine set(nil,par,x,n) c set time series from n1 to n2 to a dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) n1,n2,a do 1 j=n1,n2 1 x(j)=a write(6,'(" set ",2i10)') n1,n2 nil=.false. return 99 stop ' insufficient or wrong parameters for set' end subroutine clip(nil,par,x,n) c clip signal at level clp dimension x(n) character*56 par logical nil read(par,*,err=99,end=99) clp do 1 j=1,n x(j)=min(x(j),clp) 1 x(j)=max(x(j),-clp) write(6,'(" clip ",e13.6)') clp nil=.false. return 99 stop ' insufficient or wrong parameters for clip' end subroutine deriv(nil,par,x,n,dt) c derivative (non-recursive, non-causal, symmetric-difference) dimension x(n) character*56 par logical nil twodt=2.*dt 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 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 - in",i4," places, ",i5," values >",e13.6, & " were interpolated. 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 stop ' *** spl *** insufficient number of pivots' 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 stop ' insufficient or wrong parameters for spl' end subroutine window(nil,pli,typ,par,x,n,dt,tmin,tsec) c time window by sample number dimension x(n) character*56 par character*3 typ logical nil,pli if(pli.and.typ.ne.'tap') write(*,*) & " WARNING when using PLI, don't apply any time window" if(typ.eq.'skp'.or.typ.eq.'cos'.or.typ.eq.'tap') 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 c *** tap if(typ.ne.'tap') goto 16 pi=4.*atan(1.) if(np1.gt.n/2) write(6,*) 'np1 set to n/2' np1=min(np1,n/2) fak=pi/np1 do 14 j=1,np1 tap=0.5-0.5*cos(fak*(j-1)) x(j)=x(j)*tap 14 x(n-j+1)=x(n-j+1)*tap acor=1./sqrt(1.-np1*1.25/n) write(6,19) acor 19 format(' amplitude correction factor ',f6.3,' applied') do 18 j=1,n 18 x(j)=x(j)*acor goto 15 c *** cos 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 c *** win, swi 13 if(typ.eq.'win'.or.typ.eq.'swi') then do 11 j=1,n 11 x(j)=x(np1+j) endif c *** sin, sis if(typ.eq.'sin'.or.typ.eq.'sis') then pi=4.*atan(1.) f=pi/(n+1) w2=sqrt(2.) nex=1 if(typ.eq.'sis') nex=2 do 12 j=1,n 12 x(j)=w2*x(np1+j)*sin(f*j)**nex 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(*,*) ' 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 time 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 ' impossible time window' 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 stop ' insufficient or wrong parameters for twi' end subroutine timeset(nil,par,tmin,tsec) character*56 par logical nil read(par,*,err=99,end=99) tmin,tsec write(6,*) ' time stamp set to tmin=',tmin,', tsec=',tsec nil=.false. return 99 stop ' insufficient or wrong parameters for tim' end subroutine decim(nil,par,x,n,dt) c decimate time series by a factor up to 100 dimension x(n),g(200) character*56 par logical nil read(par,*,err=99,end=99) idez if(idez.gt.200) stop ' cannot decimate by more than 200' 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,'(" time series was decimated by ",i10)') idez nil=.false. return 99 stop ' insufficient or wrong parameters for dec' end subroutine mean(nil,par,x,n) c remove average character*56 par dimension x(n) logical nil if(trim(par).eq.'') write(par,*) 0 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 ",e13.6)') sum nil=.false. return 99 stop ' insufficient or wrong parameters for avg' end subroutine quant(nil,par,x,n) c quantization character*56 par dimension x(n) logical nil if(trim(par).eq.'') write(par,*) 1. read(par,*,err=99,end=99) a if(a.le.0.) goto 99 do 1 j=1,n 1 x(j)=a*anint(x(j)/a) write(6,'(" qnt: signal was quantized to integer *",f10.6)') a nil=.false. return 99 stop ' quantization interval must be positive' 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 if(trim(par).eq.'') write(par,*) 0 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 ",2f13.6)') a,b nil=.false. return 99 write(*,*) ' insufficient or wrong parameters for ',typ stop end subroutine shift(nil,typ,par,x,n) c shift signal left or right implicit real*8 (a-h,o-z) real*4 x(n) logical nil character typ*3,par*56 read(par,*,err=99,end=99) ns if(ns.eq.0) then nil=.false. return else if(ns.lt.0) then goto 99 endif if(typ.eq.'shl') then do i=1,n-ns x(i)=x(i+ns) enddo do i=n-ns+1,n x(i)=0. enddo else if(typ.eq.'shr') then do i=n,ns+1,-1 x(i)=x(i-ns) enddo do i=ns,1,-1 x(i)=0. enddo endif write(6,'(1x,a3,i7,"signal shifted in time")') typ,ns nil=.false. return 99 write(*,*) ' insufficient or wrong parameters for ',typ stop end subroutine zspline(nil,typ,par,x,n) c remove z-trend with a z-spline implicit real*8 (a-h,o-z) real*4 x(n),fn32,x2,x3 logical nil character typ*3,par*56 common /tre/ a,b read(par,*,err=99,end=99) n1,n2,n3,n4 fn32=n3-n2 x2=x(n2) x3=x(n3) if(n1.lt.1.or.n2.lt.n1+1.or.n3.lt.n2+1.or.n4.lt.n3+1.or.n4.gt.n) & then stop ' illegal parameters for zspline' endif call trend(nil,typ,par,x(n1),n2-n1+1) a1=a b1=b call trend(nil,typ,par,x(n3),n4-n3+1) a3=a b3=b do 2 j=1,n2 2 x(j) = x(j)-sngl(a1+b1*j) do 3 j=n2+1,n3-1 3 x(j) = x(j)-(x2+(j-n2)/fn32*(x3-x2)) do 4 j=n3,n4 4 x(j) = x(j)-sngl(a3+b3*(j-n3+1)) write(6,'(" zspline eck ",4i10,)') n1,n2,n3,n4 write(6,'(" zspline coef",4f10.6)') a1,b1,a3,b3 nil=.false. return 99 stop ' insufficient or wrong parameters for zspline' end subroutine polytrend(nil,par,x,n) c remove polynomial trend parameter(ndi=13) 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(1x,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 stop ' insufficient or wrong parameters for pol' 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 stop ' insufficient or wrong parameters for exp' end subroutine sig(nil,zeile,x,n) c signum function or absolute value real*4 x(n) logical nil character zeile*99 if(zeile(1:3).eq.'sig') then do i=1,n x(i)=sign(1.,x(i)) enddo write(6,'(a)') ' sig: signum function applied' else do i=1,n x(i)=abs(x(i)) enddo write(6,'(a)') ' abs: absolute value taken' endif nil=.false. return 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 subroutine tides(nil,par,x,n,dt) c remove tides. number of frequencies is automatically chosen according c to the total length of the record. implicit real*8 (a-h,o-z) real*4 x(n) logical nil character*56 par dimension a(13,13),rs(13),c(13),d(13),e(13),f(13), & omega(6),omeg(6),cor(6) data omega/1.93227d0,0.92954d0,2.00000d0,1.00274d0,1.89567d0, & 0.89293d0/ data zero,one,two/0.d0,1.d0,2.d0/ ndi=13 dth=dt/two read(par,*,err=99,end=99) nstep if(nstep.eq.0) nstep=300./dts nstep=max(nstep,1) step=nstep tstep2=(step-one)*dth tint=tstep2+dth c determine the number of frequencies required for a good fit dur=n*dt/3600.d0 nfreq=6 if(dur.lt.35.d0) nfreq=5 if(dur.lt.18.d0) nfreq=4 if(dur.lt.14.d0) nfreq=3 if(dur.lt. 5.d0) nfreq=2 nco=2*nfreq+1 omeg0=8.d0*datan(one)/86400.d0 do 101 i=1,nfreq 101 omeg(i)=omeg0*omega(i) c determine partial amplitudes by least-squares fit do 1 i=1,nco rs(i)=zero do 1 k=1,nco 1 a(i,k)=zero c correction for averaging over nstep samples do 102 j=1,nfreq 102 cor(j)=step*dsin(omeg(j)*dth)/dsin(step*omeg(j)*dth) c set up system of linear equations c(1)=one do 2 j=1,n-nstep+1,nstep sx=zero do 12 jj=j,j+nstep-1 12 sx=sx+x(jj) sx=sx/step t=(j-1)*dt+tstep2 do 103 k=1,nfreq c(2*k) =dcos(omeg(k)*t) 103 c(2*k+1)=dsin(omeg(k)*t) do 2 i=1,nco rs(i)=rs(i)+sx*c(i) do 2 k=1,nco 2 a(i,k)=a(i,k)+c(i)*c(k) c solve for partial amplitudes call gauss(a,nco,ndi,rs,f) do 3 j=1,n-nstep+1,nstep t=(j-1)*dt+tstep2 c remove average and tides do 104 k=1,nfreq k2=2*k k3=k2+1 c(k2)=dcos(omeg(k)*t) c(k3)=dsin(omeg(k)*t) d(k2)=-omeg(k)*c(k3) d(k3)= omeg(k)*c(k2) e(k2)=-omeg(k)*d(k3) 104 e(k3)= omeg(k)*d(k2) xgez1=f(1) xgez2=zero xgez3=zero do 105 k=1,nfreq k2=2*k k3=k2+1 xgez1=xgez1+cor(k)*(f(k2)*c(k2)+f(k3)*c(k3)) xgez2=xgez2+cor(k)*(f(k2)*d(k2)+f(k3)*d(k3)) 105 xgez3=xgez3+cor(k)*(f(k2)*e(k2)+f(k3)*e(k3))/two if(j.gt.n-2*nstep+1) nstep=n+1-j do 3 jj=j,j+nstep-1 tdif=(jj-j)*dt-tstep2 3 x(jj)=x(jj)-sngl(xgez1+xgez2*tdif+xgez3*tdif*tdif) write(6,9) f(1),(omega(j),f(2*j),f(2*j+1),j=1,nfreq) 9 format(" tid avg=",f10.0/(5x,"freq=",f10.5," cos=", & e11.3," sin=",e11.3)) nil=.false. return 99 stop ' insufficient or wrong parameters for tides' end c subroutine trans(nil,typ,par,x,n,dt) c remove transient response implicit real*8 (a-h,o-z) real*4 x(n),dts logical nil character typ*3, par*56 dimension a(3,3),rs(3),c(3),f(3) data zero,one/0.d0,1.d0/ twopi=8.d0*datan(one) if(typ.eq.'tr3') then read(par,*,err=99,end=99) n1,t0,h if(h.ge.one) then stop ' tra: h must be less than 1' endif c determine partial amplitudes by least-squares fit wh=dsqrt(one-h**2) omeg=twopi/t0 do 1 i=1,3 rs(i)=zero do 1 k=1,3 1 a(i,k)=zero c set up system of linear equations c(1)=one do 2 j=n1,n t=(j-1)*dt ex=dexp(-omeg*h*t) c(2)=dcos(omeg*wh*t)*ex c(3)=dsin(omeg*wh*t)*ex do 2 i=1,3 rs(i)=rs(i)+x(j)*c(i) do 2 k=1,3 2 a(i,k)=a(i,k)+c(i)*c(k) c solve for partial amplitudes call gauss(a,3,3,rs,f) else c use previously determined partial amplitudes read(par,*,err=99,end=99) n1,t0,h,f(1),f(2),f(3) wh=dsqrt(one-h**2) omeg=twopi/t0 endif do 5 j=1,n1-1 5 x(j)=0. do 3 j=n1,n t=(j-1)*dt ex=dexp(-omeg*h*t) omt=omeg*wh*t 3 x(j)=x(j)-sngl(f(1)+f(2)*dcos(omt)*ex+f(3)*dsin(omt)*ex) write(6,9) f 9 format(" transient: avg=",f13.3," f2,f3=",2f13.3) nil=.false. return 99 stop ' insufficient or wrong parameters for trans' end subroutine filter(nil,zeile,x,n,dt,npw,frq) implicit real*8 (a-h,o-z) character typ*3,typ1*3,typ2*3,par*56,zeile*99,warp*13 real*4 x(n) logical nil,npw,frq typ=zeile(1:3) par=zeile(5:60) pi=4.d0*datan(1.d0) 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 if(frq) t0=1./t0 it=11 else if(typ.eq.'hp1') then read(par,*,err=99,end=99) t0 if(frq) t0=1./t0 it=12 else if(typ.eq.'lp2') then read(par,*,err=99,end=99) t0,h if(frq) t0=1./t0 it=21 else if(typ.eq.'hp2') then read(par,*,err=99,end=99) t0,h if(frq) t0=1./t0 it=22 else if(typ.eq.'bp2') then read(par,*,err=99,end=99) t0,h if(frq) t0=1./t0 it=24 else if(typ.eq.'int') then it= 1 else if(typ.eq.'he1'.or.typ.eq.'le1') then read(par,*,err=99,end=99) t0s,t0 if(frq) t0s=1./t0s if(frq) t0=1./t0 it=13 else if(typ.eq.'he2'.or.typ.eq.'le2') then read(par,*,err=99,end=99) t0s,hs,t0,h if(frq) t0s=1./t0s if(frq) t0=1./t0 it=23 else if(typ.eq.'sez') then read(par,*,err=99,end=99) t0s,hs if(frq) t0s=1./t0s t0=1.e12 h=1. it=23 else goto 2 endif c prewarp and calculate filter weights if(it.eq.1) then t0p=1.d0/dt t0sp=t0p else if(npw) then t0p=t0/dt t0sp=t0s/dt warp=' no prewarp ' else t0p=pi/dtan(pi/t0*dt) t0sp=pi/dtan(pi/t0s*dt) warp=' prewarp ' endif endif call rfk(it,t0p,h,t0sp,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) nil=.false. c confirm execution if(it.eq.1) write(6,'(1x,a3)') typ if(it.eq.11.or.it.eq.12) write(6,8) typ,warp,t0 if(it.eq.13) write(6,9) typ,warp,t0s,t0 7 format(1x,a3,a13,2f10.3,f16.3,f10.3) 8 format(1x,a3,a13,5f10.3) 9 format(1x,a3,a13,f10.3,f16.3) if(it.eq.21.or.it.eq.22.or.it.eq.24) write(6,8) typ,warp,t0,h if(it.eq.23) then if(typ.eq.'sez') then write(6,8) typ,warp,t0s,hs else write(6,7) typ,warp,t0s,hs,t0,h endif endif write(*,11) " f0,f1,f2,g1,g2= ",f0,f1,f2,g1,g2 11 format(a,5f12.6) write(*,*) return c Butterworth filters 2 if(typ.eq.'lpb'.or.typ.eq.'hpb') then read(par,*,err=99,end=99) t0,m if(frq) t0=1./t0 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 c prewarp and calculate filter weights if(it.eq.1) then t0p=1.d0/dt t0sp=t0p else if(npw) then t0p=t0/dt t0sp=t0s/dt warp=' no prewarp ' else t0p=pi/dtan(pi/t0*dt) t0sp=pi/dtan(pi/t0s*dt) warp=' prewarp ' endif endif 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,warp,t0 endif pih=2.d0*datan(1.d0) do 3 j=1,mm h=dsin(pih*(2*j-1)/m) c prewarp and calculate filter weights if(it.eq.1) then t0p=1.d0/dt t0sp=t0p else if(npw) then t0p=t0/dt t0sp=t0s/dt warp=' no prewarp ' else t0p=pi/dtan(pi/t0*dt) t0sp=pi/dtan(pi/t0s*dt) warp=' prewarp ' endif endif 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,warp,t0,h 3 continue nil=.false. endif return 99 write(*,*) ' insufficient or wrong parameters for ',typ stop end subroutine extrem(nil,x,npl1,npl2,n) dimension x(n) logical nil nil=.false. xmin=1e96 xmax=-xmin amin=xmin do 1 j=npl1,npl2 if(abs(x(j)).gt.amin) goto 4 amin=abs(x(j)) namin=j 4 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 if(xmax.ge.-xmin) then amax=xmax namax=nmax else amax=-xmin namax=nmin endif write(6,*) 'ext: extrema between samples ',npl1,' and ',npl2,':' write(6,2) xmin,nmin,xmax,nmax 2 format(" xmin=",e11.3," (j=",i7,") xmax=",e11.3," (j=",i7,")") write(6,5) amin,namin,amax,namax 5 format(" amin=",e11.3," (j=",i7,") amax=",e11.3," (j=",i7,")") write(6,6) xmax-xmin 6 format(" peak-to-peak ",e11.3) return end subroutine rms(nil,x,npl1,npl2,n) dimension x(n) logical nil nil=.false. avg=0. do 3 j=npl1,npl2 3 avg=avg+x(j) avg=avg/(npl2-npl1+1) xrms=0. do 1 j=npl1,npl2 1 xrms=xrms+(x(j)-avg)**2 write(6,2) npl1,npl2,sqrt(xrms/(npl2-npl1+1)) 2 format(" rms over de-meaned samples ",i8," to ",i8,": ",e11.3) return end subroutine square(nil,x,n) dimension x(n) logical nil nil=.false. do 1 j=1,n 1 x(j)=x(j)**2 write(6,2) 2 format(" square: the signal has been squared") return end subroutine intpo(nil,par,dt1,x,np1) dimension x(np1) logical nil character*56 par nil=.false. read(par,*,err=99,end=99) dt2 if(dt2.le.dt1) then stop ' sampling inteval can be increased but not decreased' endif j2=0 1 j2=j2+1 t=(j2-1)*dt2 j10=t/dt1 frac=t/dt1-j10 j11=j10+1 j12=j11+1 if(j12.gt.np1) goto 2 x(j2)=(1.-frac)*x(j11)+frac*x(j12) goto 1 2 np1=j2-1 write(6,3) dt1,dt2 3 format(" csi ",2f10.3) dt1=dt2 return 99 stop ' insufficient or wrong parameters for intpo' end c subroutine heapsort(nil,ra,n) dimension ra(n) logical nil l=n/2+1 ir=n 10 continue if(l.gt.1) then l=l-1 rra=ra(l) else rra=ra(ir) ra(ir)=ra(1) ir=ir-1 if(ir.eq.1) then ra(1)=rra goto 30 endif endif i=l j=l+l 20 if(j.le.ir) then if(j.lt.ir) then if(ra(j).lt.ra(j+1)) j=j+1 endif if(rra.lt.ra(j)) then ra(i)=ra(j) i=j j=j+j else j=ir+1 endif goto 20 endif ra(i)=rra goto 10 30 nil=.false. write(*,*) ' series was heap-sorted, is no longer a signal' return end