IRISSeismic (version 1.0.5)

crossSpectrum: Cross-Spectral Analaysis

Description

The crossSpectrum() function is based on R's spec.pgram() function and attempts to provide complete results of cross-spectal 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

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
  cross-periodogram for ts1 and ts2
• Pxy
  periodogram for ts2

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.untion(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.

References

http://www.ltrr.arizona.edu/~dmeko/notes_6.pdf{Spectral Analysis -- Smoothed Periodogram Method}

http://octave-signal.sourcearchive.com/documentation/1.0.7/pwelch_8m-source.html{Octave pwelch() source code}

http://cran.r-project.org/web/packages/psd/vignettes/normalization.pdf{Normalization of Power Spectral Density estimates}

See Also

McNamaraPSD

Examples

# Create a new IrisClient
iris <- new("IrisClient", debug=TRUE)

# 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 smpling 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)