crossSpectrum {IRISSeismic}

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 [e-mail]

References

Octave pwelch() source code

Normalization of Power Spectral Density estimates

See Also

McNamaraPSD

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)