2
$\begingroup$

I am trying to use the periodogram to tell when a signal is periodic or not by following the tutorial for the astropy Lomb-scargle periodogram here.

http://docs.astropy.org/en/stable/stats/lombscargle.html

I simulated some data, one that is a sinusoid (period = 200) and one which is a skew gaussian (i.e, a single transient event). The hope was that periodogram would pick out the period for the periodic object and and give a period for the transient that would imply a single event occured in the window.

Unfortunately, the results don't make sense at all. I have attached code at the end and the figures generated below. Due to random noise in my simulation, each result is different and I provide two examples of the results below. I use the same method outlined in the link above where we use the function LombScargle.model() on the best fit frequency from frequency[argmax[power]]

The red line is the true function I simulated the data from. The green is the best fit from the periodogram. The right hand plots are the PSD from the periodogram.

Example 1 enter image description here

Here we can see the best fit frequency for the sinusoid (top right plot) picked out is 0.105 (i.e. a period of 9-10 days) which is not near a frequency of 0.05 I'd expect for something with a period of 200 days, yet when I feed this best fit frequency of 0.105 to the lomb-scargle model fitter, a nicely matching periodic curve comes out with a period of 200 days?

This does not make sense.

Example 2

enter image description here

Here I ran the code again and this time the results are switched around? It fit a very large period to the transient so that I can confidently say it is a single transient event yet the sinusoidal fit is terrible. The best frequency is still 0.105 (period=10) yet the lomb-scargle model fitter overlays something that appears to have a period of 60 days which is wrong?

Could I get clarification on if i'm doing something wrong? I've been told the periodogram is the de facto tool for unenvenly sampled data like this yet... the results seem awful half the time.

To clarify my questions

  1. Can it be explained how in the first plot the best-fit frequncy of 0.105 that I feed into the astropy lomb-scargle model fitter somehow creates a correctly matching sinusoid with frequency 0.05? What is the explanation?

  2. Why are there 5 strong peaks in the top right plot for the periodogram in example 1 when I only expect 1? The middle two are close to the real value of 0.05 (at 0.045 and 0.055)

Here is the short code used to simulate and plot the data

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sc
from astropy.stats import LombScargle
import math


#simulate parameters
data_range = [i for i in range(1,1001)]
number_of_samples = 50
gauss_skew = sc.skewnorm.pdf
skew = -10
period = 200
location = data_range[int(len(data_range)/2)]


y1= [(2 * (1. + np.sin(2. * np. pi * x/period)) + np.random.normal(loc =0.0, scale = 0.5)) for x in data_range]
errors1 = [np.random.normal(loc = 0.0, scale = 1) for x in data_range]
y2 = [(1000* (gauss_skew(x,skew,loc=location ,scale = 50)) + np.random.normal(loc =0.0, scale = 1)) for x in data_range]
errors2 = [np.random.normal(loc = 0.0, scale = 1) for x in data_range]

sample_rate = int(len(data_range)/number_of_samples)# To thin the data a bit
y1 = y1[::sample_rate]
y2 = y2[::sample_rate]
errors1 = errors1[::sample_rate]
errors2 = errors2[::sample_rate]
x1= x2 =data_range[::sample_rate]

truth1 = [2* (1. + np.sin(2. * np. pi * x/period)) for x in data_range]
truth2 = [1000 * (gauss_skew(x,-10,loc=location ,scale = 50)) for x in data_range]

truths = [truth1,truth2]
x= [x1,x2]
y=[y1,y2]
errors = [errors1,errors2]

fig,ax  = plt.subplots(nrows=2,ncols=2)

ax[0][0].errorbar(x1,y1,yerr=errors1,fmt='o')
ax[0][0].set_xlabel('time')
ax[1][0].errorbar(x2,y2,yerr=errors2,fmt='o')
ax[1][0].set_xlabel('time')
ax[0][1].set_xlabel('frequency')
ax[1][1].set_xlabel('frequency')


for i in range(0,2):
    frequency, power = LombScargle(x[i],y[i],errors[i]).autopower()
    #Get the best fit frequency as in the tutorial
    best_frequency = frequency[np.argmax(power)]
    print('best frequency:',best_frequency)
    t_fit = np.linspace(x[i][0], math.floor(x[i][-1]),num =50)

    #Fit the best fit frequency
    #plot the best best model based on the best fit
    y_fit = LombScargle(x[i], y[i], errors[i]).model(t_fit, best_frequency)
    ax[i][0].plot(t_fit,y_fit,'g')
    ax[i][0].plot(data_range,truths[i],'r')

    #Plot the PSD
    ax[i][1].plot(frequency,power)
    ax[i][1].axvline(x=best_frequency,color='black', ls='--')
plt.show()
$\endgroup$
3
  • $\begingroup$ Why are you making your error bar have a random size and applying a different error to the actual data? I noted on the plot that some error bars were tiny compared to the others. The data are weighted by the square of the error, so this could have unpredictable consequences. $\endgroup$
    – ProfRob
    Mar 12, 2019 at 7:17
  • $\begingroup$ Hi, I am trying to simulate telescope observational data which does not have uniform percentage errors for each data point for observational reasons. At least, the feature of tiny error bars on some data points and large ones on others is something i've seen before in astronomical data, such as the data from the Kaggle LSST Plasticc challenge. This is just an initial approximate simulation and I understand it won't be perfect, but by eye, periodicity is clear to me and I thought a well-known tool like the Lomb-Scargle would cope fine with this. $\endgroup$
    – wrahman
    Mar 12, 2019 at 9:38
  • $\begingroup$ I'm trying to determine here, if the tool is just failing here or if there is something wrong with my code or that I have a fundamental misunderstanding of the how the tool works. Any help appreciated, thanks! $\endgroup$
    – wrahman
    Mar 12, 2019 at 9:41

1 Answer 1

2
$\begingroup$

I think I see what is going on. Your best fit is actually a beat frequency between the true frequency (0.005) and double the sampling frequency (0.05). This does indeed produce a model that goes through your data points but because you have only plotted 50 points in the model you haven't been able to see it.

If you change this line t_fit = np.linspace(x[i][0], math.floor(x[i][-1]),num =1000)

so that it shows the model at more points you can see that it is indeed at a much higher frequency.

The trick is to know (roughly) the answer before you start.

If you limit your frequency range in the Lomb-Scargle command to below the Nyquist frequency:

frequency, power = LombScargle(x[i],y[i],errors[i]).autopower(nyquist_factor=1.5)

You will find the behaviour you expected. e.g. (NB I reduced the error bar size to focus in on your problem, but also fixed the error bar size to a single positive value. Your code was introducing very small and even negative errors and I do not know how the code would deal with those).

Periodograms

In unevenly spaced data, things are not so simple because the Nyquist frequency is not well defined.

$\endgroup$
2
  • $\begingroup$ Yes, thank you! This is pretty much it. I have been reading more about this in a brilliant periodogram review below arxiv.org/pdf/1703.09824.pdf and it detailed the effects of the sampling frequency being convolved with the function frequency which I hadn't considered. $\endgroup$
    – wrahman
    Mar 12, 2019 at 13:00
  • 1
    $\begingroup$ Hey Rob, do you understand what the significance of Lomb-Scargle power is? It's different from the power-spectrum obtained with FFT, and I'm trying to understand what the definition is. I've been asking this question multiple times but got no helpful answer. Please check here $\endgroup$
    – Eric Kim
    Jul 16, 2019 at 16:12

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .