# roll_stalta: Rolling STA/LTA In seismicRoll: Fast Rolling Functions for Seismology using Rcpp

## Description

Fast rolling STA/LTA using C++/Rcpp. Additional performance gains can be achieved by skipping increment values between calculations.

The STA/LTA ratio method is used for automatic detection of seismic signal arrival times.

## Usage

 1 roll_stalta(x, n_sta, n_lta, increment = 1) 

## Arguments

 x an R numeric vector n_sta integer STA window size n_lta integer LTA window size increment integer shift to use when sliding the window to the next location

## Details

The roll_stalta function described here does no preprocessing of the incoming data and merely calculates the ratio of the average value in the STA window to the average value in the LTA window. Windows are aligned so that the index is at the left edge of the STA window and at the right edge of the LTA window.

STA(x_i) = \frac{1}{ns}∑_{j=i}^{i+ns}{x_i}

LTA(x_i) = \frac{1}{nl}∑_{j=i-nl}^{i}{x_i}

r_i = \frac{STA_i}{LTA_i}

[---------- LTA --------*]........

.......................[*- STA --]

For proper use of this algorithm seismic data should be preprocessed as in the example below with:

• demean, detrend and taper the raw signal

• square the processed signal to get power

With increment=1, this function is equivalent to, eg:

sta <- roll_mean(x,3,1,"left")

lta <- roll_mean(x,30,1,"right")

r <- sta/lta

For increments greater than one, the rolling means above will not align properly, hence the need for a dedicated roll_stalta function.

Values within n_lta-1 of the beginning or n_sta-1 of the end of x are set to NA.

Setting increment to a value greater than one will result in NAs for all skipped-over indices.

## Value

A vector of values of the same length as x with each point containing the STA/LTA ratio at that point.

## Examples

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 # Contrived example x <- rep(c(1,5,3,2,1),each=20) p <- roll_stalta(x,3,6) plot(x, pch=17, cex=0.8, ylim=c(0,max(x)), main="Test of roll_stalta on artificial data") points(p,cex=1.5,col='red',type='b') legend('topleft', legend=c('data','STA/LTA'), pch=c(17,1), col=c('black','red')) # Real example requiring the 'seismic' package ## Not run: require(seismic) # Create a new IrisClient iris <- new("IrisClient") # Seismic data with a large quake starttime <- as.POSIXct("2010-02-27 06:30:00", tz="GMT") endtime <- as.POSIXct("2010-02-27 07:00:00", tz="GMT") st <- getDataselect(iris,"IU","ANMO","00","BHZ",starttime,endtime) tr <- st@traces[[1]] # Preprocess the data x <- DDT(tr)@data # Calculate the first break 'picker' n_sta <- 3 * tr@stats@sampling_rate n_lta <- 10 * n_sta p <- roll_stalta(x^2,n_sta,n_lta) first_break <- which(p == max(p,na.rm=TRUE)) plot(x,type='l', main='Test of STA/LTA first break picker on raw seismic data') abline(v=first_break,col='red') ## End(Not run) 

seismicRoll documentation built on May 29, 2017, 5:40 p.m.