Nothing
### This function is a component of astrochron: An R Package for Astrochronology
### Copyright (C) 2021 Stephen R. Meyers
###
###########################################################################
### mtmML96 function - (SRM: December 22, 2013; July 31, 2014; January 31, 2015;
### February 1, 2015; February 5, 2015; February 22, 2015;
### February 26, 2015; March 6, 2015; September 10, 2015;
### December 14-15, 2015; May 20, 2016; August 22, 2016,
### October 4, 2016; March 20, 2017; November 20, 2017;
### August 21, 2018; January 14, 2021)
###
### Perform Mann and Lees (1996) robust red noise mtm analysis, with some modifications.
### Uses multitaper library and built in functions from R.
###########################################################################
mtmML96 <- function (dat,tbw=3,ntap=NULL,padfac=5,demean=T,detrend=F,medsmooth=0.2,opt=1,linLog=2,siglevel=0.9,output=0,CLpwr=T,xmin=0,xmax=Nyq,sigID=T,pl=1,genplot=T,verbose=T)
{
if(verbose) cat("\n----- PERFORMING Mann and Lees (1996) Robust Red Noise Analysis -----\n")
dat <- data.frame(dat)
npts <- length(dat[,1])
dt <- dat[2,1]-dat[1,1]
# error checking
if(dt<0)
{
if (verbose) cat("\n * Sorting data into increasing height/depth/time, removing empty entries\n")
dat <- dat[order(dat[,1], na.last = NA, decreasing = F), ]
dt <- dat[2,1]-dat[1,1]
npts <- length(dat[,1])
}
dtest <- dat[2:npts,1]-dat[1:(npts-1),1]
epsm=1e-9
if( (max(dtest)-min(dtest)) > epsm )
{
cat("\n**** ERROR: sampling interval is not uniform.\n")
stop("**** TERMINATING NOW!")
}
if (verbose)
{
cat(" * Number of data points in stratigraphic series:",npts,"\n")
cat(" * Stratigraphic series length (space or time):",(npts-1)*dt,"\n")
cat(" * Sampling interval (space or time):",dt,"\n")
}
numtap=trunc((2*tbw)-1)
if(is.null(ntap))
{
ntap=numtap
if (verbose) cat(" * Will use default setting of",ntap,"DPSS tapers\n")
}
if(ntap > numtap)
{
ntap=numtap
if (verbose) cat("**** WARNING: The number of DPSS tapers specified is too large. ntap reset to default value of",ntap,"\n")
}
if(ntap < 2)
{
ntap=numtap
if (verbose) cat("**** WARNING: The number of DPSS tapers specified is too small. ntap reset to default value of",ntap,"\n")
}
###########################################################################
### MTM Power spectrum using 'multitaper' library
###########################################################################
# remove mean and linear trend if requested
if (demean)
{
dave <- colMeans(dat[2])
dat[2] <- dat[2] - dave
if(verbose) cat(" * Mean value subtracted=",dave,"\n")
}
if (!demean && verbose) cat(" * Mean value NOT subtracted\n")
# use least-squares fit to remove linear trend
if (detrend)
{
lm.0 <- lm(dat[,2] ~ dat[,1])
dat[2] <- dat[2] - (lm.0$coeff[2]*dat[1] + lm.0$coeff[1])
if(verbose) cat(" * Linear trend subtracted. m=",lm.0$coeff[2],"b=",lm.0$coeff[1],"\n")
}
if (!detrend && verbose) cat(" * Linear trend NOT subtracted\n")
# what is the estimated AR1 coefficient for the combined signal?
lag0 <- dat[1:(npts-1),2]
lag1 <- dat[2:npts,2]
rho_raw <- cor(lag0,lag1)
if(verbose) cat(" * Estimated Conventional AR1 coefficient =",rho_raw,"\n")
# calculate Nyquist freq
Nyq <- 1/(2*dt)
# calculate rayleigh frequency
Ray <- 1/(dt*npts)
npad=as.integer(npts*padfac)
# add another zero if we don't have an even number of data points, so Nyquist exists.
if((npad*padfac)%%2 != 0) npad = npad + 1
# padded frequency grid
df = 1/(npad*dt)
if(verbose)
{
cat(" * Nyquist frequency:",Nyq,"\n")
cat(" * Rayleigh frequency:",Ray,"\n")
cat(" * MTM Power spectrum bandwidth resolution (halfwidth):",tbw/(npts*dt),"\n")
cat(" * Padded to",npad,"points\n")
cat(" * Frequency grid spacing:",df,"\n")
}
# make dat a time series object, here with unit sampling interval
datTS <- as.ts(dat[,2])
spec <- spec.mtm(datTS,nw=tbw,k=ntap,Ftest=T,nFFT=npad,taper=c("dpss"),centre=c("none"),jackknife=F,returnZeroFreq=F,plot=F)
# assign correct frequencies to spec$freq (note: this variable is returned if output = 4)
spec$freq <- spec$freq/dt
# note: no zero frequency present, also remove Nyquist now
nfreq = length(spec$freq) - 1
freq <- spec$freq[1:nfreq]
# normalize power (divided by npts in spec.mtm)
pwrRaw <- spec$spec[1:nfreq]/npts
FtestRaw <- spec$mtm$Ftest[1:nfreq]
###########################################################################
### Robust AR(1) coefficient estimate by analytic fit to median smoothed
### spectrum
###########################################################################
### CALCULATE MEDIAN SMOOTHED SPECTRUM USING runmed FUNCTION #############
if(verbose) cat(" * Calculating median smoothed spectrum using Tukey's robust end-point rule and symmetrical medians\n")
freqSmooth = Nyq * medsmooth
# number of frequencies per each smoothing window
nptsSmooth = as.integer( round(freqSmooth/df, digits=1) )
# check to see if nptsSmooth is even
if(nptsSmooth %% 2 == 0)
{
nptsSmooth=nptsSmooth+1
if(verbose) cat(" * Median smoothing window increased by 1 point\n")
}
if(verbose) cat(" * Number of smoothing points=",nptsSmooth,"\n")
pwrMedian <- runmed(pwrRaw,nptsSmooth,endrule="median")[1:nfreq]
### determine average power (So)
So <- mean(pwrMedian)
###########################################################################
### Analytical fit of AR1 to median smoothed spectrum
###
### Use Brent's method, Gauss-Newton, or brute force grid search approach (if
### other approaches fail).
###########################################################################
# check to ensure logarthim exists
if(min(pwrMedian) < 0)
{
cat("**** ERROR IN ANALYTIC FIT TO MEDIAN SMOOTHER: Spectrum has negative power\n")
stop("**** TERMINATING NOW!")
}
if(opt == 1 || opt == 2)
{
if(opt==1)
{
# define function for Brent's method (option 1)
rednoise1 <- function (rho)
{
rhospec=So * (1-rho^2)/(1 - (2*rho*cos(pi*freq/Nyq)) + rho^2)
if(linLog == 1) sum((pwrMedian-rhospec)^2)
if(linLog == 2) sum((log(pwrMedian)-log(rhospec))^2)
}
if(verbose) cat(" * Calculating analytic fit of AR1 to median smoothed spectrum using Brent's method\n")
MLred=optim(par=rho_raw,rednoise1,method="Brent",lower=0,upper=1)
rho <- MLred$par
}
if(opt==2)
{
# define function for Gauss-Newton algorithm (option 2)
rednoise2 <- function (freq,rho,So,Nyq)
{
So * (1-rho^2)/(1 - (2*rho*cos(pi*freq/Nyq)) + rho^2)
}
if(verbose) cat(" * Calculating analytic fit of AR1 to median smoothed spectrum using Newton-Gauss algorithm\n")
if(verbose) cat(" * Gauss-Newton Convergence for AR1 fit:\n")
if(linLog == 1) MLred <- nls( pwrMedian ~ rednoise2(freq,rho,So,Nyq), data=data.frame(cbind(freq,pwrMedian)), start=list(rho=rho_raw), trace=T )
if(linLog == 2) MLred <- nls( log(pwrMedian) ~ log(rednoise2(freq,rho,So,Nyq)), data=data.frame(cbind(freq,pwrMedian)), start=list(rho=rho_raw), trace=T )
rho <- summary(MLred)$coefficient[1]
}
}
# Grid search
if(opt == 3)
{
if(verbose) cat(" * Calculating analytic fit of AR1 to median smoothed spectrum using brute force grid search\n")
min_MSE = 10^36
for (rho1 in seq(0, 0.999, by=0.001))
{
# calculate mean square error
rhospec = So * (1-rho1^2)/(1 - (2*rho1*cos(pi*freq/Nyq)) + rho1^2)
if(linLog == 1) MSE=sum((pwrMedian-rhospec)^2)
if(linLog == 2) MSE=sum((log(pwrMedian)-log(rhospec))^2)
MSE=MSE/nfreq
if (MSE < min_MSE)
{
min_MSE = MSE
rho = rho1
}
}
}
if(verbose) cat(" * Estimated Robust AR1 coefficient =",rho,"\n")
#***************************************************************
### Calculate Theoretical AR(1) spectrum (EQ's 4 and 5 of Mann and Lees, 1996)
#***************************************************************
### pwrMedAR will be the AR1 power spectrum for the composite signal
pwrMedAR = So * (1-(rho^2)) / ( 1 - (2*rho*cos(pi*freq/Nyq)) + (rho^2) )
# calculate significance levels
dof = (2*ntap)
chiMedAR <- (pwrRaw/pwrMedAR) * dof
chiCLMedAR <- pchisq(chiMedAR, df=dof)
# 90, 95 and 99% CL for power
MedAR1_90 <- pwrMedAR*qchisq(0.9, df=dof)/dof
MedAR1_95 <- pwrMedAR*qchisq(0.95, df=dof)/dof
MedAR1_99 <- pwrMedAR*qchisq(0.99, df=dof)/dof
### f-test CL
dof=2*ntap
prob <- pf(FtestRaw,2,dof-2)
### generate plots
if(genplot)
{
if(pl == 1) logxy="y"
if(pl == 2) logxy=""
if(pl == 3) logxy="xy"
if(pl == 4) logxy="x"
if(pl == 3 || pl == 4) xmin=freq[1]
par(mfrow=c(3,1))
plot(freq,pwrRaw,type="l", col="black", xlim=c(xmin,xmax), xlab="Frequency",ylab="Power",main="MTM Power (black), Robust AR1 fit (red), smoothed (green)",cex.axis=1.1,cex.lab=1.1,lwd=2,bty="n",log=logxy)
lines(freq,pwrMedAR,col="red",lwd=2)
lines(freq,pwrMedian,col="seagreen",lwd=2,lty=3)
if(CLpwr)
{
lines(freq,MedAR1_90,col="red",lwd=1,lty=3)
lines(freq,MedAR1_95,col="red",lwd=1,lty=3)
lines(freq,MedAR1_99,col="red",lwd=1,lty=3)
}
}
###########################################################################
### Identify "significant" frequencies
###########################################################################
if(verbose)
{
cat("\n * Searching for significant spectral peaks that satisfy",siglevel*100,"% CL\n")
cat(" requirements outlined in Meyers (2012):\n")
}
### identify the harmonic f-test peaks
res <- peak(cbind(freq,prob),level=siglevel,genplot=F,verbose=F)
numpeak = length(res[,1])
freqloc = res[,1]
probmax = res[,3]
# FORTRAN wrapper
peakfilter <- function (numpeak,nfreq,tbwRay,siglevel,freqloc,probmax,freq,background,pwr,cl)
{
F_dat = .Fortran('peakfilter_r',
numpeak=as.integer(numpeak),nfreq=as.integer(nfreq),tbwRay=as.double(tbwRay),
siglevel=as.double(siglevel),freqloc=as.integer(freqloc),probmax=as.double(probmax),
freq=as.double(freq),background=as.double(background),pwr=as.double(pwr),
cl=as.double(cl),
loc=integer(as.integer(numpeak)), nout=integer(1)
)
# return the results
return(F_dat)
}
# identify maxima of peaks
tbwRay=tbw*Ray
res2 <- peakfilter(numpeak,nfreq,tbwRay,siglevel,freqloc,probmax,freq,pwrMedAR,pwrRaw,chiCLMedAR)
numpeak2=res2$nout
loc=res2$loc[1:numpeak2]
Frequency <- freq[loc]
Harmonic_CL <- prob[loc]
Red_Noise_CL <- chiCLMedAR[loc]
if(verbose)
{
cat(" * Number of significant F-test peaks identified =",numpeak2,"\n")
cat("ID / Frequency / Period / Harmonic_CL / Rednoise_CL\n")
for(i in 1:numpeak2) cat(i," ", Frequency[i]," ",1/Frequency[i]," ",Harmonic_CL[i]*100," ",Red_Noise_CL[i]*100,"\n")
}
if(genplot)
{
if(pl == 1) logxy=""
if(pl == 2) logxy=""
if(pl == 3) logxy="x"
if(pl == 4) logxy="x"
if(sigID && numpeak2 > 0)
{
### plot "significant" F-test frequencies (on power plot first)
plfreq=double(numpeak2)
pltext=double(numpeak2)
for (k in 1:numpeak2)
{
plfreq[k]=Frequency[k]
pltext[k]=k
}
abline(v=plfreq,col="gray",lty=3)
mtext(pltext[seq(from=1,to=numpeak2,by=2)], side=3,line=0.25,at=plfreq[seq(from=1,to=numpeak2,by=2)],cex=0.5,font=4)
if(numpeak2 > 1) mtext(pltext[seq(from=2,to=numpeak2,by=2)], side=3,line=-0.25,at=plfreq[seq(from=2,to=numpeak2,by=2)],cex=0.5,font=4)
}
plot(freq,chiCLMedAR*100,type="l",col="red",xlim=c(xmin,xmax),ylim=c(0,100),cex.axis=1.1,cex.lab=1.1,lwd=2,xlab="Frequency",ylab="Confidence Level",main="Robust AR1 Confidence Level Estimates",bty="n",log=logxy)
abline(h=c(90,95,99),col="black",lty=3)
if(sigID && numpeak2 > 0)
{
### plot "significant" F-test frequencies (on red noise CL plot)
abline(v=plfreq,col="gray",lty=3)
mtext(pltext[seq(from=1,to=numpeak2,by=2)], side=3,line=0.25,at=plfreq[seq(from=1,to=numpeak2,by=2)],cex=0.5,font=4)
if(numpeak2 > 1) mtext(pltext[seq(from=2,to=numpeak2,by=2)], side=3,line=-0.25,at=plfreq[seq(from=2,to=numpeak2,by=2)],cex=0.5,font=4)
}
plot(freq,prob*100,type="l",col="red",xlim=c(xmin,xmax),ylim=c(80,100),cex.axis=1.1,cex.lab=1.1,xlab="Frequency",ylab="Confidence Level",main="MTM Harmonic F-test Confidence Level Estimates",bty="n",lwd=2,log=logxy)
abline(h=c(90,95,99),col="black",lty=3)
if(sigID && numpeak2 > 0)
{
### plot "significant" F-test frequencies (on probabilty plot)
abline(v=plfreq,col="gray",lty=3)
mtext(pltext[seq(from=1,to=numpeak2,by=2)], side=3,line=0.25,at=plfreq[seq(from=1,to=numpeak2,by=2)],cex=0.5,font=4)
if(numpeak2 > 1) mtext(pltext[seq(from=2,to=numpeak2,by=2)], side=3,line=-0.25,at=plfreq[seq(from=2,to=numpeak2,by=2)],cex=0.5,font=4)
}
}
if (output==1)
{
spectrum <- data.frame(cbind(freq,pwrRaw,prob*100,chiCLMedAR*100,pwrMedAR,MedAR1_90,MedAR1_95,MedAR1_99,pwrMedian))
colnames(spectrum)[1] <- 'Frequency'
colnames(spectrum)[2] <- 'Power'
colnames(spectrum)[3] <- 'Harmonic_CL'
colnames(spectrum)[4] <- 'AR1_CL'
colnames(spectrum)[5] <- 'AR1_fit'
colnames(spectrum)[6] <- 'AR1_90_power'
colnames(spectrum)[7] <- 'AR1_95_power'
colnames(spectrum)[8] <- 'AR1_99_power'
colnames(spectrum)[9] <- 'Median_Smoothed_Power'
return(spectrum)
}
if (output==2)
{
sigfreq <- data.frame(Frequency[1:numpeak2])
colnames(sigfreq) <- 'Frequency'
return(sigfreq)
}
if (output==3)
{
sigfreq <- data.frame(Frequency[1:numpeak2],Harmonic_CL[1:numpeak2])
colnames(sigfreq)[1] <- 'Frequency'
colnames(sigfreq)[2] <- 'Harmonic_CL'
return(sigfreq)
}
#if (output==4)
# {
# return(spec)
# }
#### END function mtmML96
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.