IRISSeismic (version 1.0.7-1)

psdList: Apply McNamara PSD algorithm to a seismic signal

Description

The psdList function breaks up a seismic Stream object into a series of shorter segments with 50% overlap and uses the McNamaraPSD method to return a smoothed (aka binned) Power Spectral Density (PSD) for each segment.

Usage

psdList(st)

Arguments

st

a Stream object

Value

A list of PSD objects is returned. Each element of the list is an R list object with the following elements:

freq, spec, snclq, starttime, endtime

Note: Individual PSDs have not had instrument correction applied.

Details

A Stream will be broken up into segments depending upon the snclq identifier associated with this seismic data. The binning frequencies are also channel dependent as exemplified in this code extract where Z is the segment length in seconds:

alignFreq <- 0.1
  
if (stringr::str_detect(channel,"^L")) {
  Z <- 3 * 3600
  loFreq <- 0.001
  hiFreq <- 0.5 * tr_merged@stats@sampling_rate
} else if (stringr::str_detect(channel,"M")) {
  Z <- 2 * 3600
  loFreq <- 0.0025
  hiFreq <- 0.5 * tr_merged@stats@sampling_rate
} else {
  Z <- 3600
  loFreq <- 0.005
  hiFreq <- 0.5 * tr_merged@stats@sampling_rate
}

Each new segment starts half way through the previous segment. (50% overlap)

References

Seismic Noise Analysis System Using Power Spectral Density Probability Density Functions (McNamara and Boaz 2005)

See Also

McNamaraPSD, psdList2NoiseMatrix, psdPlot, psdStatistics,

Examples

  ## Not run: 
# # Create a new IrisClient
# iris <- new("IrisClient", debug=TRUE)
# 
# # Get seismic data
# starttime <- as.POSIXct("2011-05-05", tz="GMT") # 2011.125
# endtime <- starttime + 1*24*3600
# st <- getDataselect(iris,"IU","GRFO","--","BHE",starttime,endtime)
# 
# # Generate power spectral density for each hour long segment
# psdList <- psdList(st)
# 
# # Plot uncorrected PSDs
# period <- 1/psdList[[1]]$freq
# plot(period, psdList[[1]]$spec, log='x', type='l',
#      xlab="Period (Sec)", ylab="Power (dB)",
#      main="Uncorrected PSDs")
#      
# for (i in seq(2:length(psdList))) {
#   points(period, psdList[[i]]$spec, type='l')
# }
#   ## End(Not run)