spec.lomb {spectral} | R Documentation |
Lomb-Scargle Periodigram
Description
The Lomb-Scargle periodigram represents a statistical estimator for the amplitude and phase at a given frequency. This function takes also multivariate (n-dimensional) input data.
Usage
spec.lomb(
x = NULL,
y = stop("Missing y-Value"),
f = NULL,
ofac = 1,
w = NULL,
mode = "normal",
maxMem = 8,
cl = NULL
)
Arguments
x |
sampling vector or data frame |
y |
input data vector or data frame |
f |
optional frequency vector / data frame. If not supplied |
ofac |
in case |
w |
weights for data. It must be a 1D vector. |
mode |
|
maxMem |
sets the amount of memory (in MB) to utilize, as a rough approximate. |
cl |
if |
Details
Since the given time series does not need to be evenly sampled, the data
mainly consists of data pairs x1, x2, x3, ...
(sampling points) and (one)
corresponding value y
, which stores the realisation/measurement data.
As can be seen from the data definition above, multivariate (n-dimensional)
input data is allowed and properly processed.
Two different methods are implemented: the standard Lomb-Scargle method with
y(t) = a * cos(\omega (t - \tau)) + b * sin(\omega (t - \tau))
as model function and the generalized Lomb-Scargle (after Zechmeister 2009) method with
y(t) = a * cos(\omega t) + b * sin(\omega t) + c
as model function, which investigates a floating average parameter c
as well.
Both methods can be supplied by an artifical dense frequency vector f
.
In conjunction with the resulting phase information the user might be able to
build a "Fourier"-like spectrum to reconstruct or interpolate the timeseries in equally
spaced sampling. Remind the band limitation which must be fulfilled for this.
- f
The frequencies should be stored in a 1D vector or – in case of multivariate analysis – in a
data.frame
structure to preserve variable namesofac
If the user does not provide a corresponding frequency vector, the
ofac
parameter causes the function to estimatenf = ofac*length(x)/2
equidistant frequencies.
p
-valueThe
p
-value (aka false alarm probability FAP) gives the probability, wheter the estimated amplitude is NOT significant. However, ifp
tends to zero the amplidutde is significant. The user must decide which maximum value is acceptable, until an amplitude is not valid.
If missing values NA
or NaN
appear in any column, the corresponding row
is excluded from calculation.
Value
The spec.lomb
function returns an object of the class lomb
,
which is a list
containg the following information:
A
A vector with amplitude spectrum
f
corresponding frequency vector
phi
phase vector
PSD
power spectral density normalized to the sample variance
floatAvg
floating average value only in case of
mode == "generalized"
w
if,
mode == "generalized"
contains the weighting vectorx,y
original data
- p
p-value False Alarm Probability
Speed Up
In general the function calculates everything in a vectorized manner, which
speeds up the procedure. If the memory requirement is more than maxMem
,
the calculation is split into chunks which fit in the memory (cache). Depending on the
problem size (number of frequencies and data size) a tuning of this value
enhances speed.
Please consider to replpace the BLAS library by a multithreaded version. For example https://prs.ism.ac.jp/~nakama/SurviveGotoBLAS2/binary/windows/x64/ is hosting some Windows RBlas.dll files. Refer to https://mattstats.wordpress.com/2016/02/07/r-with-gotoblas-on-windows-10/ for further information.
The parameter cl
controls a possible cluster, which can be invoked. It takes an
integer number of workers (i. e. cl = 4
), a list with node names c("localhost",...)
or an object of class 'cluster'
or similar. The first two options cause the
function to create the cluster internally. This takes time due to the initialization.
The faster way is to provide an already initialized cluster to the function.
References
A. Mathias, F. Grond, R. Guardans, D. Seese, M. Canela, H. H. Diebner, and G. Baiocchi, "Algorithms for spectral analysis of irregularly sampled time series", Journal of Statistical Software, 11(2), pp. 1–30, 2004.
J. D. Scargle, "Studies in astronomical time series analysis. II - Statistical aspects of spectral analysis of unevenly spaced data", The Astrophysical Journal, 263, pp. 835–853, 1982.
M. Zechmeister and M. Kurster, "The generalised Lomb-Scargle periodogram. A new formalism for the floating-mean and Keplerian periodograms", Astronomy & Astrophysics, 496(2), pp. 577–584, 2009.
See Also
Examples
# create two sin-functions
x_orig <- seq(0,1,by=1e-2)
y_orig <- 2*sin(10*2*pi*x_orig) + 1.5*sin(2*2*pi*x_orig)
# make a 10% gap
i <- round(length(x_orig)*0.2) : round(length(x_orig)*0.3)
x <- x_orig
y <- y_orig
x[i] <- NA
y[i] <- NA
# calculating the lomb periodogram
l <- spec.lomb(x = x, y = y,ofac = 20,mode = "normal")
# select a frequency range
m <- rbind(c(9,11))
# select and reconstruct the most significant component
l2 = filter.lomb(l, x_orig, filt = m)
# plot everything
par(mfrow=c(2,1),mar = c(4,4,2,4))
plot(x,y,"l", main = "Gapped signal")
lines(l2$x, l2$y,lty=2)
legend("bottomleft",c("gapped","10Hz component"),lty=c(1,2))
plot(l,main = "Spectrum")
summary(l)
### Multivariate -- 3D Expample ###
require(lattice)
fx <- 8.1
fy <- 5
fz <- 2
# creating frequency space
f <- expand.grid( fx = seq(-10,10,by = 0.5)
,fy = seq(-10,10,by = 0.5)
,fz = 0:3
)
# creating spatial space
pts <- expand.grid( x = seq(0,1,by = 0.02)
,y = seq(0,1,by = 0.02)
,z = seq(0,1,by = 0.02)
)
# gapping 30%
i <- sample(1:dim(pts)[1],0.7*dim(pts)[1])
pts <- pts[i,]
# caluculating function
pts$val <- cos(2*pi*( fx*pts$x
+ fy*pts$y
+ fz*pts$z
) + pi/4
) +
0.5 * cos(2*pi*( - 0.5 * fx*pts$x
+ 0.5*fy*pts$y
+ 1 * pts$z
) + pi/4
)
# display with lattice
levelplot(val~x+y,pts,subset = z == 0,main = "with z = 0")
# calculating lomb takes a while
# or we sample only a few points
# which enlarges the noise but accelerates the calculation
l <- spec.lomb(y = pts[sample(1:dim(pts)[1],2e3),]
,f = f
# ,mode = "generalized"
)
# name the stripes
l$fz_lev <- factor(x = paste("fz =",l$fz)
)
# display output
levelplot(PSD~fx+fy|fz_lev,l)
# the result is an oversampled spectrum of a non equidistant
# sampled function. We recognize a 3D analysis in all provided
# spatial directions x, y, z.
summary(l)