# crossSpectrum: Cross-Spectral Analysis In IRISSeismic: Classes and Methods for Seismic Data 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

 ```1 2 3 4``` ```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:

 ```1 2 3``` ``` 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` cross-periodogram for ts1 and ts2 `Pxy` periodogram for ts2

## Author(s)

Jonathan Callahan jonathan@mazamascience.com

## References

`McNamaraPSD`
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43``` ```## 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) ```