lsp: Lomb Scargle Periodogram

In Tape09/LSP: LSP: Lomb Scargle Periodogram

Description

Calculates the Lomb Scargle Periodogram using the algorithm described by K. Hocke and N. Kämpfer (2009).

Usage

lsp(t, y, ofac = 8, mc_sim = NULL, zero_factor = 0, thresh = 1,
  verbose = 0)

Arguments

t	Numeric vector of timepoints.
y	Numeric vector of values corresponding to t.
ofac	Oversampling factor. Recommend >=2. Higher improves granularity, but may cause artefacts.
mc_sim	monte carlo simulation object returned by function "monte_carlo_lsp". If not NULL will return p-values. If NULL will use spd for thresholding; using the "thresh" parameter
zero_factor	Integer, zero padding factor. Pads fourier series to increase resolution of approximated function. Number of zeros is equal to zero_factor*length of (unpadded) fourier series. Value of 0 means no padding.
thresh	Threshhold for p-value or Spectral Power Density. Filter fourier series based on peak strength p value if mc_sim not NULL, otherise on spd. This will affect the reconstructed series.
verbose	Verbosity level. 0=silent, 1=printed results.

Value

A list with the following elements:
t - original t.
y - original y.
t2 - approximated t.
y2 - reconstructed y. plot against t2.
frequency - frequency range used.
spectral_power_density - Plot against Frequency for the periodogram. unitless.
phase_radian - phase relative to 0 in radians. range 0:2*pi.
phase_default - phase relative to 0 in the same unit as t. This is the distance from 0 to the nearest peak after t[1].
peak_info - data.frame with index, frequency, spd, phase_default, and (optional) pvalue of each detected peak.
peak_idx - index of highest peak
peak_spd - highest spd
peak_period - period corresponding to peak_idx
peak_phase - phase_default corresponding to peak_idx

References

Hocke, K., and N. Kämpfer. "Gap filling and noise reduction of unevenly sampled data by means of the Lomb-Scargle periodogram." Atmospheric Chemistry and Physics 9.12 (2009): 4197-4206.

Examples

library(LSP)

#set up testing parameters
ofac = 8;
amp=3;
sd = 1.0;
n = 100; #initial series length
n_out = 50; #sereis length after random sampling
period = 45;
freq = 1/period;
offset = 25;
phase = pi/2;
p_thresh = 0.05;
N_mc_runs = 1000;
zero_factor = 2;
N_cores = 1;

#create testing data
x = seq(0,n-1);
x = x + offset;
y = amp*cos(2*pi*freq*x - phase);

#random sampling to simulate missing data/irregular timepoints
idx = seq(1,n);
idx = sample(idx,n_out);
idx = sort(idx)

#add some noise
y = rnorm(n,0,sd) + y;

#final data
y = y[idx];
x = x[idx];

#monte carlo
mc_sim = monte_carlo_lsp(x,N_mc_runs,ofac,N_cores);

#LSP
# result = lsp(x,y,ofac=ofac,mc_sim = mc_sim,zero_factor = zero_factor, thresh = p_thresh,verbose=1);
result = lsp(x,y,ofac=ofac,mc_sim = mc_sim,zero_factor = zero_factor, thresh = 0.05,verbose=1);

#plot input data vs reconstructed signal
plot(result$t,result$y,col="red")
lines(result$t2,result$y2,col="green")

Tape09/LSP documentation built on May 9, 2019, 4:19 p.m.