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.
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
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.
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.
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
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?
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.errorbar(x1,y1,yerr=errors1,fmt='o') ax.set_xlabel('time') ax.errorbar(x2,y2,yerr=errors2,fmt='o') ax.set_xlabel('time') ax.set_xlabel('frequency') ax.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], 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].plot(t_fit,y_fit,'g') ax[i].plot(data_range,truths[i],'r') #Plot the PSD ax[i].plot(frequency,power) ax[i].axvline(x=best_frequency,color='black', ls='--') plt.show()