Implement Lomb-Scargle Versions of Fourier algorithms for various products

[dhuppenkothen]

Including the PSD, the Crossspectrum and the Dynamical Powerspectrum. This is also related to #533.

[ranjan2829]

The Lomb-Scargle periodogram is a widely used method for analyzing unevenly sampled time series data. It is particularly useful in astronomy, where observations may not be taken at regular intervals due to the orbital motion of telescopes or other factors.

One possible project for implementing Lomb-Scargle versions of Fourier algorithms would be to adapt existing Fourier-based algorithms in various scientific products to use Lomb-Scargle periodograms instead. This would involve developing new code or modifying existing code to use Lomb-Scargle periodograms as a substitute for traditional Fourier methods.

Some examples of scientific products that could benefit from Lomb-Scargle versions of Fourier algorithms include:

Signal processing tools: Many signal processing tools use Fourier analysis to extract features from signals. By implementing Lomb-Scargle versions of these algorithms, these tools could be extended to handle unevenly sampled data.

Astronomy software: As mentioned earlier, astronomy is a field where unevenly sampled data is common. By implementing Lomb-Scargle versions of Fourier algorithms in astronomy software packages, researchers could more easily analyze astronomical time series data.

Machine learning libraries: Machine learning algorithms often require regular time series data as input. By developing Lomb-Scargle versions of Fourier algorithms for use in machine learning libraries, researchers could more easily analyze unevenly sampled time series data.

Overall, implementing Lomb-Scargle versions of Fourier algorithms in various scientific products could greatly expand the capabilities of these tools and make them more useful for analyzing real-world data.

[ranjan2829]

Here's an example implementation of the Lomb-Scargle periodogram algorithm in Python:

import numpy as np

def lomb_scargle(t, y, freq):
    omega = 2.0 * np.pi * freq
    tau = np.arctan2(np.sum(np.sin(2.0 * omega * t)), np.sum(np.cos(2.0 * omega * t))) / (2.0 * omega)
    tau = tau + np.arange(-1, 2) / (2.0 * freq)
    y_fit = np.zeros_like(y)
    for i in range(len(tau)):
        y_fit += np.interp(t + tau[i], t, y) * np.sinc(freq * (t - tau[i]))
    y_fit /= len(tau)
    pgram = np.sum((y - y_fit) ** 2)
    return pgram

'''

Here, t and y are the time and flux values of the light curve, respectively, and freq is an array of frequency values at which to compute the Lomb-Scargle periodogram. The function returns an array of periodogram values corresponding to each frequency value.

Note that this implementation assumes evenly spaced time values. To handle unevenly spaced time values, you would need to compute the Lomb-Scargle periodogram using a more general algorithm that accounts for the time spacing of the data. There are several Python packages that implement such algorithms, such as AstroML's LombScargle class and the gatspy.periodic.LombScargle class in the gatspy package.

[ranjan2829]

@dhuppenkothen would be happy to contribute for your organization !

[dhuppenkothen]

Hi @ranjan2829 Thank you for your contribution! This is great start! We are indeed mainly interested in the LSP for unevenly sampled data. For this project, we'd like to explore whether we can reuse existing implementations for example in astropy, and we also need to be careful about the normalisation of the LSP. It also needs to follow our existing class structure for our Fourier classes like Powerspectrum, and we have some plans to implement related functionality for example for the cross spectrum. So it's a bit more involved than maybe this issue suggests, but that's why we're planning it as a Google Summer of Code project! :)

For now, if you're interested in contribution, I'd suggest to take a look at the issues labelled good-first-issue.

[ranjan2829]

Hi! It sounds like an exciting project, and I'm glad to hear that my contribution was a good start. I understand that the project is more involved than initially suggested and that you have specific requirements for the implementation, such as normalizing the LSP and following the existing class structure. Reusing existing implementations, such as those in astropy, could be a good starting point, but it may also require some modifications to fit your specific needs.

I'm sure there are several ways to approach this project, and I hope that I can assist you in finding the best solution. Let me know if you have any questions or if there is anything else I can help you with!

[SethiAbhinav]

Hi @dhuppenkothen I have a simple implementation of the LombScargle library of the astropy package on PSD.

I generated a time series and got this result (was just exploring how the library works):

import numpy as np
import matplotlib.pyplot as plt
from astropy.timeseries import LombScargle

# Simulate data
np.random.seed(0)
time = np.sort(np.random.uniform(0, 10, 100))
signal = 5 * np.sin(2 * np.pi * 0.6 * time + np.pi / 4) + np.random.normal(0, 1, 100)

# Compute Lomb-Scargle periodogram
frequency, power = LombScargle(time, signal).autopower(minimum_frequency=0.1, maximum_frequency=2)

# Plot periodogram
plt.plot(frequency, power)
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.show()

After a bit of research I implemented it for PSD like this (just a single line of implementation, however I was mainly learning about what PSD actually is :P):

# Compute PSD
frequency, power = LombScargle(time, signal).autopower(normalization='psd', method='fast')

plt.plot(frequency, power)
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.show()

I wanted some more information on where exactly I should start for this project on your github repo. Also, is there a preferred place (discord/slack etc) for me to continue the discussion?

[dhuppenkothen]

Hi @SethiAbhinav welcome! Our slack is indeed a good place to join the discussion! You can find the link on our documentation website (stingray.science).

For this project, the main objective will be to implement a Lomb-Scargle periodogram in such a way that it matches the normalisations and the structure of existing Stingray functionality. This probably means implementing a new LombScargle class analogous to the existing Powerspectrum and Crossspectrum classes that can take similar inputs and has similar structure and attributes. So looking at the existing functionality in Stingray would be a good option.

[SethiAbhinav]

Got it, thank you!

[SamiraBhattacharya]

Hi @dhuppenkothen. I was writing my proposal and decided to give this code a try. This is just a rough code with fake data. I have written the code according to my proposal. Let me know how I can improve my logic in this code. How does this look?

from stingray import Lightcurve, Crossspectrum, AveragedCrossspectrum, Powerspectrum, AveragedPowerspectrum, LombScargle

# Generate some fake data
t = np.linspace(0, 10, 1000)
y1 = np.sin(2*np.pi*t)
y2 = np.sin(2*np.pi*t + np.pi/4)

# Create two light curve objects
lc1 = Lightcurve(t, y1)
lc2 = Lightcurve(t, y2)

# Compute the Lomb-Scargle periodogram for each light curve
ls1 = LombScargle(lc1.time, lc1.counts)
frequency1, power1 = ls1.compute_power()
ls2 = LombScargle(lc2.time, lc2.counts)
frequency2, power2 = ls2.compute_power()

# Compute the Lomb-Scargle version of the Fourier Cross Spectrum
fcs_ls = Crossspectrum(lc1, lc2, norm='leahy', power_type='lomb-scargle')
freq_fcs_ls, power_fcs_ls = fcs_ls.lombscargle()

# Compute the Lomb-Scargle version of the Fourier Power Spectrum
fps_ls1 = Powerspectrum(lc1, norm='leahy', power_type='lomb-scargle')
freq_fps_ls1, power_fps_ls1 = fps_ls1.lombscargle()
fps_ls2 = Powerspectrum(lc2, norm='leahy', power_type='lomb-scargle')
freq_fps_ls2, power_fps_ls2 = fps_ls2.lombscargle()

# Compute the Lomb-Scargle version of the Averaged Fourier Cross Spectrum (optional)
afcs_ls = AveragedCrossspectrum(lc1, lc2, segment_size=10, norm='leahy', power_type='lomb-scargle')
freq_afcs_ls, power_afcs_ls = afcs_ls.lombscargle()

# Compute the Lomb-Scargle version of the Averaged Fourier Power Spectrum (optional)
afps_ls1 = AveragedPowerspectrum(lc1, segment_size=10, norm='leahy', power_type='lomb-scargle')
freq_afps_ls1, power_afps_ls1 = afps_ls1.lombscargle()
afps_ls2 = AveragedPowerspectrum(lc2, segment_size=10, norm='leahy', power_type='lomb-scargle')
freq_afps_ls2, power_afps_ls2 = afps_ls2.lombscargle()

# Plot the periodograms and the Lomb-Scargle versions of the Fourier products
import matplotlib.pyplot as plt

plt.figure(figsize=(12, 8))

# Plot the Lomb-Scargle periodograms
plt.subplot(231)
plt.plot(frequency1, power1)
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.title('Lightcurve 1 (LS Periodogram)')

plt.subplot(234)
plt.plot(frequency2, power2)
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.title('Lightcurve 2 (LS Periodogram)')

# Plot the Lomb-Scargle version of the Fourier Cross Spectrum
plt.subplot(232)
plt.plot(freq_fcs_ls, power_fcs_ls)
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.title('Lomb-Scargle Cross Spectrum')

# Plot the Lomb-Scargle version```

[matteobachetti]

Also related to: #368

[GitGuardianPro]

please tell me about some good first issue

[GitGuardianPro]

please tell me about some good first issue

[teddu7891]

Hello
here is a very interesting method that I used in my course on signal processing. (Least Squares Method). Very useful for astronomy.

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import inv
import seaborn as sns
import matplotlib.pyplot as plt

def least_square(S, C):
   ls = np.dot(inv(np.dot(C, C.T)) , np.dot(C, S))
   return ls
   
def Gauss(x, mu, sigma):
    gaussienne = 1/(sigma * np.sqrt(2*np.pi)) * np.exp(-0.5*((x-mu)/sigma)**2)
    return(gaussienne) 
    
 # X-axis
x_liste =  np.linspace(90, 210, 1000)

# Creating 3 density fonction

# Density 1
mu_a1 = 120
sigma_a1 = 1.5
mu_a2 = 180
sigma_a2 = 2
gauss_1 =  Gauss(x_liste, mu_a1, sigma_a1) + Gauss(x_liste, mu_a2, sigma_a2)

# Density 2
mu_b1 = 150
sigma_b1 = 15
mu_b2 = 130
sigma_b2 = 20

gauss_2 = Gauss(x_liste, mu_b1, sigma_b1) + Gauss(x_liste, mu_b2, sigma_b2)

# Density 3
mu_c1 = 110
sigma_c1 = 1.9
mu_c2 = 140
sigma_c2 = 10
gauss_3 = Gauss(x_liste, mu_c1, sigma_c1) + Gauss(x_liste, mu_c2, sigma_c2)

# Normalization:
component_1 = gauss_1/np.max(gauss_1)
component_2 = gauss_2/np.max(gauss_2)
component_3 = gauss_3/np.max(gauss_3)

# SHow
sns.set(style="whitegrid")

fig, ax = plt.subplots(figsize=(8, 6))

sns.lineplot(x=x_liste, y=component_1, label='density 1', ax=ax)
sns.lineplot(x=x_liste, y=component_2, label='density 2', ax=ax)
sns.lineplot(x=x_liste, y=component_3, label='density 3', ax=ax)

ax.set_title('3 density fonction', fontsize=15)
ax.set_xlabel('x - axis', fontsize=15)
ax.set_ylabel('Density', fontsize=15)
ax.legend()

plt.show()

# Mix weight
c_1 = 0.3
c_2 = 0.4
c_3 = 0.1

# The mixture
query_spectra = c_1 * component_1 + c_2 * component_2 + c_3 *component_3

# Let's add it some noise for a bit of realism:
query_spectra = query_spectra +  np.random.normal(0, 0.01, len(x_liste))

sns.set(style="whitegrid")

sns.lineplot(x=x_range, y=query_spectra, color='blue', label='Mixture spectrum with noise')

plt.title('Mix', fontsize=12)
plt.xlabel('x-axis', fontsize=12)
plt.ylabel('Total Intensity', fontsize=12)
plt.show()

# Components
components = np.array([component_1, component_2, component_3])

# Apply Least squares
cs = least_sq(query_spectra, components)
# And plot the result:
plt.plot(x_range, query_spectra, color = 'blue', label = 'Mix spectrum' ) # Plot the unknown spectrum
plt.plot(x_range, np.dot(cs,components), color = 'red', linewidth = 2, label = 'Mehtod') # Plot the calculated spectrum

plt.title('Least Squares Method ', fontsize = 12)
plt.xlabel('x-axis', fontsize = 12)
plt.ylabel('total density', fontsize = 12)
plt.legend()
plt.show()

[Oliver-Graz]

Hi, @dhuppenkothen. Our team is trying to make power spectrum out of my data from Swift/XRT, which is an X-ray afterglow lightcurve from a long Gamma-ray Burst. Things went pretty well since we initially used Lomb-Scargle Periodogram (LSP) and Weighted Wavelet Z-transform (WWZ) to make power spectrum and dynamical power spectrum. Then we try to use Stingray to do the same work, but we are stucked this time. Whenever we make power spectrum using Powerspectrum.from_lightcurve as https://docs.stingray.science/index.html has shown, 'AssertionError: No GTIs are equal to or longer than segment_size' will be received.

I have read as much discussion between you and others around this issue as I can find, and I suspect this error occurs
cause my data is unevenly binned. Our team would greatly appreciate if you could provide some help and advice! (The code and data we are using are attached below)

import numpy as np
import matplotlib.pyplot as plt
from stingray import Lightcurve
from stingray import Powerspectrum

data=np.loadtxt("data.txt")
t=data[:,0]#first column in 'data.txt' is time series (units:s)
y=data[:,1]#second column in 'data.txt' is count rate (units:erg/cm^2/s)
#the 3 & 4 column will not be used in this work

gti_1=[[93.66592081,102.21053878],[110.4575,517.4115]]#the XRT is not working between 101.924742938783s & 110.37s
#we take the beginning and the end of gti_1 by -/+ dt/2
lc=Lightcurve(t,y,gti=gti_1)
lc.plot(marker = 'k',labels=('Time', "Flux"))
plt.loglog()
plt.show()

ps = Powerspectrum.from_lightcurve(lc, norm="leahy")#here is where the very error occurs
print(ps)

data.txt

[matteobachetti]

@Oliver-Graz indeed, Stingray has very limited support for unbinned data, and only for the multi-taper periodogram: https://docs.stingray.science/notebooks/Multitaper/multitaper_example.html#Time-series-with-uneven-temporal-sampling:-Multitaper-Lomb-Scargle
We are in the process of implementing the Lomb-Scargle algorithms for the periodogram and the cross spectrum as part of this year's Google Summer of Code

[Oliver-Graz]

@matteobachetti thank you very much for your valuable reply. I have just made a simple test on our data using multi-taper periodogram, and got result similar to what we've got from our own Lomb-Scargle Periodogram code, which makes me feel at ease a little. Just for curious, is there another way to make PSD out of the unevenly-binned data (like the data I've sent in my last comment ) using Stingray (I mean, a Stingray method without using Lomb-Scargle) ? Looking forward to your next reply!

[matteobachetti]

@Oliver-Graz not at the moment, but @pupperemeritus is working on it (#737) and we should have something ready in the next release (by the end of the summer).

[Oliver-Graz]

@matteobachetti , thank you for your help, I will keep on focusing on the project you mentioned, best wishes.

[matteobachetti]

Resolved by #737