---
title: "Lomb Scargle Aliases in R"
format: 
  html:
    embed-resources: true
---
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:
```{r}
#| output: false
library(tidyverse)
library(lomb)
library(pracma)
library(glue)
```

Snagging data:
```{r}
#| output: false
df <- read_csv("./lombscargle_demo.csv")
head(df)
```

Creating initial power spectrum:
```{r}
power <- lsp(df, from = 1 / 100, to = 1 / 0.01, ofac = 2)
```

Grabbing out results to data frame for easier viewing in ggplot:
```{r}
lc1 <- tibble(freq = power$scanned, power = power$power)


ggplot(data = lc1, aes(x = freq, y = power)) +
    geom_line() +
    xlim(0, 2)
```
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`.
```{r}
f1 = power$peak.at[1]
glue("f1 = {f1} hz")
```


Now, to look for aliases, we need to create the spectral window:
```{r}
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.

```{r}
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.
```{r}
f_samp = spec_win$peak.at[1]
glue("f_sample = {f_samp} Hz")
```

And now we can visualize the aliases using our equation.
```{r}
#| code-fold: true
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)
```

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.