How to improve the signal-to-noise of a periodogram?

Asteroseismology is the study of stellar oscillations. To see them, we usually want to transfer from the time domain that light curves are in to the frequency domain. We can do that with methods such as Fourier Transforms or Lomb-Scargle Periodograms.

Lightkurve has built in methods for working with data in the frequency domain in the LombScarglePeriodogram class.

Below we demonstrate some of these functionalities. In particular, we will demonstrate how one may optimize the signal-to-noise of a periodogram by varying the detrending, varying the aperture mask, or combining data from multiple observing periods.

import lightkurve as lk
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u

As an example, we can use a red giant star from Campaign 5. You can read more about asteroseismology with giants in K2 in Stello et al 2015 here.

First we’ll use search_targetpixelfile to download the Target Pixel File for the target.

TARGET = "EPIC 211416749"
tpf = lk.search_targetpixelfile(TARGET, author="K2", campaign=5, cadence='short').download()

Be aware that this is a short cadence light curve, which means that this dataset is much larger than others we’ve used in previous tutorials. If you’re rerunning this tutorial locally, it might take a few minutes to rerun the later steps.

Let’s plot the target below to see if the aperture is centered. We’re plotting in ‘log’ scale here so that we can see the wings of the PSF.

tpf.plot(scale='log', aperture_mask=tpf.pipeline_mask);

We can now create a KeplerLightCurve using Simple Aperture Photometry with the to_lightcurve method. We can then normalize, remove NaN values and remove outliers in one easy line. We can also use the new fill_gaps method to fill in any gaps that are in our data using linear interpolation of the nearest neighbours. This creates an almost perfectly sampled time array with no gaps.

lc = tpf.to_lightcurve()
lc = lc.normalize(unit='ppm').remove_nans().remove_outliers().fill_gaps()

We can plot our light curve easily with the plot method.

lc.plot();

Let’s take a look at the periodogram for this light curve. You can create a periodogram object using lc.periodogram(). This can be passed any array of frequencies that the user wants, but by default will create an array of frequency for you. You can also change the frequency units, depending on the range of frequencies you’re looking for. In asteroseismology we usually use \(\mu Hz\).

pg = lc.to_periodogram(freq_unit=u.microHertz, maximum_frequency=400, minimum_frequency=10)

With this new periodogram object we can now plot the power spectrum of the data. Let’s do that now!

ax = pg.plot();

1. Varying the SFF Motion Correction

Unfortunately there is some ringing going on…there is a periodicity in the data due to the K2 motion. We can see below when we plot the data in “Period” space instead of “Frequency” space there is a significant periodicity at ~ 6 hours.

ax = pg.plot(view='Period', unit=u.hour)
ax.text(6, 2000, '6 Hour Roll Motion');

We can apply our own Self Flat Fielding (SFF) correction to the light curve to remove the K2 roll motion. Below we correct with SFF using the default correct() method, with windows=10. (You can read more about SFF in our other tutorials.)

lc = tpf.to_lightcurve().normalize().remove_nans().remove_outliers()
clc = lc.to_corrector("sff").correct(windows=10).remove_outliers().fill_gaps()

clc.plot();

Our new, corrected lightcurve looks much flatter and has a much smaller CDPP. Let’s try plotting the periodogram again.

pg_clean = clc.to_periodogram(freq_unit=u.microHertz, maximum_frequency=400, minimum_frequency=100)
ax = pg_clean.plot()

Now we can see the oscillation modes of the target! We can even find out what effect our detrending parameters have on our target. Below we change the detrending windows from \(w=2\; ... \; 10\) and plot the periodogram in each case. We can see different detrending parameters alter the strength of each mode.

# Loop over several windows
from tqdm import tqdm
for windows in tqdm([2, 4, 8]):
    # Create the light curve
    lc =  tpf.to_lightcurve().normalize().remove_nans().remove_outliers()
    clc = lc.to_corrector('sff').correct(windows=windows).remove_outliers().fill_gaps()
    # Create the periodogram
    p_clean = clc.to_periodogram(freq_unit=u.microHertz, maximum_frequency=400, minimum_frequency=100)
    # Plot the periodogram
    if windows == 2:
        ax = p_clean.plot(alpha=0.4, label='Windows: {}'.format(windows))
    else:
        p_clean.plot(ax=ax, alpha=0.4, label='Windows: {}'.format(windows))
ax.legend();

100%|█████████████████████████████████████████████| 3/3 [00:11<00:00,  3.83s/it]

It looks like there is a significant effect on the oscillation modes as we vary the window size. When using K2 data for detecting these oscillations it is important to understand and mitigate the effects of detrending to ensure the modes have the highest possible signal to noise ratio.

2. Varying the aperture mask size

We can also easily change the aperture size. If the aperture is too small, increasing the aperture size should allow us to include more of the target flux. If there is a contaminant nearby, decreasing the aperture may increase our singal to noise.

Below we create two alternative apertures for the target. In this case we’ll create contiguous apertures where pixels have a value greater than 50 and 2\(\;\sigma\) above the median, respectively.

aperture_small = tpf.create_threshold_mask(threshold=50)
aperture_large = tpf.create_threshold_mask(threshold=2)

We can plot up our new aperture against the pipeline aperture for comparison.

# Two subplots
fig, axs = plt.subplots(1, 2, figsize=(12, 4))

# Plot pipeline aperture mask
tpf.plot(axs[0], scale='log', aperture_mask=aperture_small)
axs[0].set_title('Small aperture')

# Plot larger aperture mask
tpf.plot(axs[1], scale='log', aperture_mask=aperture_large)
axs[1].set_title('Large aperture');

Now let’s create our two light curves, one with the pipeline aperture and one with the new aperture.

tpf.to_lightcurve() creates our SAP flux with the pipeline aperture as a default. To use a new aperture we simple use tpf.to_lightcurve(aperture_mask=NEW_APERTURE).

# Create the light curve with the pipelien aperture.
lc_small = tpf.to_lightcurve(aperture_mask=aperture_small).normalize().remove_nans().remove_outliers()
lc_small = lc_small.to_corrector('sff').correct().remove_outliers().fill_gaps()

# Create a light curve with a slightly larger aperture
lc_large = tpf.to_lightcurve(aperture_mask=aperture_large).normalize().remove_nans().remove_outliers()
lc_large = lc_large.to_corrector('sff').correct().remove_outliers().fill_gaps()

When we plot these two light curves we can see that the larger aperture is much less noisy.

#Plot the pipeline and large aperture light curves
ax = lc_small.plot(label='Small aperture')
lc_large.plot(ax=ax, label='Large aperture');

Finally, when we plot the periodogram we can see we’ve increased the signal to noise ratio of our stellar oscillations.

# Create the periodograms
pg_small = lc_small.to_periodogram(freq_unit=u.microHertz, maximum_frequency=400, minimum_frequency=100)
pg_large = lc_large.to_periodogram(freq_unit=u.microHertz, maximum_frequency=400, minimum_frequency=100)


# Plot the periodograms
ax = pg_small.plot(c='k', alpha=0.5, label='Small aperture')
pg_large.plot(ax=ax, c='C3', alpha=0.5, label='Large aperture')

ax.legend();

By increasing our aperture size and including more pixels here we have increased the signal to noise of the oscillation modes of this red giant.

3. Using data from multiple campaigns

Finally, this target was actually observed twice by K2, once in April 2015 in Campaign 5 and once in December 2017 in Campaign 16. We can get both of these data sets and compare the results.

# Download the C16 TPF
tpf_c16 = lk.search_targetpixelfile(TARGET, author="K2", campaign=16, cadence='short').download()

# Create a new light curve for C16
aperture_mask_c16 = tpf_c16.create_threshold_mask(threshold=2)
lc_c16 = tpf_c16.to_lightcurve(aperture_mask=aperture_mask_c16).normalize().remove_nans().remove_outliers()
lc_c16 = lc_c16.to_corrector('sff').correct(windows=10).remove_outliers().fill_gaps()

# Create a periodogram for c16 data
pg_c16 = lc_c16.to_periodogram(freq_unit=u.microHertz, maximum_frequency=400, minimum_frequency=100)

# Create subplots to plot into
ax = pg_large.plot(c='k', alpha=0.5, label='C5')
pg_c16.plot(ax=ax, c='C3', alpha=0.5, label='C16')

ax.set_xlim(100, 400)
ax.legend();

It looks like the two data sets provide similar modes, however the two campaigns have very different instrument systematics. To find the true answer, users should iterate over many different detrending parameters and aperture sizes, and combine datasets to increase signal to noise.

previous

Verifying the location of a periodic signal in the pixel data

next

How to estimate a star’s mass and radius using asteroseismology