eha: Evolutive Harmonic Analysis & Evolutive Power Spectral Analysis

Description

Evolutive Harmonic Analysis & Evolutive Power Spectral Analysis using the Thomson multitaper method (Thomson, 1982)

Usage

eha(dat,tbw=2,pad,fmin,fmax,step,win,demean=T,detrend=T,siglevel=0.90,
    sigID=F,ydir=1,output=0,pl=1,palette=6,centerZero=T,ncolors=100,xlab,ylab,
    genplot=2,verbose=T)

Arguments

dat: Stratigraphic series to analyze. First column should be location (e.g., depth), second column should be data value.
tbw: MTM time-bandwidth product (<=10)
pad: Pad with zeros to how many points? Must not factor into a prime number >23. Maximum number of points is 200,000.
fmin: Smallest frequency for analysis and plotting.
fmax: Largest frequency for analysis and plotting.
step: Step size for EHA window, in units of space or time.
win: Window size for EHA, in units of space or time.
demean: Remove mean from data series? (T or F)
detrend: Remove linear trend from data series? (T or F)
siglevel: Significance level for peak identification/filtering (0-1)
sigID: Identify signficant frequencies on power, amplitude, and probabilty plots. Only applies when one spectrum is calculated. (T or F)
ydir: Direction for y-axis in EHA plots (depth,height,time). -1 = values increase downwards (slower plotting), 1 = values increase upwards
output: Return output as new data frame? 0=no; 1=all results; 2=power; 3=amplitude; 4=probability; 5=significant frequencies (only for one spectrum); 6=significant frequencies and their probabilities (only for one spectrum)
pl: Plot logarithm of spectral power (1) or linear spectral power (2)?
palette: What color palette would you like to use? (1) rainbow, (2) grayscale, (3) blue, (4) red, (5) blue-white-red (if values are negative and positive, white is centered on zero), (6) viridis
centerZero: Center color scale on zero (use an equal number of postive and negative color divisions)? (T or F)
ncolors: Number of colors steps to use in palette.
xlab: Label for x-axis. Default = "Frequency"
ylab: Label for y-axis. Default = "Location"
genplot: Plotting options. 0= no plots; 1= power, amplitude, f-test, probability; 2=data series, power, amplitude, probability; 3= data series, power, normalized amplitude (maximum in each window normalized to unity), normalized amplitude filtered at specified siglevel; 4= data series, normalized power (maximum in each window normalized to unity), normalized amplitude (maximum in each window normalized to unity), normalized amplitude filtered at specified siglevel
verbose: Verbose output? (T or F)

Details

The power spectrum normalization approach applied here divides the Fourier coefficients by the number of points (npts) in the stratigraphic series, which is equivalent to dividing the power by (npts*npts). The (npts*npts) normalization has the convenient property whereby -- for an unpadded series -- the sum of the power in the positive frequencies is equivalent to half of variance; the other half of the variance is in the negative frequencies.

Note that the 'spec.mtm' function in package 'multitaper' (Rahim et al., 2014) is used for MTM spectrum estimation.

References

Thomson, D. J., 1982, Spectrum estimation and harmonic analysis, Proc. IEEE, 70, 1055-1096, doi:10.1109/PROC.1982.12433.

Examples

Run this code

## as an example, evaluate the modelA
data(modelA)

## interpolate to even sampling interval of 0.075 m
ex1=linterp(modelA, dt=0.075)
  
## perform EHA with a time-bandwidth parameter of 2, using an 7.95 meter window, 0.15 m step, 
## and pad to 1000 points
## set labels for plots (optional)
eha(ex1,tbw=2,win=7.95,step=0.15,pad=1000,xlab="Frequency (cycles/m)",ylab="Height (m)")

## for comparison generate spectrum for entire record, using time-bandwidth parameter of 3, and 
## pad to 5000 points
## start by making a new plot
pl(1)
eha(ex1,tbw=3,win=38,pad=5000,xlab="Frequency (cycles/m)")