# Rolling STA/LTA 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. ```r roll_stalta(x, n_sta, n_lta, increment = 1) ``` ## Arguments - `x`: an 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 ## Returns A vector of values of the same length as `x` with each point containing the STA/LTA ratio at that point. ## 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}\sum_{j=i}^{i+ns}{x_i} $$ $$ LTA(x_i) = \frac{1}{nl}\sum_{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 `NA`s for all skipped-over indices. ## Examples ```r # 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) ``` ## References [First Break Picking](https://en.wikipedia.org/wiki/First_break_picking)

roll_stalta function