Nothing
### This function is a component of astrochron: An R Package for Astrochronology
### Copyright (C) 2017 Stephen R. Meyers
###
###########################################################################
### integratePower: determine the variance within a given bandwidth, and
### ratio it to the total power in the spectrum
### (SRM: December 13, 2013; December 14, 2013; December 17, 2013;
### December 18, 2013; December 24, 2013; January 15-16, 2015;
### March 6, 2015; March 19, 2015; February 14, 2017;
### August 30, 2021)
###
###########################################################################
integratePower <- function (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)
{
if(verbose) cat("\n----- INTEGRATING POWER SPECTRUM -----\n")
if(unity && verbose) cat("\n * NOTE: Spectrum will be normalized such that total variance (up to fmax) is unity.\n")
if(!unity && is.null(npts))
{
cat("\n**** ERROR: Please specify the number of points (npts) in the processed time series window.\n")
stop("**** TERMINATING NOW!")
}
if(!unity && is.null(pad))
{
cat("\n**** ERROR: Please specify the total padding length used (pad).\n")
stop("**** TERMINATING NOW!")
}
if(pad < npts)
{
cat("\n**** ERROR: pad must be larger than, or equal to, npts.\n")
stop("**** TERMINATING NOW!")
}
if( palette != 1 && palette != 2 && palette != 3 && palette != 4 && palette != 5 && palette != 6)
{
cat("\n**** WARNING: palette option not valid. Will use palette = 6.\n")
palette = 6
}
if(!unity) fac=0.5*(pad/npts)
if(unity) fac=0.5
# ensure we have a data frame
spec=data.frame(spec)
# assign frequencies from first column of spec
freq=spec[,1]
cols=length(spec)
if(cols > 2)
{
# assign locations for each spectrum (column headers)
loc=suppressWarnings(as.numeric(substr(colnames(spec[2:cols]),start=2,stop=100)))
# for negative depth/height/time values, "-" has been changed to "."
# this will create NAs. implement modification of fix recommended by Mathieu Martinez
neg=grepl(".",substr(colnames(spec[2:cols]), start=2,stop=2),fixed=T)
fixloc=which(neg)
if(any(neg)) {loc[fixloc]=-1*as.numeric(substr(colnames(spec[(fixloc+1)]),start=3,stop=100))}
# assign specta (amplitude, power, or probability)
sp=as.matrix( spec[2:cols] )
numrec=length(loc)
if(is.null(ymin)) ymin = min(loc)
if(is.null(ymax)) ymax = max(loc)
}
if(cols == 2)
{
numrec=1
sp=spec[,2]
}
spPlot=sp
if(ln) spPlot=log(spPlot)
numfreq=length(freq)
if(verbose) cat("\n * Number of spectra to analyze =",numrec,"\n")
if(verbose) cat(" * Number of frequencies per spectrum =",numfreq,"\n")
if(is.null(xmin)) xmin = min(freq)
if(is.null(xmax)) xmax = max(freq)
if(is.null(fmax)) fmax = max(freq)
if(genplot) plotit=T
if(!genplot)
{
if(is.null(flow) || is.null(fhigh))
{
plotit=T
cat("\n **** NOTE: flow and/or fhigh have not been specified. Interactive plotting will be activated.\n")
}
if(!is.null(flow) && !is.null(fhigh)) plotit=F
}
if(plotit && numrec>1)
{
# set color palette
# rainbow colors
if(palette == 1) colPalette = tim.colors(ncolors)
# grayscale
if(palette == 2) colPalette = gray.colors(n=ncolors,start=1,end=0,gamma=1.75)
# dark blue scale (from larry.colors)
if(palette == 3) colPalette = colorRampPalette(c("white","royalblue"))(ncolors)
# red scale
if(palette == 4) colPalette = colorRampPalette(c("white","red2"))(ncolors)
# blue to red plot
if(palette == 5) colPalette = append(colorRampPalette(c("royalblue","white"))(ncolors/2),colorRampPalette(c("white","red2"))(ncolors/2))
# viridis colormap
if(palette == 6) colPalette = viridis(ncolors, alpha = 1, begin = 0, end = 1, direction = 1, option = "D")
# spectrogram plot
dev.new(height=h,width=w)
mat <- matrix(c(1, 2, 3, 4), nrow = 1, ncol = 4)
layout(mat, widths = c(3.5, 1.5, 1.5, 1.5))
# use fields library for access to 'tim.colors'
# useRaster=T results in a faster plotting time.
xlimset=c(xmin,xmax)
if (ydir == -1)
{
# in this case, reset ylim range.
# note that useRaster=T is not a viable option, as it will plot the results backwards, even though the
# y-axis scale has been reversed! This option will result in a slower plotting time.
ylimset=c(ymax,ymin)
image(freq,loc,spPlot,xlim=xlimset,ylim=ylimset,col = tim.colors(ncolors),xlab="Frequency",ylab="Location",main="")
if(is.null(flow) || is.null(fhigh)) mtext("Select frequency band by clicking twice on plot",side=3,line=1)
}
if (ydir == 1)
{
# useRaster=T results in a faster plotting time.
ylimset=c(ymin,ymax)
image(freq,loc,spPlot,xlim=xlimset,ylim=ylimset,col = colPalette,useRaster=T,xlab="Frequency",ylab="Location",main="")
if(is.null(flow) || is.null(fhigh)) mtext("Select frequency band by clicking twice on plot",side=3,line=1)
}
if(is.null(flow) || is.null(fhigh))
# interactive plot selection
{
cat("\n PLEASE CLICK ON PLOT TO SELECT TWO FREQUENCIES TO DEFINE POWER INTEGRATION BAND\n")
# Now overlay x-y plot for graphical interface
par(new=T)
plot(-1,-1,xlim=xlimset,ylim=ylimset,xaxs="i",yaxs="i",yaxt='n',bty='n',ylab="",xlab="")
ff = locator(n = 2, type = "p" ,col="white",lwd=2,pch=3)
flow=min(ff$x)
fhigh=max(ff$x)
}
abline(v=c(flow,fhigh),col="white",lty=3,lwd=2)
}
if(plotit && numrec==1)
{
# spectrogram plot
dev.new(height=5,width=7)
par(mfrow=c(1,1))
plot(freq,spPlot,type="l",xlim=c(xmin,xmax),xlab="Frequency",ylab="",bty = "n", lwd = 2, cex.axis = 1.1, cex.lab = 1.1)
if(is.null(flow) || is.null(fhigh))
# interactive plot selection
{
mtext("Select frequency band by clicking twice on plot",side=3,line=1)
cat("\n PLEASE CLICK ON PLOT TO SELECT TWO FREQUENCIES TO DEFINE POWER INTEGRATION BAND\n")
ff = locator(n = 2, type = "p" ,col="red",lwd=1,pch=3)
flow=min(ff$x)
fhigh=max(ff$x)
}
abline(v=c(flow,fhigh),col="red",lty=3,lwd=2)
}
# isolate spectrum up to fmax
ifreqMax= which( (freq <= fmax) )
# isolate the frequency band that you want to integrate
ifreq= which( (freq >= flow) & (freq <= fhigh) )
# only count f0 once (use 0.5 since fac will double others to get power from neg freqs)
if(f0) sp[1]=sp[1]*0.5
# perform integration
bandPwr=double(numrec)
totalPwr=double(numrec)
for (i in 1:numrec)
{
if(numrec>1)
{
if(verbose) cat(" * Integrating window=",i,"; Location=",loc[i],"\n")
sp1=sp[ifreqMax,i]
sp2 = sp[ifreq,i]
}
if(numrec==1)
{
sp1=sp[ifreqMax]
sp2 = sp[ifreq]
}
# integrate up to fmax
totalPwr[i]=sum(sp1)/fac
# integrate band of interest
bandPwr[i]=sum(sp2)/fac
if(unity) bandPwr[i]=bandPwr[i]/totalPwr[i]
if(unity) totalPwr[i] = 1
}
if(verbose)
{
cat("\n Frequency Band selected for Power Integration=",flow,",",fhigh,"\n")
cat("\n Equivalent periods=",1/flow,",",1/fhigh,"\n\n")
}
if(numrec > 1)
{
out <- data.frame(cbind(loc,bandPwr,totalPwr,bandPwr/totalPwr))
colnames(out)<-c("Location","Band Power","Total Power","Band/Total Power")
if(plotit)
{
plot(bandPwr,loc,type="l",ylim=ylimset,main="Band Power",ylab="Location",xlab="Cumulative Power")
plot(totalPwr,loc,type="l",ylim=ylimset,main="Total Power",ylab="Location",xlab="Cumulative Power")
plot(bandPwr/totalPwr,loc,type="l",ylim=ylimset,main="Band/Total Power",ylab="Location",xlab="Fraction of Power")
}
}
if(numrec==1)
{
out <- data.frame(cbind(bandPwr,totalPwr,bandPwr/totalPwr))
colnames(out)<-c("Band Power","Total Power","Band/Total Power")
}
return(out)
### END function integratePower
}
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.