Apply McNamara PSD algorithm to a seismic signal
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.
psdList(st)
st |
a |
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)
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.
Jonathan Callahan jonathan@mazamascience.com
Seismic Noise Analysis System Using Power Spectral Density Probability Density Functions (McNamara and Boaz 2005)
## 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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.