psdList: Apply McNamara PSD algorithm to a seismic signal

In IRISSeismic: Classes and Methods for Seismic Data Analysis

Apply McNamara PSD algorithm to a seismic signal

Description

The psdList function subsets 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

Details

A Stream will be subset into segments depending upon the channel identifier (@stats@channel) 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)

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.

Author(s)

Jonathan Callahan [e-mail]

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)

IRISSeismic documentation built on June 8, 2025, 10:30 a.m.