# Lomb-Scargle Periodograms for Multiband Data#

The Lomb-Scargle periodogram (after Lomb [1], and Scargle [2]) is a commonly
used statistical tool designed to detect periodic signals in unevenly spaced
observations. The base `LombScargle`

provides an
interface for several implementations of the Lomb-Scargle periodogram. However,
`LombScargle`

only handles a single band of data.
The `LombScargleMultiband`

class adapts this
interface to handle multiband data (where multiple bands/filters are present).

The code here is adapted from the astroml package ([3], [4]) and the
gatspy package ([5], [6]), but conforms closely to the design paradigms
established in `LombScargle`

. For a detailed
practical discussion of the Multiband Lomb-Scargle periodogram, which guided
the development of this class, see
*Periodograms for Multiband Astronomical Time Series* [6].

## Basic Usage#

Note

As in `LombScargle`

, frequencies in
`LombScargleMultiband`

are **not**
angular frequencies, but rather frequencies of oscillation (i.e., number of
cycles per unit time).

The Lomb-Scargle Multiband periodogram is designed to detect periodic signals in unevenly spaced observations with multiple bands of data present.

### Example#

To detect periodic signals in unevenly spaced observations, consider the following multiband data, where 5 bands (u, g, r, i, and z) have 60 datapoints each.

```
>>> import numpy as np
>>> t = []
>>> y = []
>>> bands = []
>>> dy = []
>>> N=60
>>> for i, band in enumerate(['u','g','r','i','z']):
... rng = np.random.default_rng(i)
... t_band = 300 * rng.random(N)
... y_band = 3 + 2 * np.sin(2 * np.pi * t_band)
... dy_band = 0.01 * (0.5 + rng.random(N)) # uncertainties
... y_band += dy_band * rng.standard_normal(N)
... t += list(t_band)
... y += list(y_band)
... dy += list(dy_band)
... bands += [band] * N
```

The Lomb-Scargle periodogram, evaluated at frequencies chosen
automatically based on the input data, can be computed as follows
using the `LombScargleMultiband`

class,
with the `bands`

argument being the sole difference in comparison
to the `LombScargle`

interface:

```
>>> from astropy.timeseries import LombScargleMultiband
>>> frequency,power = LombScargleMultiband(t, y, bands, dy).autopower()
```

Plotting the result with Matplotlib gives:

The periodogram shows a clear spike at a frequency of 1 cycle per unit time,
as we would expect from the data we constructed. The resulting power is a
single array, with combined input from each of the bands dependent upon the
implementation chosen in the `method`

keyword.

## Periodograms from `TimeSeries`

objects#

`LombScargleMultiband`

is able to operate on
`TimeSeries`

objects, provided the
`TimeSeries`

object meets a formatting requirement.
The requirement is that the flux (or magnitudes) and errors for each band are
provided in separate columns. If instead, your
`TimeSeries`

object has a singular flux column with
an associated band label column, these columns may be passed directly to
`LombScargleMultiband`

as 1-d arrays.

### Example#

Consider the following generator code for a
`TimeSeries`

object where timeseries data is
populated for three photometric bands (g,r,i).

```
>>> from astropy.timeseries import LombScargleMultiband, TimeSeries
>>> from astropy.table import MaskedColumn
>>> import numpy as np
>>> import astropy.units as u
```

```
>>> rng = np.random.default_rng(1)
>>> deltas = 240 * rng.random(180)
>>> ts1 = TimeSeries(time_start="2011-01-01T00:00:00",
... time_delta=deltas*u.minute)
```

```
>>> # g band fluxes
>>> g_flux = [0] * 180 * u.mJy
>>> g_err = [0] * 180 * u.mJy
>>> y_g = np.round(3 + 2 * np.sin(10 * np.pi * ts1['time'].mjd[0:60]),3)
>>> dy_g = np.round(0.01 * (0.5 + rng.random(60)), 3) # uncertainties
>>> g_flux.value[0:60] = y_g
>>> g_err.value[0:60] = dy_g
>>> ts1["g_flux"] = MaskedColumn(g_flux, mask=[False]*60+[True]*120)
>>> ts1["g_err"] = MaskedColumn(g_err, mask=[False]*60+[True]*120)
>>> # r band fluxes
>>> r_flux = [0] * 180 * u.mJy
>>> r_err = [0] * 180 * u.mJy
>>> y_r = np.round(3 + 2 * np.sin(10 * np.pi * ts1['time'].mjd[60:120]),3)
>>> dy_r = np.round(0.01 * (0.5 + rng.random(60)), 3) # uncertainties
>>> r_flux.value[60:120] = y_r
>>> r_err.value[60:120] = dy_r
>>> ts1['r_flux'] = MaskedColumn(r_flux, mask=[True]*60+[False]*60+[True]*60)
>>> ts1['r_err'] = MaskedColumn(r_err, mask=[True]*60+[False]*60+[True]*60)
>>> # i band fluxes
>>> i_flux = [0] * 180 * u.mJy
>>> i_err = [0] * 180 * u.mJy
>>> y_i = np.round(3 + 2 * np.sin(10 * np.pi * ts1['time'].mjd[120:]),3)
>>> dy_i = np.round(0.01 * (0.5 + rng.random(60)), 3) # uncertainties
>>> i_flux.value[120:] = y_i
>>> i_err.value[120:] = dy_i
>>> ts1["i_flux"] = MaskedColumn(i_flux, mask=[True]*120+[False]*60)
>>> ts1["i_err"] = MaskedColumn(i_err, mask=[True]*120+[False]*60)
>>> ts1
<TimeSeries length=180>
time g_flux g_err r_flux r_err i_flux i_err
mJy mJy mJy mJy mJy mJy
Time float64 float64 float64 float64 float64 float64
----------------------- ------- ------- ------- ------- ------- -------
2011-01-01T00:00:00.000 3.0 0.012 ——— ——— ——— ———
2011-01-01T02:02:50.231 3.891 0.009 ——— ——— ——— ———
2011-01-01T05:50:56.909 4.961 0.007 ——— ——— ——— ———
2011-01-01T06:25:32.807 4.697 0.014 ——— ——— ——— ———
2011-01-01T10:13:13.359 4.451 0.005 ——— ——— ——— ———
2011-01-01T11:28:03.732 4.283 0.008 ——— ——— ——— ———
2011-01-01T13:09:39.633 1.003 0.015 ——— ——— ——— ———
2011-01-01T16:28:18.550 3.833 0.008 ——— ——— ——— ———
2011-01-01T18:06:31.018 1.02 0.013 ——— ——— ——— ———
... ... ... ... ... ... ...
2011-01-15T13:01:17.603 ——— ——— ——— ——— 1.054 0.008
2011-01-15T16:03:17.207 ——— ——— ——— ——— 4.656 0.014
2011-01-15T17:29:38.139 ——— ——— ——— ——— 1.423 0.01
2011-01-15T20:03:35.935 ——— ——— ——— ——— 4.805 0.008
2011-01-15T21:35:02.069 ——— ——— ——— ——— 3.042 0.007
2011-01-15T23:06:35.567 ——— ——— ——— ——— 1.162 0.01
2011-01-16T01:07:30.330 ——— ——— ——— ——— 4.99 0.009
2011-01-16T01:11:31.138 ——— ——— ——— ——— 5.0 0.011
2011-01-16T03:09:58.569 ——— ——— ——— ——— 1.314 0.01
2011-01-16T07:03:09.586 ——— ——— ——— ——— 3.383 0.005
```

Our timeseries data is set up to be asynchronous, where a given timestamp
corresponds to a measurement in a single band. However, if your data instead
has one timestamp per multiple band measurements, or a mixture,
`LombScargleMultiband`

will still be able to
operate on it.

To operate on the example `TimeSeries`

,
`LombScargleMultiband`

has a loader function, as
follows:

```
>>> ls = LombScargleMultiband.from_timeseries(ts1, signal_column=['g_flux', 'r_flux', 'i_flux'],
... uncertainty_column=['g_err', 'r_err', 'i_err'],
... band_labels=['g', 'r', 'i'])
```

`signal_column`

requires a list of columns that correspond to the flux
or magnitude measurements in each band. `uncertainty_column`

and
`band_labels`

are optional, but if specified must be lists of equal size to
`signal_column`

. `uncertainty_column`

specifies the columns containing the
associated errors per band, while `band_labels`

provides the labels to use
for each photometric band. From here,
`LombScargleMultiband`

can be worked with as normal.
For example:

```
>>> frequency,power = ls.autopower()
```

## Consistencies with `LombScargle`

#

`LombScargleMultiband`

is an inherited class of
`LombScargle`

, and was developed to provide as
similar of an interface to `LombScargle`

as
possible. From this, there are several core aspects of
`LombScargle`

that remain true for
`LombScargleMultiband`

.

### Measurement Uncertainties#

The `LombScargleMultiband`

interface can also handle
data with measurement uncertainties. As shown in the example above.

### Periodograms and Units#

The `LombScargleMultiband`

interface properly
handles `Quantity`

objects with units attached,
and will validate the inputs to make sure units are appropriate.

### Specifying the Frequency Grid#

As shown above, the `autopower()`

method automatically determines a frequency grid, using
`autofrequency()`

. The tunable parameters
are identical to those shown for `LombScargle`

. And
likewise, a custom frequency grid may be supplied directly to the
`power()`

function.

#### Example#

```
>>> frequency = np.linspace(0, 2, 1000)
>>> power = LombScargleMultiband(t, y, bands, dy).power(frequency)
```

### Periodogram Implementations#

Two implementations of the Multiband Lomb-Scargle Periodogram are available
within `LombScargleMultiband`

, `flexible`

and
`fast`

, which are selectable via the
`power()`

method’s
`method`

parameter. `flexible`

is a direct port of the LombScargleMultiband
algorithm used in the gatspy gatspy package. It constructs a common model,
and an offset model per individual band. It then applies regularization to the
resulting model to constrain complexity, resulting in a flexible model for any
given multiband timeseries dataset. As it’s name implies, `fast`

is
potentially quicker alternative that fits each band independently and combines
them by weight. The independent band-by-band fits leverage
`LombScargle`

. As a result the `sb_method`

parameter is available in
`power()`

to choose the
single-band method used in `power()`

for
each band. Keep in mind that the speed of `fast`

is dependent on the
underlying speed of the choice of `sb_method`

.

#### Example#

`flexible`

:

```
>>> frequency, power = LombScargleMultiband(t,y,bands,dy).autopower(method='flexible')
```

`fast`

, with `fast`

also chosen as the
`power()`

method:

```
>>> frequency, power = LombScargleMultiband(t,y,bands,dy).autopower(method='fast', sb_method='fast')
```

### The Multiband Lomb-Scargle Model#

The `model()`

method fits a
sinusoidal model to the data at a chosen frequency. The sinusoidal model
complexity is tunable via the `nterms_base`

and `nterms_band`

parameters.
These control the number of sinusoidal terms available to the base model
(common to all bands) and the number of sinusoidal terms available to each
bands offset model.

Note

Either of `nterms_base`

and `nterms_band`

may be set to 0, though not
both. The case when `nterms_base`

=0 and `nterms_band`

=1 is a special
case referred to as the *multi-phase model*, where the base model is reduced
to a simple offset, and therefore the bands are solved independently (a
single-band fit). Further discussed in
*Periodograms for Multiband Astronomical Time Series* [6]

#### Example#

The following example uses the same data as above.
`autopower()`

is used to return
the periodogram, and we can select the frequency at which the power is maximum
for our model:

```
>>> model = LombScargleMultiband(t, y, bands, dy, nterms_base=1, nterms_band=1)
>>> frequency, power = model.autopower(method='flexible')
>>> freq_maxpower = frequency[np.argmax(power)]
```

We can then model based on the found frequency, and time (phased by the frequency):

```
>>> t_phase = np.linspace(0, 1/freq_maxpower, 100)
>>> y_fit = model.model(t_phase, freq_maxpower)
```

The resulting fit is then of shape (number of bands, number of timesteps), or (5,100) in this particular case. By plotting the result, we see the model has recovered the expected sinusoid recovered at the correct frequency:

### False Alarm Probabilities#

Unlike `LombScargle`

,
`LombScargleMultiband`

does not have False Alarm
Probabilities implemented. The algorithms available for
`LombScargle`

are valid only for single term
periodograms, which is rarely valid for models in the Multiband case.