Lomb Scargle Aliases in R

I didn’t have time in class to show the corresponding 1-to-1 R analogs to what we did in Python, so as promised I have tried to put that all together here. With minimal, but some comments to remind you of what was happening along the way.

Grabbing our libraries:

library(tidyverse)
library(lomb)
library(pracma)
library(glue)

Snagging data:

df <- read_csv("./lombscargle_demo.csv")
head(df)

Creating initial power spectrum:

power <- lsp(df, from = 1 / 100, to = 1 / 0.01, ofac = 2)

Grabbing out results to data frame for easier viewing in ggplot:

lc1 <- tibble(freq = power$scanned, power = power$power)


ggplot(data = lc1, aes(x = freq, y = power)) +
    geom_line() +
    xlim(0, 2)

Warning: Removed 19572 rows containing missing values or values outside the scale range
(`geom_line()`).

It looks like the main peak that we’d want is also the largest, so we could just use the peak.at value from power. Or we could use pracma’s findpeaks.

f1 = power$peak.at[1]
glue("f1 = {f1} hz")

f1 = 0.12517775240842 hz

Now, to look for aliases, we need to create the spectral window:

ones <- rep(1, length(df$time))
ones[1] <- 2 # R's lomb scargle breaks on pure 1s
spec_win <- lsp(times = df$time, x = ones, ofac=2)

spec_win_df <- tibble(freq = spec_win$scanned, power = spec_win$power)

Note that we couldn’t use pure one’s for our y-values here. It turns out that R’s Lomb-Scargle approach insists on subtracting out the mean and then dividing to normalize the data, which always results in NA powers. Changing one value to a 2 seems to still get us a weak signal that we can find the peaks of.

ggplot(data=spec_win_df, aes(x=freq, y=power)) + geom_line()

It is pretty clearly the largest peak that we want here, so we can rely on spec_win$peak.at, or we could have used findpeaks.

f_samp = spec_win$peak.at[1]
glue("f_sample = {f_samp} Hz")

f_sample = 0.440625688477638 Hz

And now we can visualize the aliases using our equation.

Code

ggplot(data=lc1, aes(x=freq, y=power)) +
  geom_line() + 
  geom_vline(aes(xintercept=f1), color='red') +
  geom_vline(
    aes(xintercept=x), 
    data=tibble(x=abs(f1 - seq(1,5)*f_samp)),
    color='darkgreen', linetype='dashed'
    ) +
  geom_vline(
    aes(xintercept=x), 
    data=tibble(x=abs(f1 + seq(1,5)*f_samp)),
    color='darkgreen', linetype='dashed'
    ) +
  xlim(0,2)

Warning: Removed 19572 rows containing missing values or values outside the scale range
(`geom_line()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_vline()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_vline()`).

I am noting that my aliases seem to just slightly miss some of the later peaks in a systematic way, which I suspect is from the noisier spectral window that we had to generate here. It should be fine though to establish that things generally line up. If I come up with a more genious way of creating the spectral window in R, I’ll certainly let you know. This issue caught me unawares.