crossSpectrum {IRISSeismic} | R Documentation |
Cross-Spectral Analysis
Description
The crossSpectrum() function is based on R's spec.pgram() function and attempts to provide complete results of cross-spectral FFT analysis in a programmer-friendly fashion.
Usage
crossSpectrum(x, spans = NULL, kernel = NULL, taper = 0.1,
pad = 0, fast = TRUE,
demean = FALSE, detrend = TRUE,
na.action = stats::na.fail)
Arguments
x |
multivariate time series |
spans |
vector of odd integers giving the widths of modified Daniell smoothers to be used to smooth the periodogram |
kernel |
alternatively, a kernel smoother of class "tskernel" |
taper |
specifies the proportion of data to taper. A split cosine bell taper is applied to this proportion of the data at the beginning and end of the series |
pad |
proportion of data to pad. Zeros are added to the end of the series to increase its length by the proportion pad |
fast |
logical. if TRUE, pad the series to a highly composite length |
demean |
logical. If TRUE, subtract the mean of the series |
detrend |
logical. If TRUE, remove a linear trend from the series. This will also remove the mean |
na.action |
NA action function |
Details
The multivariate timeseries passed in as the first argument should be a union of two separate timeseries with the same sampling rate created in the following manner:
ts1 <- ts(data1,frequency=sampling_rate) ts2 <- ts(data2,frequency=sampling_rate) x <- ts.union(ts1,ts2)
The crossSpectrum() function borrows most of its code from R's spec.pgram() function. It omits any plotting functionality and returns a programmer-friendly dataframe of all cross-spectral components generated during Fourier analysis for use in calculating transfer functions.
The naming of cross-spectral components is borrowed from the Octave version of MATLAB's pwelch() function.
Value
A dataframe with the following columns:
freq |
spectral frequencies |
spec1 |
'two-sided' spectral amplitudes for ts1 |
spec2 |
'two-sided' spectral amplitudes for ts2 |
coh |
magnitude squared coherence between ts1 and ts2 |
phase |
cross-spectral phase between ts1 and ts2 |
Pxx |
periodogram for ts1 |
Pyy |
periodogram for ts2 |
Pxy |
cross-periodogram for ts1 and ts2 |
Author(s)
Jonathan Callahan jonathan@mazamascience.com
References
Normalization of Power Spectral Density estimates
See Also
Examples
## Not run:
# Create a new IrisClient
iris <- new("IrisClient")
# Get seismic data
starttime <- as.POSIXct("2011-05-01", tz="GMT")
endtime <- starttime + 3600
st1 <- getDataselect(iris,"CI","PASC","00","BHZ",starttime,endtime)
st2 <- getDataselect(iris,"CI","PASC","10","BHZ",starttime,endtime)
tr1 <- st1@traces[[1]]
tr2 <- st2@traces[[1]]
# Both traces have a sampling rate of 40 Hz
sampling_rate <- tr1@stats@sampling_rate
ts1 <- ts(tr1@data,frequency=sampling_rate)
ts2 <- ts(tr2@data,frequency=sampling_rate)
# Calculate the cross spectrum
DF <- crossSpectrum(ts.union(ts1,ts2),spans=c(3,5,7,9))
# Calculate the transfer function
transferFunction <- DF$Pxy / DF$Pxx
transferAmp <- Mod(transferFunction)
transferPhase <- pracma::mod(Arg(transferFunction) * 180/pi,360)
# 2 rows
layout(matrix(seq(2)))
# Plot
plot(1/DF$freq,transferAmp,type='l',log='x',
xlab="Period (sec)",
main="Transfer Function Amplitude")
plot(1/DF$freq,transferPhase,type='l',log='x',
xlab="Period (sec)", ylab="degrees",
main="Transfer Function Phase")
# Restore default layout
layout(1)
## End(Not run)