program stalta dimension cmp1(30000), cmp2(30000) integer npts,nenvfil, count real del, caz1, caz2, sta,lta,cutoff,hanfil,trans character*60 f1, f2 character*16 kstnm character*60 tempdir, rd, tr real*4 picks1(2,1000), picks2(2,1000), picks3(2,1000) integer npks common b,e,cmpi,stla,stlo,evla,evlo,evdp,gcarc,baz, *ny,nd,nh,nm,ns,nms,kstnm,kval,dt c c READ INPUT PARAMETERS c read(5,'(A)')f1 read(5,'(A)')f2 read(5,'(A)')tempdir read(5,'(F6.3)')sta read(5,'(F6.3)')lta read(5,'(F6.3)')cutoff read(5,'(I3)')nenvfil read(5,'(F6.3)')hanfil c call sacin(f1,cmp1,npts,del,caz1) call sacin(f2,cmp2,npts,del,caz2) c print*,'delta= ',dt c print*,'b= ',b c print*,'e= ',e c print*,'cmpi= ',cmpi c print*,'cmpaz= ',cmpaz c print*,'stla= ',stla c print*,'stlo= ',stlo c print*,'evla= ',evla c print*,'evlo= ',evlo c print*,'gcarc= ',gcarc c print*,'baz= ',baz c print*,'year= ',ny c print*,'day= ',nd c print*,'hour= ',nh c print*,'minute= ',nm c print*,'second= ',ns c print*,'msec= ',nms c print*,'npts= ',npts c print*,'stnm= ',kstnm c c sta = 1. c lta = 7. c cutoff = 2. c nenvfil = 50. c hanfil = 0.33 c c print*,'file= ',f1 c print*,'file= ',f2 c print*,'sta= ',sta c print*,'lta= ',lta c print*,'CUTOFF= ',cutoff c print*,'nenvfil= ',nenvfil c print*,'hanfil= ',hanfil c c ROTATE INTO BAZ c call rotbaz(cmp1,cmp2,baz,npts,caz1,caz2) c c OUTPUT RADIAL + TRANSVERSE c rd = tempdir(1:lnblnk(tempdir))//'file.r' tr = tempdir(1:lnblnk(tempdir))//'file.t' trans = baz + 90. call sacout(rd,cmp1,npts,del,baz) call sacout(tr,cmp2,npts,del,trans) call pickseis(cmp1,npts,del,b,sta,lta,cutoff, & nenvfil,hanfil,picks3,npks) do 300 count=1,1000,1 if (picks3(1,count) .ne. 0.0) then print*,picks3(1,count),picks3(2,count) endif 300 continue c print*,picks c stop c end c ********************************************************************* c subroutine sacin(fle,dat,npts,del,cmpaz) c parameter(max=30000) integer ny,nd,nh,nm,ns,nm real b,e,cmpi,stla,stlo,evlo, evdp,evla,gcarc,baz character*16 kstnm,kval character*60 fle common b,e,cmpi,stla,stlo,evla,evlo,evdp,gcarc,baz, *ny,nd,nh,nm,ns,nms,kstnm,kval,dt c dimension dat(1) c call rsac1(fle,dat,npts,beg,del,max,nerr) c if(nerr.eq.0) goto 999 c write(6,*)' ATTENTION MAX DIMENSION IS EXCEDEED,/, * TERMINATION .... ' return c 999 call getfhv('delta',dt,nerr) call getfhv('b',b,nerr) call getfhv('e',e,nerr) call getfhv('cmpinc',cmpi,nerr) call getfhv('cmpaz',cmpaz,nerr) call getfhv('stla',stla,nerr) call getfhv('stlo',stlo,nerr) call getfhv('evla',evla,nerr) call getfhv('evlo',evlo,nerr) call getfhv('evdp',evdp,nerr) call getfhv('gcarc',gcarc,nerr) call getfhv('baz',baz,nerr) call getnhv('nzyear',ny,nerr) call getnhv('nzjday',nd,nerr) call getnhv('nzhour',nh,nerr) call getnhv('nzmin',nm,nerr) call getnhv('nzsec',ns,nerr) call getnhv('nzmsec',nms,nerr) call getnhv('npts',npts,nerr) call getkhv('kstnm',kstnm,nerr) c call getihv('iztype',kval,nerr) return end c c c ************************* SACOU ************************************* c c WRITE DATA IN SAC FORMAT FROM A FORTRAN CODE c c ********************************************************************* c subroutine sacout(fil,dat,npts,del,caz) c integer ny,nd,nh,nm,ns,nm real b,e,cmpi,stla,stlo,evlo,evdp,evla,gcarc,baz,caz, del character*16 kstnm,kval character*60 fil c common b,e,cmpi,stla,stlo,evla,evlo,evdp,gcarc,baz, *ny,nd,nh,nm,ns,nms,kstnm,kval,dt c dimension dat(1) c print*,'cmpaz= ',cmpaz c call newhdr call setfhv('delta',del,nerr) call setfhv('b',b,nerr) call setfhv('e',e,nerr) call setnhv('nzyear',ny,nerr) call setnhv('nzjday',nd,nerr) call setnhv('nzhour',nh,nerr) call setnhv('nzmin',nm,nerr) call setnhv('nzsec',ns,nerr) call setnhv('nzmsec',nms,nerr) call setnhv('npts',npts,nerr) call setfhv('cmpinc',cmpi,nerr) call setfhv('cmpaz',caz,nerr) call setfhv('stla',stla,nerr) call setfhv('stlo',stlo,nerr) call setfhv('evla',evla,nerr) call setfhv('evlo',evlo,nerr) call setfhv('evdp',evdp,nerr) call setfhv('gcarc',gcarc,nerr) call setfhv('baz',baz,nerr) call setkhv('kstnm',kstnm,nerr) c call setihv('iztype',kval,nerr) c call wsac0(fil,xdummy,dat,nerr) c return end c c c ************************* ROTBAZ ************************************ c c ROTATE TWO ORTHOGONAL COMPONENTS TO BACKAZIMUTH c c ********************************************************************* c subroutine rotbaz(x1,x2,baz,npts,az1,az2) c c ... Rotates horizontal components,x1&x2 into radial&transverse components c ... given backazimuth : baz c ... and their orientations : az1 , az2 c ... Radial is positive away from the source c ... Transverse is positive clockwise from + radial direction c dimension x1(1),x2(1) c rad(deg)=deg / 57.295779 azck = az1 + 90.0 diff = abs(azck-az2) if(diff.lt.0.01) goto 110 if(diff.lt.179.)then write(6,*)' PROBLEM IN ROTATION ',az1,az2,azck return else endif 110 aa=sin(rad(baz)-rad(az1)) bb=cos(rad(baz)-rad(az1)) do 4 i = 1,npts radial = - x1(i)*bb - x2(i)*aa transv = x1(i)*aa - x2(i)*bb x1(i)=radial x2(i)=transv 4 continue return end ************************************************************************** * Version 1.0 April 2, 1993 * Version 1.1 Nov 15, 2002 * Paul Earle, pearle@.usgs.gov * * * PICKSEIS generates a list of picks for an input seismogram using an * STA/LTA algorithm. * * INPUT: seis = seismogram * npts = number of points in seismogram * dt = sample interval (s) * tstart = start time of record (s after origin time) * sta = short term average (s) * lta = long term average (s) * cutoff = minimum value to register a pick (threshold value) * nenvfil = half-length of filter for envelope function (samples) * = 0 for simple absolute value function * hanfil = half length of the Hanning filter (s) * * OUTPUT: picks(2,1000) pick array * first column pick times (s) * second column pick confidences (magnitude * of the smoothed ratio function following * the trigger point) * * npks number of picks * * * NOTES: * * If there is a large dc offset in the seismogram it is a good idea * to demean it. * * This routine is by no means optimized for speed, and may give * peculiar results for some combinations of input parameters. * * Picks are not made close to the beginning and end of the seismogram. * The earliest that a pick can be made is:(nenvfil*dt) + sta + lta + hanfil (s) * after tstart. * * This algorithm was designed mainly to detect weak arrivals, not to give * precise picks for strong arrivals. However, tweaking of the input * parameters can result in reasonable timing accuracy. The spread * of the auto-pick residuals is less sensitive to changes in input values * then is the mean. * * The following parameter settings have been found to give * satisfactory results: * sta lta cutoff nenvfil hanfil * Global long-period (~25 s waveforms) 8 30 2 100 10 * Global short-period (~1 s waveforms) .5 7 2 50 3 * Local short-period (~10 Hz waveforms) .1 1.5 2 50 0.33 * * However, all data sets are different and it is recommended that the user * test to find the optimal parameters. The algorithm will run faster for * small values of nenvfil and hanfil but too small a value of hanfil will * cause problems. If picks appear too late, increase hanfil; if they are * too early, decrease hanfil. * * To negate problems arising from round-off errors STAs and LTAs smaller * then smallnum (defined by a parameter statement in STALTA) are considered * to be zero. * * Reference: Earle, P., and P.M. Shearer, Characterization of Global * Seismograms Using An Automatic Picking Algorithm, * BSSA, Vol 84, No 2, p. 366-376, Apr 1994 * * Revised: May 11 93 earle * Revised: Oct 28 94 earle changed STALAT to return correct # of 0 values * Revised v1.1: Nov 15 earle * fixed loop range in makepicks * npts-ii+1 changed to npts-2 * * *************************************************************************** subroutine pickseis(seis,npts,dt,tstart,sta,lta,cutoff, & nenvfil,hanfil,picks,npks) parameter (maxpts=30000) !maximum points in seismogram real*4 seis(*),dt,tstart,sta,lta,cutoff,hanfil integer npts,nenvfil real*4 picks(2,1000) integer npks real*4 env(maxpts),ratios(maxpts),smratios(maxpts) C Convert from seconds to samples ista = int(sta/dt) ilta = int(lta/dt) nhanfil = int(hanfil/dt) C Check to see if there are enough points to process seismogram minpts = 2*(nenvfil + ista + ilta + nhanfil) nonzero = 0 do 20 ii = 1,npts if(seis(ii) .NE. 0.) nonzero = nonzero + 1 20 continue if(nonzero .LT. minpts) then print*,'Not non-zero enough points in seismogram to process' npks = 0 goto 999 endif C calculate hilbert transform and envelope function call TILBERT(seis,dt,npts,nenvfil,env) c open(2,file='output.txt') c c do 300 count=1,npts,1 c c if (env(count) .ne. 0.0) then c write(2,'(F10.3,2X,F12.3)'),count,env(count) c endif c c 300 continue C Calculate (long term average)/(short term average) call STALTA(env,npts,ilta,ista,ratios) C Smooth the time series by convolving with a cos**2 filter call HANFILT(ratios,npts,nhanfil,smratios) C Pick phase arrivals call MAKEPICKS(smratios,tstart,dt,npts,cutoff,picks,npks) 999 return end subroutine TILBERT(a,dt,npts,nfil,e) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c c TILBERT computes Hilbert transform and envelope time function c of time series in the time domain. c Inputs: a = time series c dt = time spacing c npts = number of points in a c nfil = half-number of points in filter c = 0 to return simple absolute value function c Returns: b = Hilbert transform of time series NOT RETURNED c e = envelope time function c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC real a(*),b(30000),e(*),h(-1000:1000) save dt0,npts0,h pi=3.141592654 c set e to absolute value if nfil=0 if (nfil.eq.0) then do 5 i=1,npts 5 e(i)=abs(a(i)) return end if C Make filter if new parameters if (dt.ne.dt0.and.npts.ne.npts0) then do 10 i=-nfil,nfil t=float(i)*dt if (i.ne.0) then h(i)=-1./(pi*t) else h(i)=0. end if 10 continue nfiltot=2*nfil+1 call TAPER1(h(-nfil),nfiltot,1.,h(-nfil)) dt0=dt npts0=npts end if C filter series 25 do 30 i=1,npts b(i)=0. e(i)=0. 30 continue i1=1+nfil i2=npts-nfil do 50 i=i1,i2 do 40 j=-nfil,nfil ii=i-j b(i)=b(i)+a(ii)*h(j) 40 continue b(i)=b(i)*dt 50 continue C Calculate envelope function do 60 i=i1,i2 60 e(i)=sqrt(a(i)**2+b(i)**2) return end subroutine TAPER1(x,n,frac,xtap) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c c TAPER1 multiplies a time series by a cosine squared taper c c Inputs: x = input time series c n = number of time points c frac = fraction cosine squared taper c = 1. for Hanning taper c = 0. for no taper c Returns: xtap = tapered time series c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC real x(*),xtap(*) pihalf=3.1415927/2. do 10 i=1,n 10 xtap(i)=x(i) if (frac.eq.0.) go to 999 xlen=float(n-1) taplen=frac*xlen/2. i2=nint(taplen) do 20 i=1,i2 angle=pihalf*float(i-1)/(i2-1) fact=sin(angle)**2 xtap(i)=x(i)*fact xtap(n+1-i)=x(n+1-i)*fact 20 continue 999 return end subroutine STALTA(aa,npts,n1,n2,b) ****************************************************** * * stalta4 finds the ratio of the short term average (STA) * to the long term average (LTA) of a time series and * stores the result in a vector * There is no overlap in the STA and LTA * * INPUT * aa time series * npts number of points in the time series * n1 length of long term window (samples) * n2 length of short term window (samples) * OUTPUT * b vector containing STA/LTA data * the first non zero entry of this * vector will be b(n1+n2) * ****************************************************** parameter (smallnum=1.e-4) !threshold for valid lta and sta real*4 aa(*),b(*) real*4 lta,sta,lta_tot,sta_tot do i= 1,n1+n2 b(i) = 0. enddo C find first non zero do 20 ii = 1,npts if(aa(ii) .NE. 0.) then nz = ii goto 25 endif 20 continue 25 i1 = nz+n1-1 C take the first long term average lta_tot = 0. do j= nz,i1 lta_tot = lta_tot + abs(aa(j)) enddo lta = lta_tot/float(n1) C take the first short term average sta_tot = 0. do j = i1+1,i1+n2 sta_tot = sta_tot + abs(aa(j)) enddo sta= sta_tot/float(n2) if(lta .GE. smallnum .AND. sta .GE. smallnum) then b(i1+n2) = sta/lta else b(i1+n2) = 0. endif C loop over time series do i= i1+1,npts-n2 lta_tot = lta_tot + abs(aa(i)) - abs(aa(i-n1)) sta_tot = sta_tot + abs(aa(i+n2)) - abs(aa(i)) lta = lta_tot/float(n1) sta= sta_tot/float(n2) if(lta .GE. smallnum .AND. sta .GE. smallnum) then b(i+n2) = sta/lta else b(i+n2) = 0. endif enddo return end subroutine HANFILT(a,n,nfil,b) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c c HANFILT convolves a time series with a cos**2 filter c c Inputs: a = input time series c n = number of points in a c nfil = filter half-length (number of points, max=200) c Returns: b = filtered version of a c c NOTE: b is set to zero within nfil of ends c cCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC real a(n),b(n),filt(-800:800) save nfilold,filt,tfilt C Construct filter if (nfil.ne.nfilold) then pihalf=3.1415927/2. tfilt=0. do 10 k=-nfil,nfil s1=float(abs(k))/float(nfil+1) filt(k)=cos(s1*pihalf)**2 tfilt=tfilt+filt(k) 10 continue nfilold=nfil end if C Find first nonzero entry in matrix do 30 ii = 1,n if(a(ii) .NE. 0.) then goto 35 endif 30 continue C Filter data 35 do 50 i=ii+nfil,n-nfil b(i)=0. do 40 k=-nfil,nfil b(i)=b(i)+a(i+k)*filt(k) 40 continue b(i)=b(i)/tfilt 50 continue do 60 i=1,nfil+ii b(i)=0. 60 continue do 70 i=1,nfil b(n+1-i)=0. 70 continue return end subroutine MAKEPICKS(a,tstart,dt,npts,cutoff,p,npk) ******************************************************************* * * MAKEPICKS makes estimates of phase arrival times * * Inputs: a = STA/LTA vector * tstart = time of first element in a(*) (sec) * dt = Time increment (sec) * npts = number of points in a * cutoff = minimum amplitude of local max of a(*) to * be considered a pick. * * Returns: p = array holding pick information * npk = number of picks * ******************************************************************** * real*4 a(*),p(2,*),cutoff,tstart,dt integer npts,npk,idx_min integer idx_inflex real*4 curr_slope,prev_slope real*4 curr_inflex,prev_inflex npk = 0 idx_min = 1 idx_inflex = 0 C Find first nonzero entry in a do 50 ii = 1,npts if (a(ii) .NE. 0.) goto 75 50 continue 75 prev_slope = a(ii+1) - a(ii) prev_inflex = (a(ii+2)+a(ii))/2 - a(ii+1) do 100 i = ii+1,npts-ii+1 curr_slope = a(i+1) - a(i) curr_inflex = (a(i+2)+a(i))/2 - a(i+1) C look for inflection point concave up to concave down if(prev_inflex .GE. 0. .AND. curr_inflex .LE. 0.) & idx_inflex = i+1 C look for local max if(prev_slope .GE. 0. .AND. curr_slope .LT. 0. .AND. & a(i+1) .GE. cutoff .AND. idx_inflex .NE. 0) then npk = npk + 1 p(1,npk) = tstart + dt*idx_inflex p(2,npk) = a(i) endif prev_slope = curr_slope prev_inflex = curr_inflex 100 continue return end