IRISSeismic (version 1.0.5)

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 Rlist 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

http://pubs.usgs.gov/of/2005/1438/pdf/OFR-1438.pdf{Seismic Noise Analysis System Using Power Spectral Density Probability Density Functions} (McNamara and Boaz 2005)

See Also

McNamaraPSD, psdList2NoiseMatrix, psdPlot, psdStatistics,

Examples

# 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')
}