Package 'astrochron'

Description

This software provides routines for astrochronologic testing, astronomical time scale construction, and time series analysis <doi:10.1016/j.earscirev.2018.11.015>. Also included are a range of statistical analysis and modeling routines that are relevant to time scale development and paleoclimate analysis.

Details

Note

Development of the 'astrochron' package is partially supported by the U.S. National Science Foundation and the Heising-Simons Foundation:

Leveraging the Geologic Record to Constrain Solar System Evolution, Earth-Moon Dynamics, Paleoclimate Change and Geological Time (Heising-Simons Foundation Award 2021-2797)

Collaborative Research: Improving the Late Cretaceous-Eocene geomagnetic polarity time scale by integrating the global magnetic anomaly record and astrochronology (U.S. National Science Foundation Award OCE 2051616)

CAREER: Deciphering the Beat of a Timeless Rhythm - The Future of Astrochronology (U.S. National Science Foundation Award EAR 1151438)

Collaborative Research: Evolution of the Climate Continuum - Late Paleogene to Present (U.S. National Science Foundation Award OCE 1003603)

TO CITE THIS PACKAGE IN PUBLICATIONS, PLEASE USE:

Meyers, S.R. (2014). Astrochron: An R Package for Astrochronology. https://cran.r-project.org/package=astrochron

Also cite the original research papers that document the relevant algorithms, as referenced on the help pages for specific functions.

Some examples (including R scripts) demonstrating the application of the astrochron package can be found in Meyers (2019), Ma et al. (2017), and Crampton et al. (2018), among numerous other scientific publications. Please also see examples in the help pages for specific functions.

Author(s)

Stephen Meyers

Maintainer: Stephen Meyers <[email protected]>

References

J.S. Campton, S.R. Meyers, R.A. Cooper, P.M Sadler, M. Foote, D. Harte, 2018, Pacing of Paleozoic macroevolutionary rates by Milankovitch grand cycles: Proceedings of the National Academy of Sciences, doi:10.1073/pnas.1714342115.

C. Ma, S.R. Meyers, B.B. Sageman, 2017 Theory of chaotic orbital variations confirmed by Cretaceous geological evidence: Nature v.542, 468-470, doi:10.1038/nature21402.

S.R. Meyers, 2019, Cyclostratigraphy and the problem of astrochronologic testing: Earth-Science Reviews v.190, 190-223, doi:10.1016/j.earscirev.2018.11.015.

Description

Calculate geochemical proxy accumulation rates using defined sedimentation rate and dry bulk density. Sedimentation rate and dry bulk density data are piecewise linearly interpolated to data locations, to estimate bulk accumulation rate.

Usage

accum(dat,sedrate=NULL,density=NULL,genplot=T,verbose=T)

Arguments

Examples

# Generate an example geochemical data set, using a stochastic model
    ex=ar1(npts=1001)
# the first column in 'ex' is depth (m), and we will take the second column to be CaCO3
#  rescale second column to range from 0-100% CaCO3
    m=100/(max(ex[2])-min(ex[2]))
    b=100-(m*max(ex[2]))
    ex[2]= m*ex[2] +b
    autoPlot(ex)
# generate example sedimentation rate history, with a stepwise tripling 
#  from 1m/kyr to 3m/kyr at 500 m
    sr=ex
    sr[1:500,2]=1
    sr[501:1001,2]=3
    autoPlot(sr)

# calculate accumulation rates
    res=accum(dat=ex,sedrate=sr,density=2.65)

Description

Anchor a floating astrochronology to a radioisotopic age. The floating astrochronology is centered on a given ('floating') time datum and assigned the 'anchored' age.

Usage

anchorTime(dat,time,age,timeDir=1,flipOut=F,genplot=T,verbose=T)

Arguments

Description

Generate AR(1) surrogates. Implement shuffling algorithm of Meyers (2012) if desired.

Usage

ar1(npts=1024,dt=1,mean=0,sdev=1,rho=0.9,shuffle=F,nsim=1,genplot=T,verbose=T)

Arguments

Details

These simulations use the random number generator of Matsumoto and Nishimura [1998]. If shuffle = T, the algorithm from Meyers (2012, pg. 11) is applied: (1) two sets of random sequences of the same length are generated, (2) the first random sequence is then sorted, and finally (3) the permutation vector of the sorted sequence is used to reorder the second random number sequence. This is done to guard against potential shortcomings in random number generation that are specific to spectral estimation.

References

S.R. Meyers, 2012, Seeing red in cyclic stratigraphy: Spectral noise estimation for astrochronology: Paleoceanography, v. 27, PA3328.

Description

Simulate a combined AR(1) + ETP signal, plot spectrum and confidence levels

Usage

ar1etp(etpdat=NULL,nsim=100,rho=0.9,wtAR=1,sig=90,tbw=2,padfac=5,ftest=F,fmax=0.1,
       speed=0.5,pl=2,graphfile=0)

Arguments

Details

Note: Setting wtAR=1 will provide equal variance contributions from the etp model and the ar1 model. More generally, set wtAR to the square root of the desired variance contribution (wtAR=0.5 will generate an AR1 model with variance that is 25% of the etp model). If you would like to exclusively evaluate the noise (no etp), set wtAR < 0.

Note: You may use the function etp to generate eccentricity-tilt-precession models.

References

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A.C.M., Levrard, B., 2004, A long term numerical solution for the insolation quantities of the Earth: Astron. Astrophys., Volume 428, 261-285.

See Also

getLaskar, and etp

Examples

## Not run: 

# run simulations using the default settings
ar1etp()

# compare with a second model:
# generate etp model spanning 0-2000 ka, with sampling interval of 5 ka.
ex1=etp(tmin=0,tmax=2000,dt=5)
# run simulations, with rho=-.7, and scaling noise to have 50
ar1etp(etpdat=ex1,rho=0.7,wtAR=sqrt(0.5))

 
## End(Not run)

Description

Arcsine transformation of stratigraphic series

Usage

arcsinT(dat,genplot=T,verbose=T)

Arguments

See Also

demean, detrend, divTrend, logT, prewhiteAR, and prewhiteAR1

Description

Generate an autoregressive moving-average time series model

Usage

armaGen(npts=1024,dt=1,m=0,std=1,rhos=c(0.9),thetas=c(0),genplot=T,verbose=T)

Arguments

Description

Calculate Average Spectral Misfit with Monte Carlo spectra simulations, as updated in Meyers et al. (2012).

Usage

asm(freq,target,fper=NULL,rayleigh,nyquist,sedmin=1,sedmax=5,numsed=50,
    linLog=1,iter=100000,output=F,genplot=T)

Arguments

Details

This function will caculate the Average Spectral Misfit between a data spectrum and astronomical target spectrum, following the approach outlined in Meyers and Sageman (2007), and the improvements of Meyers et al. (2012).

Value

A data frame containing: Sedimentation rate (cm/ka), ASM (cycles/ka), Null hypothesis significance level (0-100 percent), Number of astronomical terms fit.

References

S.R. Meyers and B.B. Sageman, 2007, Quantification of Deep-Time Orbital Forcing by Average Spectral Misfit: American Journal of Science, v. 307, p. 773-792.

S.R. Meyers, B.B. Sageman and M.A. Arthur, 2012, Obliquity forcing of organic matter accumulation during Oceanic Anoxic Event 2: Paleoceanography, 27, PA3212, doi:10.1029/2012PA002286.

See Also

eAsm, eAsmTrack, testPrecession, timeOpt, and timeOptSim

Examples

## These frequencies are from modelA (type '?astrochron' for more information). 
## They are for an 8 meter window, centered at 22 meters height. Units are cycles/m . 
freq <- c(0.1599833,0.5332776,1.5998329,2.6797201,3.2796575,3.8795948,5.5194235,6.5459830)
freq <- data.frame(freq)

## Rayleigh frequency in cycles/m
rayleigh <- 0.1245274

## Nyquist frequency in cycles/m
nyquist <- 6.66597

## orbital target in 1/ky. Predicted periods for 94 Ma (see Meyers et al., 2012)
target <- c(1/405.47,1/126.98,1/96.91,1/37.66,1/22.42,1/18.33)

## percent uncertainty in orbital target
fper=c(0.023,0.046,0.042,0.008,0.035,0.004)

asm(freq=freq,target=target,fper=fper,rayleigh=rayleigh,nyquist=nyquist,sedmin=0.5,sedmax=3,
    numsed=100,linLog=1,iter=100000,output=FALSE)

Description

Automatically plot and smooth specified stratigraphic data, versus location. Data are smoothed with a Gaussian kernel if desired.

Usage

autoPlot(dat,cols=NULL,dmin=NULL,dmax=NULL,vertical=T,ydir=NULL,nrows=NULL,plotype=1,
        smooth=0,xgrid=1,output=F,genplot=T,verbose=T)

Arguments

Description

Bandpass filter stratigraphic series using rectangular, Gaussian or tapered cosine (a.k.a. Tukey) window. This function can also be used to notch fiter a record (see examples).

Usage

bandpass(dat,padfac=2,flow=NULL,fhigh=NULL,win=0,alpha=3,p=0.25,demean=T,
         detrend=F,addmean=T,output=1,xmin=0,xmax=Nyq,genplot=T,verbose=T)

Arguments

Value

bandpassed stratigraphic series.

See Also

lowpass, noKernel, noLow, prewhiteAR, prewhiteAR1, and taner

Examples

# generate example series with periods of 405 ka, 100 ka, 40ka, and 20 ka, plus noise
ex=cycles(freqs=c(1/405,1/100,1/40,1/20),end=1000,dt=5,noisevar=.1)

# bandpass precession term using cosine-tapered window 
bandpass_ex <- bandpass(ex,flow=0.045,fhigh=0.055,win=2,p=.4)

# notch filter (remove) obliquity term using cosine-tapered window
#  if you'd like the final notch filtered record to be centered on the mean proxy 
#  value, set addmean=FALSE
notch_ex <- bandpass(ex,flow=0.02,fhigh=0.03,win=2,p=.4,addmean=FALSE)
notch_ex[2] <- ex[2]-notch_ex[2]
pl(2)
plot(ex,type="l",main="Eccentricity+Obliquity+Precession")
plot(notch_ex,type="l",main="Following application of obliquity notch filter")

Description

Determine the predicted precession and obliquity periods based on Berger et al. (1992). Values are determined by piecewise linear interpolation.

Usage

bergerPeriods(age,genplot=T)

Arguments

References

A. Berger, M.F. Loutre, and J. Laskar, 1992, Stability of the Astronomical Frequencies Over the Earth's History for Paleoclimate Studies: Science, v. 255, p. 560-566.

Description

bicoherence: Calculate bispectrum and bicoherence using WOSA method as detailed in Choudhury et al. (2008). Kim and Powers (1979) bicoherence normalization is used, with a conditioning factor to reduce the potential for spurious bicoherence peaks.

Usage

bicoherence(dat,overlap=50,segments=8,CF=0,CL=95,padfac=2,demean=T,detrend=T,
            taper=T,maxF=Nyq,output=0,genplot=T,color=1,id=NULL,logpwr=F,logbis=F,
            check=T,verbose=T)

Arguments

Details

This function accompanies the publication Sullivan et al. (in press 2023, PNAS): "Bicoherence is a higher-order statistic, which quantifies coupling of individual frequencies (as combination and difference tones). For three frequencies, f1, f2 and f3, high bicoherence indicates that f1 + f2 = f3, with the additional property that their phases are also coupled in the data series. We employ the bicoherence approach of Choudhury, Shah and Thornhill (2008), which uses Welch Overlapping Spectral Analysis (WOSA), as this method includes improvements to the standard approach that reduce the potential for spurious bicoherence peaks. Note that reliable bicoherence estimates require a substantial increase in bandwidth resolution ("smoothing")".

References

Choudhury, A.A.S., Shah, S.L., and Thornhill, N.F. (2008), Diagnosis of Process Nonlinearities and Valve Stiction: Data Driven Approaches: Springer, 284 pp.

Kim, Y.C., and Powers, E.J. (1979), Digital Bispectral Analysis and Its Applications to Nonlinear Wave Interactions: IEEE Transactions on Plasma Science, v. PS-7, 120-131.

Sullivan, N.B., Meyers, S.R., Levy, R.H., McKay, R.M., Golledge, N.R., Cortese, G. (in press 2023), Millennial-scale variability of the Antarctic Ice Sheet during the Early Miocene: Proceedings of the National Academy of Sciences.

Welch, P.D. (1967), The use of Fast Fourier Tranform for the estimation of power spectra: A method based on time averaging over short, modified periodograms: IEEE Transactions on Audio and Electroacoustics, AU-5 (2): 70-73.

Examples

## Not run: 
# Generate an example quadratic phase-coupled signal as in Choudhury et al. (2008, pg. 79)
    n = rnorm(500)
    t=1:500
    signal = sin(2*pi*0.12*t + pi/3) + sin(2*pi*0.18*t + pi/12) + sin(2*pi*0.3*t + 5*pi/12) + n
    ex=data.frame(cbind(t,signal))
    bicoherence(ex)
 
## End(Not run)

Description

'bioturb' is a function for simulating the bioturbed time series when bioturbation is modeled as a diffusive process. It implements the methodology outlined in Liu et al. (2021), which builds on the approaches of (Guinasso and Schink, 1975), Goreau (1977), and Goreau (1980). Given the bioturbation parameters, an impulse response function is first calculated using function 'impulseResponse'. This function is then convolved with the input time series (true signal) to output the bioturbed time series. Note that the input true signal 'dat' and impulse response function are interpolated to the same resolution before convolution.

Usage

bioturb(dat, G, ML, v, output=1, genplot=TRUE, verbose=TRUE)

Arguments

References

Guinasso, N.L. and Schinck, D.R., 1975, Quantitative estimates of biological mixing rates in abyssal sediments, J. Geophys. Res., 80, 3032-3043.

Goreau, T.J., 1977, Quantitative effects of sediment mixing on stratigraphy and biogeochemistry: a signal theory approach, Nature, 256, 730-732.

Goreau, T.J., 1980, Frequency sensitivity of the deep-sea climatic record, Nature, 287, 620-622.

Liu, H., Meyers, S.R., and Marcott, S.A., 2021, Unmixing dee-sea paleoclimate records: A study on bioturbation effects through convolution and deconvolution, Earth and Planetary Science Letters, 564, 116883.

Examples

# as a test series, use the three dominant precession terms from Berger et al. (1992)
ex1=cycles()

# mix it
res1 <- bioturb(ex1, G=4, ML=10, v=1, genplot = TRUE)

# un-mix it
res2=unbioturb(res1, G=4, ML=10, v=1, genplot = TRUE)

pl()
plot(ex1,type="l",main="black=signal, blue=bioturbated, red=unbioturbated",lwd=3)
lines(res2,col="red")
lines(res1,col="blue")

Description

Calculate eccentricity, obliquity and precession periods in ka, given g, s and k in arcsec/yr.

Usage

calcPeriods(g,s=NULL,k,output=1)

Arguments

Examples

# calculate eccentricity and precession periods for one set of g's and k
gVal= c(5.579378,7.456665,17.366595,17.910194,4.257564)
kVal= 50.475838
calcPeriods(g=gVal, k=kVal)

# calculate eccentricity and precession periods for three sets of g's and k
gVal= matrix(c(5.579378,7.456665,17.366595,17.910194,4.257564,5.494302,7.452619,
      17.480760,18.348310,4.257451,5.531285,7.456848,17.320480,17.912240,4.257456),
      nrow=3,ncol=5,byrow=TRUE)

kVal= c(50.475838,51.280910,85.790450)

calcPeriods(g=gVal, k=kVal)

Description

Bind two vectors together and return result as a data frame. Alternatively, extract specified columns from a data frame, bind them together, and return result as a data frame.

Usage

cb(a,b)

Arguments

Examples

# example dataset
x<-rnorm(100)
dim(x)<-c(10,10)
x<-data.frame(x)

# bind two columns
cb(x[1],x[5])

# bind five columns
cb(x,c(1,2,4,7,9))

Description

Create non-linear response by clipping stratigraphic series below a threshold value. Alternatively, mute response below a threshold value using a contant divisor. Both approaches will enhance power in modulator (e.g., eccentricity) and diminish power the carrier (e.g., precession).

Usage

clipIt(dat,thresh=NULL,clipval=NULL,clipdiv=NULL,genplot=T,verbose=T)

Arguments

Description

Adjust spectrum confidence levels for multiple comparisons, using the Bonferroni correction

Usage

confAdjust(spec,npts,dt,tbw=3,ntap=5,flow=NULL,fhigh=NULL,output=T,
    xmin=df,xmax=NULL,pl=1,genplot=T,verbose=T)

Arguments

Details

Multiple testing is a common problem in the evaluation of power spectrum peaks (Vaughan et al., 2011; Crampton et al., PNAS). To address the issue of multiple testing, a range of approaches have been advocated. This function will conduct an assessment using the Bonferroni correction, which is the simplest, and also the most conservative, of the common approaches (it is overly pessimistic).

If one is exclusively concerned with particular frequency bands a priori (e.g., those associated with Milankovitch cycles), the statistical power of the method can be improved by restricting the analysis to those frequency bands (use options 'flow' and 'fhigh').

Application of multiple testing corrections does not guarantee that the spectral background is appropriate. To address this issue, carefully examine the fit of the spectral background, and also conduct simulations with the function testBackground.

References

J.S. Campton, S.R. Meyers, R.A. Cooper, P.M Sadler, M. Foote, D. Harte, 2018, Pacing of Paleozoic macroevolutionary rates by Milankovitch grand cycles: Proceedings of the National Academy of Sciences, doi:10.1073/pnas.1714342115.

S. Vaughan, R.J. Bailey, and D.G. Smith, 2011, Detecting cycles in stratigraphic data: Spectral analysis in the presence of red noise. Paleoceanography 26, PA4211, doi:10.1029/2011PA002195.

See Also

testBackground,multiTest, lowspec, and periodogram

Examples

# generate example series with periods of 400 ka, 100 ka, 40 ka and 20 ka
ex = cycles(freqs=c(1/400,1/100,1/40,1/20),start=1,end=1000,dt=5)

# add AR1 noise
noise = ar1(npts=200,dt=5,sd=.5)
ex[2] = ex[2] + noise[2]

# first, let's use mtm with conventional AR1 background
spec=mtm(ex,padfac=1,ar1=TRUE,output=1)

# when blindly prospecting for cycles, it is necessary to consider all of the 
#  observed frequencies in the test
confAdjust(spec,npts=200,dt=5,tbw=3,ntap=5,output=FALSE)

# if, a priori, you are only concerned with the Milankovitch frequency bands, 
#  restrict your analysis to those bands (as constrained by available sedimentation
#  rate estimates and the frequency resolution of the spectrum). in the example below, 
#  the mtm bandwidth resolution is employed to search frequencies nearby the 
#  Milankovitch-target periods.
flow=c((1/400)-0.003,(1/100)-0.003,(1/41)-0.003,(1/20)-0.003)
fhigh=c((1/400)+0.003,(1/100)+0.003,(1/41)+0.003,(1/20)+0.003)
confAdjust(spec,npts=200,dt=5,tbw=3,ntap=5,flow=flow,fhigh=fhigh,output=FALSE)

# now try with the lowspec method. this uses prewhitening, so it has one less data point.
spec=lowspec(ex,padfac=1,output=1)
flow=c((1/400)-0.003015075,(1/100)-0.003015075,(1/41)-0.003015075,(1/20)-0.003015075)
fhigh=c((1/400)+0.003015075,(1/100)+0.003015075,(1/41)+0.003015075,(1/20)+0.003015075)
confAdjust(spec,npts=199,dt=5,tbw=3,ntap=5,flow=flow,fhigh=fhigh,output=FALSE)

# for comparison...
confAdjust(spec,npts=199,dt=5,tbw=3,ntap=5,output=FALSE)

Description

Apply a constant sedimentation rate model to transform a spatial series to temporal series.

Usage

constantSedrate(dat,sedrate,begin=0,timeDir=1,genplot=T,verbose=T)

Arguments

Description

Convolution through Fast Fourier Transform

Usage

conv_fft(x, y, index, dt)

Arguments

Details

Function 'conv_fft' is used by function 'bioturb'. x and y do not need to be of the same length.

Value

z Convolved output series. length(z) = length(x)

Description

Apply a "percent-tapered" cosine taper (a.k.a. Tukey window) to a stratigraphic series.

Usage

cosTaper(dat,p=.25,rms=T,demean=T,detrend=F,genplot=T,verbose=T)

Arguments

See Also

dpssTaper, gausTaper, and hannTaper

Description

Make a time series with specified harmonic components and noise

Usage

cycles(freqs=NULL,phase=NULL,amp=NULL,start=0,end=499,dt=1,noisevar=0,genplot=T,
        verbose=T)

Arguments

Value

modeled time series.

Examples

## test signal on pg 38 of Choudhury, Shah, and Thornhill (2008)
freqs=c(0.12,0.18,0.30,0.42)
phase=c(-pi/3,-pi/12,-pi/4,-3*pi/8)
amp=c(1,1,1,1)

cycles(freqs,phase,amp,start=0,end=4095,dt=1,noisevar=0.2)

Description

Wiener Deconvolution through Fast Fourier Transform.

Usage

deconv(x, y, index, dt, pt = 0.2, wiener = TRUE)

Arguments

Details

function 'deconv' is used by function 'unbioturb'. A cosine taper is applied to remove edge effects. The signal-to-noise ratio is chosen to be 0.05, and gamma is also chosen to be 0.05. x and y do not need to be of the same length. For additional information see Liu et al. (2021)

Value

z deconvolved/un-bioturbed series. length(z) = length(x)

References

Liu, H., Meyers, S.R., and Marcott, S.A., 2021, Unmixing dee-sea paleoclimate records: A study on bioturbation effects through convolution and deconvolution, Earth and Planetary Science Letters, 564, 116883.

Description

Interactively delete points in x,y plot.

Usage

delPts(dat,del=NULL,cols=c(1,2),ptsize=1,xmin=NULL,xmax=NULL,ymin=NULL,ymax=NULL,
       plotype=1,genplot=T,verbose=T)

Arguments

See Also

idPts, iso, trim and trimAT

Description

Remove mean value from stratigraphic series

Usage

demean(dat,genplot=T,verbose=T)

Arguments

See Also

arcsinT, detrend, divTrend, logT, prewhiteAR, and prewhiteAR1

Description

Remove linear trend from stratigraphic series

Usage

detrend(dat,output=1,genplot=T,verbose=T)

Arguments

See Also

arcsinT, demean, divTrend, logT, prewhiteAR, and prewhiteAR1

Description

Model differential accumulation. The input variable (e.g., insolation, proxy value) is rescaled to sedimentation rate curve varying from sedmin to sedmax. Input series must be evenly sampled in time.

Usage

diffAccum(dat,sedmin=0.01,sedmax=0.02,dir=1,genplot=T,verbose=T)

Arguments

Examples

# generate model with one 20 ka cycle
ex <- cycles(1/20)

diffAccum(ex)

Description

Divide data series value by linear trend observed in stratigraphic series

Usage

divTrend(dat,genplot=T,verbose=T)

Arguments

See Also

arcsinT, demean, detrend, logT, prewhiteAR, and prewhiteAR1

Description

Apply a single Discrete Prolate Spheroidal Sequence (DPSS) taper to a stratigraphic series

Usage

dpssTaper(dat,tbw=1,num=1,rms=T,demean=T,detrend=F,genplot=T,verbose=T)

Arguments

See Also

cosTaper, gausTaper, and hannTaper

Description

Calculate Evolutive Average Spectral Misfit with Monte Carlo spectra simulations, as updated in Meyers et al. (2012).

Usage

eAsm(spec,siglevel=0.9,target,fper=NULL,rayleigh,nyquist,sedmin=1,sedmax=5,
      numsed=50,linLog=1,iter=100000,ydir=1,palette=2,output=4,genplot=F)

Arguments

Details

Please see function asm for details.

References

S.R. Meyers and B.B. Sageman, 2007, Quantification of Deep-Time Orbital Forcing by Average Spectral Misfit: American Journal of Science, v. 307, p. 773-792.

S.R. Meyers, 2012, Seeing Red in Cyclic Stratigraphy: Spectral Noise Estimation for Astrochronology: Paleoceanography, 27, PA3228, doi:10.1029/2012PA002307.

S.R. Meyers, B.B. Sageman and M.A. Arthur, 2012, Obliquity forcing of organic matter accumulation during Oceanic Anoxic Event 2: Paleoceanography, 27, PA3212, doi:10.1029/2012PA002286.

See Also

asm, eAsmTrack, eha, testPrecession, timeOpt, and timeOptSim

Examples

## Not run: 
# use modelA as an example
data(modelA)

# interpolate to even sampling interval
modelAInterp=linterp(modelA)

# perform EHA analysis, save harmonic F-test confidence level results to 'spec'
spec=eha(modelAInterp,win=8,step=2,pad=1000,output=4)

# perform Evolutive Average Spectral Misfit analysis, save results to 'res'
res=eAsm(spec,target=c(1/405.47,1/126.98,1/96.91,1/37.66,1/22.42,1/18.33),rayleigh=0.1245274,
         nyquist=6.66597,sedmin=0.5,sedmax=3,numsed=100,siglevel=0.8,iter=10000,output=4)

# identify minimum Ho-SL in each record and plot
pl(1)
eAsmTrack(res[1],threshold=0.05)

# extract Ho-SL result at 18.23 m
HoSL18.23=extract(res[1],get=18.23,pl=1)

# extract ASM result at 18.23 m
asm18.23=extract(res[2],get=18.23,pl=0)
 
## End(Not run)

Description

Track ASM null hypothesis significance level minima in eASM results.

Usage

eAsmTrack(res,threshold=.5,ydir=-1,genplot=T,verbose=T)

Arguments

Details

Please see function eAsm for details.

See Also

asm, eAsm, and eha

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

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.

See Also

extract, lowspec, mtmAR, mtmML96, periodogram, trackFreq and traceFreq

Examples

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

Description

eTimeOpt: Evolutive implementation of TimeOpt (Meyers, 2015; Meyers, 2019).

Usage

eTimeOpt(dat,win=dt*100,step=dt*10,sedmin=0.5,sedmax=5,numsed=100,linLog=1,
    limit=T,fit=1,fitModPwr=T,flow=NULL,fhigh=NULL,roll=NULL,targetE=NULL,targetP=NULL,
    detrend=T,ydir=1,palette=6,ncolors=100,output=1,genplot=T,check=T,verbose=1)

Arguments

References

S.R. Meyers, 2015, The evaluation of eccentricity-related amplitude modulations and bundling in paleoclimate data: An inverse approach for astrochronologic testing and time scale optimization: Paleoceanography, v.30, 1625-1640.

S.R. Meyers, 2019, Cyclostratigraphy and the problem of astrochronologic testing: Earth-Science Reviews v.190, 190-223.

See Also

tracePeak,trackPeak,timeOpt,timeOptSim, and eTimeOptTrack

Examples

## Not run: 
# generate a test signal with precession and eccentricity
ex=cycles(freqs=c(1/405.6795,1/130.719,1/123.839,1/98.86307,1/94.87666,1/23.62069, 
1/22.31868,1/19.06768,1/18.91979),end=4000,dt=5)

# convert to meters with a linearly increasing sedimentation rate from 0.01 m/kyr to 0.03 m/kyr
ex=sedRamp(ex,srstart=0.01,srend=0.03)

# interpolate to median sampling interval
ex=linterp(ex)

# evaluate precession & eccentricity power, and precession modulations
res=eTimeOpt(ex,win=20,step=1,fit=1,output=1)

# extract the optimal fits for the power optimization
sedrates=eTimeOptTrack(res[2])

# extract the optimal fits for the envelope*power optimization
sedrates=eTimeOptTrack(res[3])

# you can also interactively track the results using functions 'trackPeak' and 'tracePeak'
#  evaluate the results from the power optimization
sedrates=tracePeak(res[2])
sedrates=trackPeak(res[2])

#  evaluate the results from the envelope*power optimization
sedrates=tracePeak(res[3])
sedrates=trackPeak(res[3])

# evaluate precession & eccentricity power, and short-eccentricity modulations
eTimeOpt(ex,win=20,step=1,fit=2,output=0)

## End(Not run)

Description

Track eTimeOpt r2 maxima.

Usage

eTimeOptTrack(res,threshold=0,ydir=-1,genplot=T,verbose=T)

Arguments

Details

Please see function eTimeOpt for details.

References

S.R. Meyers, 2015, The evaluation of eccentricity-related amplitude modulations and bundling in paleoclimate data: An inverse approach for astrochronologic testing and time scale optimization: Paleoceanography, v.30, 1625-1640.

S.R. Meyers, 2019, Cyclostratigraphy and the problem of astrochronologic testing: Earth-Science Reviews v.190, 190-223.

See Also

timeOpt, and eTimeOpt

Examples

## Not run: 
# generate a test signal with precession and eccentricity
ex=cycles(freqs=c(1/405.6795,1/130.719,1/123.839,1/98.86307,1/94.87666,1/23.62069, 
1/22.31868,1/19.06768,1/18.91979),end=4000,dt=5)

# convert to meters with a linearly increasing sedimentation rate from 0.01 m/kyr to 0.03 m/kyr
ex=sedRamp(ex,srstart=0.01,srend=0.03)

# interpolate to median sampling interval
ex=linterp(ex)

# evaluate precession & eccentricity power, and precession modulations
res=eTimeOpt(ex,win=20,step=1,fit=1,output=1)

# extract the optimal fits for the power optimization
sedrates=eTimeOptTrack(res[2])

# extract the optimal fits for the envelope*power optimization
sedrates=eTimeOptTrack(res[3])

# you can also interactively track the results using functions 'trackPeak' and 'tracePeak'
#  evaluate the results from the power optimization
sedrates=trackPeak(res[2])
sedrates=tracePeak(res[2])

#  evaluate the results from the envelope*power optimization optimization
sedrates=trackPeak(res[3])
sedrates=tracePeak(res[3])

## End(Not run)

Description

Calculate eccentricity-tilt-precession time series using the theoretical astronomical solutions. By default, the Laskar et al. (2004) solutions will be downloaded. Alternatively, one can specify the astronomical solution.

Usage

etp(tmin=NULL,tmax=NULL,dt=1,eWt=1,oWt=1,pWt=1,esinw=T,solution=NULL,standardize=T,
     genplot=T,verbose=T)

Arguments

Details

Note: If you plan to repeatedly execute the etp function, it is advisable to download the astronomical solution once using the function getLaskar.

Note: It is common practice to construct ETP models that have specified variance ratios (e.g., 1:1:1 or 1:0.5:0.5) for eccentricity, obliquity and precession. In order to construct such models, it is necessary to choose 'standardize=T', and to set the individual weights (eWt, oWt, pWt) to the square root of the desired variance contribution.

Value

Eccentricity + tilt + precession.

References

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A.C.M., Levrard, B., 2004, A long term numerical solution for the insolation quantities of the Earth: Astron. Astrophys., Volume 428, 261-285.

Laskar, J., Fienga, A., Gastineau, M., Manche, H., 2011, La2010: A new orbital solution for the long-term motion of the Earth: Astron. Astrophys., Volume 532, A89.

Laskar, J., Gastineau, M., Delisle, J.-B., Farres, A., Fienga, A.: 2011, Strong chaos induced by close encounters with Ceres and Vesta: Astron. Astrophys., Volume 532, L4.

See Also

getLaskar

Examples

## Not run: 
# create an ETP model from 10000 ka to 20000 ka, with a 5 ka sampling interval
# this will automatically download the astronomical solution
ex=etp(tmin=10000,tmax=20000,dt=5)

# alternatively, download the astronomical solution first
ex2=getLaskar()
ex=etp(tmin=10000,tmax=20000,dt=5,solution=ex2)
 
## End(Not run)

Description

Extract record from EHA time-frequency output or eAsm output: Use interactive graphical interface to identify record.

Usage

extract(spec,get=NULL,xmin=NULL,xmax=NULL,ymin=NULL,ymax=NULL,h=6,w=4,ydir=1,pl=0,
        ncolors=100,genplot=T,verbose=T)

Arguments

See Also

eha

Description

Flip the stratigraphic order of your data series (e.g., convert stratigraphic depth series to height series, relative to a defined datum.)

Usage

flip(dat,begin=0,genplot=T,verbose=T)

Arguments

Description

Convert record of local spatial frequency (from EHA) to sedimentation rate curve

Usage

freq2sedrate(freqs,period=NULL,ydir=1,genplot=T,verbose=T)

Arguments

Description

Apply a Gaussian taper to a stratigraphic series

Usage

gausTaper(dat,alpha=3,rms=T,demean=T,detrend=F,genplot=T,verbose=T)

Arguments

References

Harris, 1978, On the use of windows for harmonic analysis with the discrete Fourier transform: Proceedings of the IEEE, v. 66, p. 51-83.

See Also

cosTaper, dpssTaper, and hannTaper

Description

Query R for color information.

Usage

getColor(color)

Arguments

Description

Download data file from astrochron server.

Usage

getData(dat="1262-a*")

Arguments

Description

Download Laskar et al. (2004, 2011a, 2011b) astronomical solutions.

Usage

getLaskar(sol="la04",verbose=T)

Arguments

Details

la04 : three columns containing precession angle, obliquity, and eccentricity of Laskar et al. (2004)

la10a : one column containing the la10a eccentricity solution of Laskar et al. (2011a)

la10b : one column containing the la10b eccentricity solution of Laskar et al. (2011a)

la10c : one column containing the la10c eccentricity solution of Laskar et al. (2011a)

la10d : one column containing the la10d eccentricity solution of Laskar et al. (2011a)

la11 : one column containing the la11 eccentricity solution of Laskar et al. (2011b; please also cite 2011a)

insolation : one column containing insolation at 65 deg North (W/m^2) during summer solstice, from Laskar et al. (2004)

References

J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, and B. Levrard, B., 2004, A long term numerical solution for the insolation quantities of the Earth: Astron. Astrophys., Volume 428, 261-285.

Laskar, J., Fienga, A., Gastineau, M., Manche, H., 2011a, La2010: A new orbital solution for the long-term motion of the Earth: Astron. Astrophys., Volume 532, A89.

Laskar, J., Gastineau, M., Delisle, J.-B., Farres, A., Fienga, A.: 2011b, Strong chaos induced by close encounters with Ceres and Vesta, Astron: Astrophys., Volume 532, L4.

Description

Apply a Hann (Hanning) taper to a stratigraphic series

Usage

hannTaper(dat,rms=T,demean=T,detrend=F,genplot=T,verbose=T)

Arguments

See Also

cosTaper, dpssTaper, and gausTaper

Description

Execute 'head' function, with column numbers indicated for each variable. (useful for functions such as 'autopl')

Usage

headn(dat)

Arguments

Description

Calculate instantaneous amplitude (envelope) and phase via Hilbert Transform of stratigraphic series

Usage

hilbert(dat,phase=F,padfac=2,demean=T,detrend=F,output=T,addmean=F,genplot=T,check=T,
        verbose=T)

Arguments

Examples

# generate example series with 3 precession terms and noise
ex <- cycles(noisevar=.0004,dt=5)
# bandpass precession terms using cosine-tapered window 
res_ex <- bandpass(ex,flow=0.038,fhigh=0.057,win=2,p=.4)
# hilbert transform
hil_ex <- hilbert(res_ex)

Description

Interactively identify points in x,y plot.

Usage

idPts(dat1,dat2=NULL,ptsize=0.5,xmin=NULL,xmax=NULL,ymin=NULL,ymax=NULL,
       logx=F,logy=F,plotype=1,annotate=1,iso=F,output=1,verbose=T)

Arguments

See Also

delPts, iso, trim and trimAT

Description

An implementation of the Imbrie and Imbrie (1980) ice sheet model

Usage

imbrie(insolation=NULL,Tm=17,b=0.6,times=NULL,initial=0,burnin=100,standardize=T,
       output=T,genplot=1,verbose=T)

Arguments

Details

This function will implement the ice volume model of Imbrie and Imbrie (1980), following the conventions of Paillard et al. (1996).

When using the 'times' vector, consider the following example:

times= c(500,1000)

Tm=c(15,5)

b=c(0.6,0.3)

In this case, a Tm of 15 (b of 0.6) will be applied to model from 0-500 ka, and a Tm of 5 (b of 0.3) will be applied to model 500-1000 ka.

References

Imbrie, J., and Imbrie, J.Z., (1980), Modeling the Climatic Response to Orbital Variations: Science, v. 207, p. 943-953.

Lisiecki, L. E., and M. E. Raymo, 2005, A Pliocene-Pleistocene stack of 57 globally distributed benthic d18O records, Paleoceanography, 20, PA1003, doi:10.1029/2004PA001071.

Paillard, D., L. Labeyrie and P. Yiou, (1996), Macintosh program performs time-series analysis: Eos Trans. AGU, v. 77, p. 379.

Examples

## Not run: 
# make a very simple forcing (on/off)
forcing=cycles(0,end=300)
forcing[50:150,2]=1
plot(forcing,type="l")

# use this forcing to drive the imbrie ice model
# set b=0, Tm = 1
imbrie(forcing,b=0,Tm=1,output=F)

# let's view the evolution of the ice sheet
imbrie(forcing,b=0,Tm=1,output=F,genplot=2)

# now increase the response time
imbrie(forcing,b=0,Tm=10,output=F,genplot=2)

# now model slow growth, fast decay
imbrie(forcing,b=0.5,Tm=10,output=F,genplot=2)

# now make a 100 ka cyclic forcing
forcing=cycles(1/100,end=300)
imbrie(forcing,b=0,Tm=1,output=F,genplot=2)
imbrie(forcing,b=0,Tm=10,output=F,genplot=2)
imbrie(forcing,b=0.5,Tm=10,output=F,genplot=2)
# show burn-in
imbrie(forcing,b=0.5,Tm=10,output=F,genplot=2,burnin=0)

# now examine Malutin Milankovitch's hypothesis: 65 deg N, summer solstice
imbrie(b=0.5,Tm=10,output=F,genplot=2,burnin=900)

# use the ice model output to make a synthetic stratigraphic section
res=imbrie(b=0.5,T=10,output=T,genplot=1,burnin=100)
synthStrat(res,clip=F)

# generate ice model for last 5300 ka, using 65 deg. N insolation, 21 June
# allow b and Tm values to change as in Lisiecki and Raymo (2005):
insolation=getLaskar("insolation")
insolation=iso(insolation,xmin=0,xmax=5300)
#  b is 0.3 from 5300 to 3000 ka, then linearly increases to 0.6 between 3000 and 1500 ka.
#  b is 0.6 from 1500 ka to present.
set_b=linterp(cb(c(1500,3000),c(0.6,0.3)),dt=1)
set_b=rbind(set_b,c(5400,0.3))
#  Tm is 5 ka from 5300 to 3000 ka, then linearly increases to 15 ka between 3000 and 1500 ka.
#  Tm is 15 ka from 1500 ka to present.
set_Tm=linterp(cb(c(1500,3000),c(15,5)),dt=1)
set_Tm=rbind(set_Tm,c(5400,5))
# now run model
ex=imbrie(insolation=insolation,Tm=set_Tm[,2],b=set_b[,2],times=set_b[,1])
# time-frequency analysis of model result
eha(ex,fmax=0.1,win=500,step=10,pad=5000,genplot=4,pl=2)

 
## End(Not run)

Description

Calculate the analytical response function from an impulse forcing using the 1-D advection diffusion model proposed in Schink and Guinasso (1975) to model the bioturbation impact on climate time series.

Usage

impulseResponse(G, ML = NULL, v = NULL, nt = 500, genplot = FALSE,
  verbose = FALSE)

Arguments

Value

fc Impulse response function

References

Guinasso, N.L. and Schinck, D.R., 1975, Quantitative estimates of biological mixing rates in abyssal sediments, J. Geophys. Res., 80, 3032-3043.

Goreau, T.J., 1977, Quantitative effects of sediment mixing on stratigraphy and biogeochemistry: a signal theory approach, Nature, 256, 730-732.

Goreau, T.J., 1980, Frequency sensitivity of the deep-sea climatic record, Nature, 287, 620-622.

Liu, H., Meyers, S.R., and Marcott, S.A., 2021, Unmixing dee-sea paleoclimate records: A study on bioturbation effects through convolution and deconvolution, Earth and Planetary Science Letters, 564, 116883.

Examples

G <- 4
ML <- 10
v <- 1
# take a look at the IRF 
impulseResponse(G=4, ML = 10, v = 1, genplot = TRUE)

Description

Generate map that show difference in insolation between two times using Laskar et al. (2004) insolation calculations. This is a wrapper function to calculate insolation using the palinsol package, following astrochron syntax.

Usage

insoDiff(t1,t2,S0=1365,verbose=T)

Arguments

Details

This wrapper function generates maps with insolation calculations based on the Laskar et al. (2004) astronomical parameters, using the 'palinsol' package (Crucifix, 2016).

References

Crucifix, M., 2016, palinsol: Insolation for Palaeoclimate Studies. R package version 0.93: https://CRAN.R-project.org/package=palinsol

J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, and B. Levrard, B., 2004, A long term numerical solution for the insolation quantities of the Earth: Astron. Astrophys., Volume 428, 261-285.

Description

Generate insolation maps using Laskar et al. (2004) insolation calculations. This is a wrapper function to calculate insolation using the palinsol package, following astrochron syntax.

Usage

insoMap(t=0,pick=0,S0=1365,output=F,verbose=T)

Arguments

Details

This wrapper function generates maps with insolation calculations based on the Laskar et al. (2004) astronomical parameters, using the 'palinsol' package (Crucifix, 2016).

Option 'pick' will allow you to extract insolation by (1) latitude or (2) season.

References

Crucifix, M., 2016, palinsol: Insolation for Palaeoclimate Studies. R package version 0.93: https://CRAN.R-project.org/package=palinsol

J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, and B. Levrard, B., 2004, A long term numerical solution for the insolation quantities of the Earth: Astron. Astrophys., Volume 428, 261-285.

Description

Laskar et al. (2004) insolation calculations. This is a wrapper function to calculate insolation using the palinsol package, following astrochron syntax.

Usage

insoSeries(tmin=0,tmax=1000,dt=1,opt=1,long=90,lat=65,threshold=400,l1=0,l2=70,S0=1365,
         genplot=TRUE,verbose=TRUE)

Arguments

Details

This wrapper function performs four different types of insolation calculations based on the Laskar et al. (2004) astronomical parameters, using the 'palinsol' package (Crucifix, 2016).

The following 'palinsol' functions are used:

opt 1: Insol - Calculate Mean Insolation for Given Day and Latitude

opt 2: calins - Calculate Caloric Insolation for Given Latitude (Milankovitch; Berger, 1978)

opt 3: Insol_l1l2 - Calculate Integrated Insolation for Given Latitude (Berger et al., 2010)

opt 4: thrins - Calculate Integrated Insolation for Given Latitude, with Threshold (Huybers and Tziperman, 2008)

la04 - Calculate astronomical parameters (Laskar et al., 2004)

Please see those functions for further details.

If opt 1 is selected: specify 'long' and 'lat'

If opt 2 is selected: specify 'lat'

If opt 3 is selected: specify 'lat', 'l1' and 'l2'

If opt 4 is selected: specify 'lat' and 'threshold'

APPROXIMATE CONVENTIONAL DATES FOR A GIVEN MEAN LONGITUDE (long):

0: 21 March (equinox)

30: 21 April

60: 21 May

90: 21 June (solstice)

120: 21 July

150: 21 August

180: 21 September (equinox)

210: 21 October

240: 21 November

270: 21 December (solstice)

300: 21 January

330: 21 February

References

Berger, A., 1978, Long-term variations of caloric insolation resulting from the earth's orbital elements: Quaternary Research, 9, 139 - 167.

Berger, A., Loutre, M.F., Yin Q., 2010, Total irradiation during any time interval of the year using elliptic integrals: Quaternary Science Reviews, 29, 1968 - 1982, doi:10.1016/j.quascirev.2010.05.007

Crucifix, M., 2016, palinsol: Insolation for Palaeoclimate Studies. R package version 0.93: https://CRAN.R-project.org/package=palinsol

Huybers, P., Tziperman, E., 2008, Integrated summer insolation forcing and 40,000-year glacial cycles: The perspective from an ice-sheet/energy-balance model: Paleoceanography, 23.

J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, and B. Levrard, B., 2004, A long term numerical solution for the insolation quantities of the Earth: Astron. Astrophys., Volume 428, 261-285.

Description

Determine the total power within a given bandwidth, and also the ratio of this power to the total power in the spectrum (or up to a specified frequency). If bandwidth is not specified, generate interactive plots for bandwidth selection. For use with the function eha, integratePower can process spectrograms (time-frequency) or single spectra.

Usage

integratePower(spec,flow=NULL,fhigh=NULL,fmax=NULL,unity=F,f0=T,xmin=NULL,
               xmax=NULL,ymin=NULL,ymax=NULL,npts=NULL,pad=NULL,ydir=1,
               palette=6,ncolors=100,h=6,w=9,ln=F,genplot=T,verbose=T)

Arguments

Details

Depending on the normalization used, you may want to preprocess the power spectra prior to integration.

See Also

eha

Examples

# generate etp signal over past 10 Ma
ex=etp(tmax=10000)

# evolutive power
pwr=eha(ex,win=500,fmax=.1,pad=2000,output=2,pl=2)

# integrate power from main obliquity term
integratePower(pwr,flow=0.02,fhigh=0.029,npts=501,pad=2000)

Description

Isolate a section of a uni- or multi-variate stratigraphic data set for further analysis

Usage

iso(dat,xmin,xmax,col=2,logx=F,logy=F,genplot=T,verbose=T)

Arguments

See Also

delPts, idPts, trim and trimAT

Description

Tune stratigraphic series to an astronomical target using graphical interface similar to Analyseries 'Linage' routine (Paillard et al, 1996).

Usage

linage(dat,target,extrapolate=F,xmin=NULL,xmax=NULL,tmin=NULL,tmax=NULL,size=1,plotype=1,
       output=1,genplot=T)

Arguments

References

Paillard, D., L. Labeyrie and P. Yiou, 1996), Macintosh program performs time-series analysis: Eos Trans. AGU, v. 77, p. 379.

Examples

## Not run: 
# Check to see if this is an interactive R session, for compliance with CRAN standards.
# YOU CAN SKIP THE FOLLOWING LINE IF YOU ARE USING AN INTERACTIVE SESSION.
if(interactive()) {

# generate example series with 3 precession terms and noise using function 'cycles'
# then convert from time to space using sedimentation rate that increases from 1 to 7 cm/ka
ex=sedRamp(cycles(start=1,end=400, dt=2,noisevar=.00005),srstart=0.01,srend=0.07)

# create astronomical target series
targ=cycles(start=1,end=400,dt=2)

## manually tune
tuned=linage(ex,targ)

## should you need to flip the direction of the astronomical target series, use function 'cb':
tuned=linage(ex,cb(targ[1]*-1,targ[2]))

}

 
## End(Not run)

Description

Interpolate stratigraphic series onto a evenly sampled grid, using piecewise linear interpolation

Usage

linterp(dat,dt,start,genplot=T,check=T,verbose=T)

Arguments

Description

Linearly interpolate stratigraphic series using the approach of Laepple and Huybers (2013), including anti-alias filtering.

Usage

linterpLH13(dat,dt=NULL,start=NULL,dtMin=NULL,thresh=0.7,beta=1,tbw=20,smooth=0.05,
            logF=F, antialias=T,roll=10^10,nsim=1,ncores=2,output=T,maxF=NULL,pl=1,
            genplot=T,verbose=1)

Arguments

Details

This function seeks to estimate the optimal sampling interval for piecewise linear interpolation of an unevenly sampled stratigraphic series ("D"), so as to reduce bias when reconstructing the spectral background. The algorithm is slightly modified from Laepple and Huybers (2013), and consists of the following steps:

(1) Generate a 1/f stochastic noise signal, N. A sampling interval equivalent to the finest step of data set D is used for the noise model if dtMin is NULL.

(2) Subsample the 1/f noise (N) on the actual data set D sampling grid, yielding N2.

(3) Interpolate the subsampled 1/f noise (N2) to dtMin, yielding N3.

(4) Calculate the MTM power spectra for N (Spec_N) and N3 (Spec_N3), using a low resolution spectrum to provide a relatively smooth estimate of the continuum.

(5) Divide the resampled power spectrum (Spec_N3) by unresampled spectrum (Spec_N), yielding the ratio R.

(6) Smooth the ratio R using a Gaussian kernel, yielding R2.

(7) Estimate the highest reliable frequency (F) as the lowest frequency at which the ratio R2 < 0.7 (or specified thresh value).

(8) Calculate the optimal interpolation resolution, dt, as 0.5/F.

If anti-alias filtering is selected:

(9) Interpolate the data set D to 0.1*dt, yielding D2.

(10) Filter D2 using a Taner lowpass (anti-alias) filter with a cutoff frequency of 1.2/dt, yielding D3.

(11) Resample D3 at the optimal resolution dt, yielding D4.

If anti-alias filtering is not selected, instead of steps 9-11, resample D at the optimal resolution dt.

When anti-alias filtering is applied, the first and last value of the interpolated series (D4) have a tendency to be strongly biased if they closely align with the original sampling grid of D. To address this issue, the first and last interpolated points are removed if they are within 0.5*dt (the optimal interpolation interval) of the first or last sampling location of D.

NOTE: A finer (smaller increment) interpolation interval than the 'optimal' one proposed by linterpLH13 may be appropriate for astronomical signal detection, although the magnitude of the observed astronomical variability may be substantially underestimated.

References

T. Laepple and P. Huybers, 2013, Reconciling discrepancies between Uk37 and Mg/Ca reconstructions of Holocene marine temperature variability: Earth and Planetary Science Letters, v. 375, p. 418-429.

Examples

# generate example series from ar1 noise, 5 kyr sampling interval
ex = ar1(npts=501,dt=10)

# jitter sampling times
ex[1]=ex[1]+rnorm(501,sd=1)
# sort
ex = ex[order(ex[,1],na.last=NA,decreasing=FALSE),]

# view sampling statistics
strats(ex)

# run linterpLH13
linterpLH13(ex,output=FALSE)

# run linterpLH13 with multiple simulations
linterpLH13(ex,nsim=500,output=FALSE)

Description

Log transformation of stratigraphic series.

Usage

logT(dat,c=0,opt=1,genplot=T,verbose=T)

Arguments

See Also

arcsinT, demean, detrend, divTrend, prewhiteAR, and prewhiteAR1

Description

Lowpass filter stratigraphic series using rectangular, Gaussian or tapered cosine window. This function can also be used to highpass filter a record (see examples).

Usage

lowpass(dat,padfac=2,fcut=NULL,win=0,demean=T,detrend=F,addmean=T,alpha=3,p=0.25,
        xmin=0,xmax=Nyq,genplot=T,verbose=T)

Arguments

See Also

bandpass, noKernel, noLow, prewhiteAR, prewhiteAR1, and taner

Examples

# generate example series with periods of 405 ka, 100 ka, 40ka, and 20 ka, plus noise
ex=cycles(freqs=c(1/405,1/100,1/40,1/20),end=1000,dt=5,noisevar=.1)

# lowpass filter eccentricity terms using cosine-tapered window
lowpass_ex=lowpass(ex,fcut=.02,win=2,p=.4)

# highpass filter obliquity and precession terms using cosine-tapered window
#  if you'd like the final notch filtered record to be centered on the mean proxy 
#  value, set addmean=FALSE
highpass_ex=lowpass(ex,fcut=.02,win=2,p=.4,addmean=FALSE)
highpass_ex[2] <- ex[2]-highpass_ex[2]
pl(2)
plot(ex,type="l",main="Eccentricity+Obliquity+Precession")
plot(highpass_ex,type="l",main="Obliquity+Precession highpassed signal")

Description

LOWSPEC: Robust Locally-Weighted Regression Spectral Background Estimation (Meyers, 2012)

Usage

lowspec(dat,decimate=NULL,tbw=3,padfac=5,detrend=F,siglevel=0.9,setrho,
         lowspan,b_tun,output=0,CLpwr=T,xmin,xmax,pl=1,sigID=T,genplot=T,
         verbose=T)

Arguments

Details

LOWSPEC is a 'robust' method for spectral background estimation, designed for the identification of potential astronomical signals that are imbedded in red noise (Meyers, 2012). The complete algoritm implemented here is as follows: (1) initial pre-whitening with AR1 filter (default) or other filter as appropriate (e.g., see function prewhiteAR), (2) power spectral estimation via the multitaper method (Thomson, 1982), (3) robust locally weighted estimation of the spectral background using the LOWESS-based (Cleveland, 1979) procedure of Ruckstuhl et al. (2001), (4) assignment of confidence levels using a Chi-square distribution.

NOTE: If you choose to pre-whiten before running LOWSPEC (rather than using the default AR1 pre-whitening), specify setrho=0.

Candidiate astronomical cycles are subsequently identified via isolation of those frequencies that achieve the required (e.g., 90 percent) LOWSPEC confidence level and MTM harmonic F-test confidence level. Allowance is made for the smoothing inherent in the MTM power spectral estimate as compared to the MTM harmonic spectrum. That is, an F-test peak is reported if it achieves the required MTM harmonic confidence level, while also achieving the required LOWSPEC confidence level within +/- half the power spectrum bandwidth resolution. One additional criterion is included to further reduce the false positive rate, a requirement that significant F-tests must occur on a local power spectrum high, which is parameterized as occurring above the local LOWSPEC background estimate. See Meyers (2012) for futher information on the algorithm.

In this implementation, the 'robustness criterion' ('b' in EQ. 6 of Ruckstuhl et al., 2001) has been optimized for 2 and 3 pi DPSS, using a 'span' of 1. By default the robustness criterion will be estimated. Both 'b' and the 'span' can be expliclty set using parameters 'b_tun' and 'lowspan'. Note that it is permissible to decrease 'lowspan' from its default value, but this will result in an overly conservative false positive rate. However, it may be necessary to reduce 'lowspan' to provide an approporiate background fit for some stratigraphic data. Other options include decimating the data series prior to spectral estimation, or using an alternative pre-whitening filter (e.g., see function prewhiteAR).

Note that the 'rfbaseline' function in package 'IDPmisc' (Locher, 2024) is used for implementation of the procedure of Ruckstuhl et al. (2001).

Value

If option 1 is selected, a data frame containing the following is returned: Frequency, Prewhitened power, harmonic F-test CL, LOWSPEC CL, LOWSPEC background, 90%-99% LOWSPEC power levels. NOTE: as of version 0.8, the order of the columns in the output data frame has been changed, for consistency with functions mtm, mtmML96, and mtmPL.

If option 2 is selected, the 'significant' frequencies are returned (as described above).

References

W.S. Cleveland, 1979, Locally weighted regression and smoothing scatterplots: Journal of the American Statistical Association, v. 74, p. 829-836.

Locher R (2024). IDPmisc: 'Utilities of Institute of Data Analyses and Process Design (www.zhaw.ch/idp)'. R package version 1.1.21, https://CRAN.R-project.org/package=IDPmisc.

S.R. Meyers, 2012, Seeing Red in Cyclic Stratigraphy: Spectral Noise Estimation for Astrochronology: Paleoceanography, 27, PA3228, doi:10.1029/2012PA002307.

A.F. Ruckstuhl, M.P Jacobson, R.W. Field, and J.A. Dodd, 2001, Baseline subtraction using robust local regression estimation: Journal of Quantitative Spectroscopy & Radiative Transfer, v. 68, p. 179-193.

D.J. Thomson, 1982, Spectrum estimation and harmonic analysis: IEEE Proceedings, v. 70, p. 1055-1096.

See Also

eha, mtm, mtmAR, mtmML96, periodogram

Examples

# generate example series with periods of 400 ka, 100 ka, 40 ka and 20 ka
ex = cycles(freqs=c(1/400,1/100,1/40,1/20),start=1,end=1000,dt=5)

# add AR1 noise
noise = ar1(npts=200,dt=5,sd=.5)
ex[2] = ex[2] + noise[2]

# LOWSPEC analysis
pl(1, title="lowspec")
lowspec(ex)

# compare to MTM spectral analysis, with conventional AR1 noise test
pl(1,title="mtm")
mtm(ex,ar1=TRUE)

# compare to ML96 analysis
pl(1, title="mtmML96")
mtmML96(ex)

# compare to amplitudes from eha
pl(1,title="eha")
eha(ex,tbw=3,win=1000,pad=1000)

Description

Generate noise surrogates from a theoretical power spectrum.

Usage

makeNoise(S,dt=1,mean=0,sdev=1,addPt=F,nsim=1,genplot=T,verbose=T)

Arguments

Details

These simulations use the random number generator of Matsumoto and Nishimura [1998]. The algorithm of Timmer and Konig (1995) is employed to generate surrogates from any arbitrary theoretical power spectrum. See examples section below.

References

M. Matsumoto, and T. Nishimura, (1998), Mersenne Twister: A 623-dimensionally equidistributed uniform pseudo-random number generator, ACM Transactions on Modeling and Computer Simulation, 8, 3-30.

J. Timmer and K. Konig (1995), On Generating Power Law Noise, Astronomy and Astrophysics: v. 300, p. 707-710.

Examples

# create theoretical AR1 spectrum, using rho of 0.8
rho=0.8
freq=seq(0,.5,by=0.005)
Nyq=max(freq)
AR1 = (1-(rho^2)) / (  1 - (2*rho*cos(pi*freq/Nyq)) + (rho^2) )
plot(freq,AR1,type="l")

# make noise surrogates from the theoretical AR1 spectrum
makeNoise(AR1)

Description

Example stratigraphic model series.

Usage

data(modelA)

Format

Height (meters), weight percent CaCO3

Description

Multitaper method (MTM) spectral analysis (Thomson, 1982)

Usage

mtm(dat,tbw=3,ntap=NULL,padfac=5,demean=T,detrend=F,siglevel=0.9,ar1=T,output=0,
     CLpwr=T,xmin,xmax,pl=1,sigID=T,genplot=T,verbose=T)

Arguments

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.

If ar1=T, candidiate astronomical cycles are identified via isolation of those frequencies that achieve the required (e.g., 90 percent) "red noise" confidence level and MTM harmonic F-test confidence level. Allowance is made for the smoothing inherent in the MTM power spectral estimate as compared to the MTM harmonic spectrum. That is, an F-test peak is reported if it achieves the required MTM harmonic confidence level, while also achieving the required red noise confidence level within +/- half the power spectrum bandwidth resolution. One additional criterion is included to further reduce the false positive rate, a requirement that significant F-tests must occur on a local power spectrum high, which is parameterized as occurring above the local red noise background estimate. See Meyers (2012) for futher information.

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

References

S.R. Meyers, 2012, Seeing Red in Cyclic Stratigraphy: Spectral Noise Estimation for Astrochronology: Paleoceanography, 27, PA3228, doi:10.1029/2012PA002307.

Rahim, K.J. and Burr W.S. and Thomson, D.J., 2014, Applications of Multitaper Spectral Analysis to Nonstationary Data. PhD thesis, Queen's University. R package version 1.0-17, https://CRAN.R-project.org/package=multitaper.

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

See Also

eha, lowspec, mtmAR, mtmML96, periodogram

Examples

# generate example series with periods of 400 ka, 100 ka, 40 ka and 20 ka
ex = cycles(freqs=c(1/400,1/100,1/40,1/20),start=1,end=1000,dt=5)

# add AR1 noise
noise = ar1(npts=200,dt=5,sd=.5)
ex[2] = ex[2] + noise[2]

# MTM spectral analysis, with conventional AR1 noise test
pl(1,title="mtm")
mtm(ex,ar1=TRUE)

# compare to ML96 analysis
pl(1, title="mtmML96")
mtmML96(ex)

# compare to analysis with LOWSPEC
pl(1, title="lowspec")
lowspec(ex)

# compare to amplitudes from eha
pl(1,title="eha")
eha(ex,tbw=3,win=1000,pad=1000)

Description

Perform the 'intermediate spectrum test' of Thomson et al. (2001).

Paraphrased from Thomson et al. (2001): Form an intermediate spectrum by dividing MTM by AR estimate. Choose an order P for a predictor. A variety of formal methods are available in the literature, but practically, one keeps increasing P (the order) until the range of the intermediate spectrum Si(f) (equation (C4) of Thomson et al., 2001) stops decreasing rapidly as a function of P. If the intermediate spectrum is not roughly white, as judged by the minima, the value of P should be increased.

Usage

mtmAR(dat,tbw=3,ntap=NULL,order=1,method="mle",CItype=1,padfac=5,demean=T,detrend=F,
      output=1,xmin=0,xmax=Nyq,pl=1,genplot=T,verbose=T)

Arguments

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

Rahim, K.J. and Burr W.S. and Thomson, D.J., 2014, Applications of Multitaper Spectral Analysis to Nonstationary Data. PhD thesis, Queen's University. R package version 1.0-17, https://CRAN.R-project.org/package=multitaper.

Thomson, D. J., L. J. Lanzerotti, and C. G. Maclennan, 2001, The interplanetary magnetic field: Statistical properties and discrete modes, J. Geophys.Res., 106, 15,941-15,962, doi:10.1029/2000JA000113.

See Also

eha, lowspec, mtm, mtmML96, and periodogram

Examples

# generate example series with periods of 400 ka, 100 ka, 40 ka and 20 ka
ex = cycles(freqs=c(1/400,1/100,1/40,1/20),start=1,end=1000,dt=5)

# add AR1 noise
noise = ar1(npts=200,dt=5,sd=.5)
ex[2] = ex[2] + noise[2]

# MTM spectral analysis, with conventional AR1 noise test
pl(1,title="mtmAR")
mtmAR(ex)