pwrLawFit: Estimate power law (1/f) fit to power spectrum

In astrochron: A Computational Tool for Astrochronology

Estimate power law (1/f) fit to power spectrum

Description

Estimate power law (1/f) fit to power spectrum, following the algorithm of Vaughan (2005).

Usage

pwrLawFit(spec,dof=2,flow=NULL,fhigh=NULL,output=1,genplot=T,verbose=T)

Arguments

spec	Power spectrum. First column is frequency, second column is raw power (linear). Do not include the zero frequency and Nyquist.
dof	Degrees of freedom for power spectral estimate. Default is 2, for a simple periodogram.
flow	Lowest frequency to include in 1/f fit
fhigh	Highest frequency to include in 1/f fit
output	Output results of 1/f fit? (0=none; 1=Frequency,Power,Power Law CL,Unbiased Power Law fit,CL_90,CL_95,CL_99; 2=beta, unbiased log10N, biased log10N)
genplot	generate summary plots (T or F)
verbose	verbose output (T or F)

References

Vaughan, S. (2005), A simple test for periodic signals in red noise, Astronomy & Astrophysics.

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]

# calculate periodogram
res=periodogram(ex,output=1,padfac=1)

# extract power and remove the Nyquist frequency
resPwr=cb(res,c(1,3))
resPwr=resPwr[-length(resPwr[,1]),]

pwrLawFit(resPwr)

astrochron documentation built on Sept. 30, 2024, 9:14 a.m.