In this post, we will reconstruct a time-varying signal that we can't measure at consistent intervals. A lot of the tools in time-domain signal analysis are derived for continuous signals, with their discrete (digital) forms dependent on consistent sampling. Engineers talk about signaling in terms of frequencies and sample rates, assuming that in the ideal case that you'll have one measurement every \(\Delta t\) seconds. (If a measurement is missing, that means your equipment isn't working right.) The well known discrete Fourier transform (and the FFT algorithm) are based on the idea that you're sampling your signal at regular intervals, and can give you unexpected results when you violate that assumption.

On the other hand, when you're looking for signals from outer space, you might set out to measure once a night. But then you have a cloudy night, and miss a data point. And then your allotted telescope time runs out, so you have a gap before you can get more data.

Or a space telescope has to fire its thrusters mid-observation, throwing off a few measurements in your time series. (Which then might inspire you to write a blog post.)

# The Lomb-Scargle periodogram

The discrete Fourier transform has the undesirable characteristic of introducing
noise at longer periods when the data contains gaps. This can be mitigated
slightly by manipulating the data set to "fill in the gaps", but ultimately the
Fourier transform is the wrong tool for the job. The technique I'm going to use
here comes under the broader category of "least-squares spectral
analysis", and
is called the **Lomb-Scargle periodogram**.

As anyone who has fit functions to data knows, the least-squares regression
technique is a powerful way to determine the appropriate values for your free
parameters. If you had a single sine wave signal, measured at a bunch of
different times, you could use a least-squares fit to recover its amplitude and
frequency. The Lomb-Scargle technique effectively does this for a spectrum of
frequencies of interest, telling you how much power in the signal exists at each
frequency. Fortunately, we don't have to implement the specifics ourselves: an
optimized implementation exists in *SciPy*.

Let's set up our notebook.

```
%pylab inline --no-import-all
import scipy.signal
```

The implementation we'll be using lives in
`scipy.signal.lombscargle`

. The documentation
provides references to the literature, if you want to read up. The
implementation itself is only a few lines of Cython.

# Our input signal

At least to start with, we should use a signal we synthesize ourselves. Then we can be sure we're getting the same thing out as we're putting in, giving us an important sanity check on our code.

Our input signal will have a period of three days, and a phase offset \(\phi\) of \(\pi/2\). Below, I define these parameters for our input signal.

```
# Conversion factors
days_to_seconds = 60 * 60 * 24
seconds_to_days = 1.0 / days_to_seconds
# Signal parameters
A = 2. # amplitude
period = 3 * days_to_seconds # seconds
frequency = 1.0 / period # Hertz
omega = 2. * np.pi * frequency # radians per second
phi = 0.5 * np.pi # radians
```

Since this is a discretely sampled signal, we also have to define the length of our signal in terms of a number of data points and the spacing between them. (This may seem at odds with the idea that we have gaps in our data. You're right! The gaps are introduced a little later on.)

We make an array of `timesteps`

and a corresponding array `signal`

defined as
\(\text{signal}_i = A * \sin ( \omega t_i + \phi )\) for each element *i* in
`timesteps`

and `signal`

.

```
N = 1000 # number of samples we're dealing with
dt = 30 * 60 / 1.0 # 30 minutes each sample, ~5e-4 samples / sec
timesteps = np.linspace(0.0, N*dt, N)
signal = A * np.sin(omega * timesteps + phi)
```

```
# make some prettier time-axis labels, since timesteps is in seconds
timestep_day_ticks = days_to_seconds * np.arange(22)
timestep_day_labels = ['{}'.format(i) for i in np.arange(22)]
# plot the signal
plt.figure(figsize=(14,4))
plt.plot(timesteps, signal)
plt.xticks(timestep_day_ticks, timestep_day_labels)
plt.xlabel('Time (days)')
plt.ylabel('Signal Value')
```

# A moth-eaten sine wave

To get data that demonstrates the strengths of the Lomb-Scargle periodogram, we
have to take our regularly-sampled signal array and throw out points from it. To
do so, we use `numpy.random.rand`

to generate an array the same length as our
input with uniform random values, and throw out a set fraction of the input
based on the random array element and a threshold.

Here we set `frac_points`

to 0.05, meaning we want to retain 5% of the input
samples for our next step.

```
frac_points = 0.05 # fraction of points to select
r = np.random.rand(signal.shape[0])
timesteps_subset = timesteps[r >= (1.0 - frac_points)]
signal_subset = signal[r >= (1.0 - frac_points)]
```

```
plt.figure(figsize=(14,4))
plt.plot(timesteps, signal, 'k', alpha=0.3, label="Original signal")
plt.plot(timesteps_subset, signal_subset, 'b+', label="Points kept in new signal")
plt.xticks(timestep_day_ticks, timestep_day_labels)
plt.xlabel('Time (days)')
plt.ylabel('Signal Value')
plt.legend()
```

This perhaps makes it a bit too easy to see the underlying sine wave. If you
were taking these data, and didn't know *a priori* what you were looking for,
you'd probably have a plot that looked like this instead.

```
plt.figure(figsize=(14,4))
plt.plot(timesteps_subset, signal_subset, 'b+')
plt.xticks(timestep_day_ticks, timestep_day_labels)
plt.xlabel('Time (days)')
plt.ylabel('Signal Value')
plt.ylim(-4, 4)
```

Perhaps you can guess that there's some periodicity there by looking, but it's
much less obvious. This data set is also unnaturally clean by construction. The
*only* signal present above is a sinusoid; real data have systematic errors that
will further confound your analysis.

# Computing the periodogram

We call `scipy.signal.lombscargle`

with three arguments: `x`

, `y`

, and `freqs`

.
In this case, `x`

and `y`

are given by the timestamps and measurements from our
gap-filled signal array. We generate `freqs`

, the array of angular frequencies
over which the periodogram should be evaluated, by picking 10000 periods in the
range from 0.1 to 10 days and using the relation \(\text{Angular Frequency} =
\frac{2 \pi}{\text{Period}}\).

```
# Choosing points in frequency space depends on what frequencies we think we'll detect.
# There's no point looking for 100 Hz signals when we're looking at T = 3 days periods.
nout = 10000 # number of frequency-space points at which to calculate the signal strength (output)
periods = np.linspace(0.1 * days_to_seconds, 10 * days_to_seconds, nout)
freqs = 1.0 / periods
angular_freqs = 2 * np.pi * freqs
```

```
pgram = scipy.signal.lombscargle(timesteps_subset, signal_subset, angular_freqs)
```

From the documentation, we have one caveat about the output given by
`scipy.signal.lombscargle`

:

The computed periodogram is unnormalized; it takes the value

`(A**2) * N/4`

for a harmonic signal with amplitude A for sufficiently large N.

Therefore, if we want to get out the same amplitude \(A\) that we put in, we will have to divide by the number of samples, multiply by four, and take the square- root. In practice, one generally cares more about the frequency than the detected amplitude, so this is not always strictly necessary.

```
normalized_pgram = np.sqrt(4 * (pgram / signal_subset.shape[0]))
```

Let's take a look at the normalized periodogram. Below, I'm plotting the
normalized periodogram output against the *periods* it was evaluated at, instead
of the *angular frequencies* that were passed to the `lombscargle`

function.

For reference, I've also added a vertical line marking the 3 day period of the input signal.

```
plt.figure(figsize=(14,4))
plt.plot(periods, normalized_pgram)
plt.xlabel('Period $T$ (days)')
plt.axvline(3 * days_to_seconds, lw=2, color='red', alpha=0.4)
day_ticks, day_labels = np.arange(10) * days_to_seconds, ['{:2.1f}'.format(t) for t in np.arange(10)]
plt.xticks(day_ticks, day_labels)
plt.tight_layout()
```

# Visualizing the detected period

Once you have a periodogram that hints at a periodicity, you will most likely want to "fold" or "phase" the curve (pick your favorite term) to get a better idea of how each observation fits into your attempted fit.

We take a copy of the `timesteps_subset`

array and use the modulus operator to
limit the values to the range [0, 3] days.

```
timesteps_phased = timesteps_subset.copy()
timesteps_phased = timesteps_phased % (3 * days_to_seconds)
phase_angles = 2 * np.pi * timesteps_phased / (3 * days_to_seconds)
```

Now we plot the original and phased signals. The plot is color coded so that you can see how the points are rearranged on the (former) time axis to assume the appropriate phase position.

```
fig = plt.figure(figsize=(14,4))
# Gap-filled data
ax1 = plt.subplot2grid((1,3), (0,0), colspan=2)
ax1.scatter(timesteps_subset, signal_subset, c=timesteps_subset, cmap='rainbow')
ax1.set_xticks(timestep_day_ticks)
ax1.set_xticklabels(timestep_day_labels)
ax1.set_xlabel('Time (days)')
ax1.set_ylabel('Signal Value')
ax1.set_xlim(0, 21 * days_to_seconds)
# Phased data
ax2 = plt.subplot2grid((1,3), (0,2), colspan=1)
ax2.scatter(phase_angles, signal_subset, c=timesteps_subset, cmap='rainbow')
ax2.set_xticks([0, 0.5 * np.pi, np.pi, 1.5 * np.pi, 2 * np.pi], ['0', '$\pi/2$', '$\pi$', '$3\pi/2$', '$2\pi$'])
ax2.set_xlim(0, 2*np.pi)
ax2.set_xlabel('Phase $\phi$ (radians)')
ax2.set_ylabel('Signal Value')
```

How pretty! Much easier to interpret than that sparsely populated graph of a moth-eaten sine wave.

Hopefully this helped explain how to accomplish something useful to your research. Some more high level descriptions of Lomb-Scargle spectral analysis, as well as other periodograms and when to use them, is online as part of the NASA Exoplanet Archive Periodogram Service documentation. (Bear in mind that they have their own implementation of Lomb-Scargle which may follow different conventions than the SciPy one.)

If you desire more mathematical rigor, in addition to the original paper by Scargle (and the older one by Lomb), there's a treatment (with some useful generalizations) by G. Larry Bretthorst linked below. Thanks to Nick Earl from Space Telescope for sharing it with me.

## Further Reading and References

`scipy.signal.lombscargle`

(SciPy documentation)- Least Squares Spectral Analysis (Wikipedia)
- "Frequency Estimation and Generalized Lomb-Scargle Periodograms" (G. Larry Bretthorst, from Statistical Challenges in Astronomy)