Download notebook Run in Google Colab

Periodogram¶

The Periodogram feature extractor finds peaks in the Lomb–Scargle periodogram of an unevenly sampled time series, returning period and signal-to-noise ratio for each peak:

$$\mathrm{S/N}(\omega_\mathrm{peak}) = \frac{P(\omega_\mathrm{peak}) - \langle P(\omega) \rangle}{\sigma_{P(\omega)}}$$

This tutorial covers parameter selection, reading the full power spectrum, multi-period detection, and using spectrum/phase features.

In [1]:

# %pip install light-curve matplotlib

# %pip install light-curve matplotlib

In [2]:

import light_curve as licu
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
PERIOD = 23.5  # days
t = np.sort(rng.uniform(0, 300, 150))
m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)
print(f'n_obs: {len(t)}, time span: {t[-1] - t[0]:.0f} d, true period: {PERIOD} d')

import light_curve as licu import numpy as np import matplotlib.pyplot as plt rng = np.random.default_rng(0) PERIOD = 23.5 # days t = np.sort(rng.uniform(0, 300, 150)) m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150) print(f'n_obs: {len(t)}, time span: {t[-1] - t[0]:.0f} d, true period: {PERIOD} d')

n_obs: 150, time span: 298 d, true period: 23.5 d

Basic usage¶

Periodogram(peaks=N) returns the N strongest periods and their S/N ratios. With peaks=1 you get two numbers: [period, snr].

In [3]:

import numpy as np

rng = np.random.default_rng(0)
PERIOD = 23.5
t = np.sort(rng.uniform(0, 300, 150))
m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)

pg = licu.Periodogram(peaks=1, fast=True)
result = pg(t, m)
for name, val in zip(pg.names, result):
    print(f'  {name}: {val:.4f}')

import numpy as np rng = np.random.default_rng(0) PERIOD = 23.5 t = np.sort(rng.uniform(0, 300, 150)) m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150) pg = licu.Periodogram(peaks=1, fast=True) result = pg(t, m) for name, val in zip(pg.names, result): print(f' {name}: {val:.4f}')

  periodogram_period_0: 23.6491
  periodogram_period_s_to_n_0: 12.5051

Key parameters¶

Parameter	Default	What it controls
peaks	1	Number of period/S/N pairs returned
resolution	10	Frequency grid oversampling factor
max_freq_factor	1.0	Upper limit multiplier on the Nyquist frequency
nyquist	'average'	How the Nyquist frequency is estimated ('average', 'median', or a float quantile)
fast	True	NFFT-based O(N log N) algorithm; False uses exact O(N²)
normalization	'psd'	Power normalisation (see below)
freqs	None	Explicit angular frequency grid; overrides resolution, max_freq_factor, nyquist

resolution controls how densely the frequency grid is sampled relative to the minimum resolvable frequency (1/T_baseline). Higher values find weaker signals but increase runtime linearly. Values of 10–20 work well for most survey data.

normalization options:

Value	Range	Notes
'psd' (default)	[0, ∞)	Raw power; good for S/N comparisons
'standard'	[0, 1]	$P_\text{std} = P \cdot 2/(n-1)$; matches astropy
'model'	[0, ∞)	$P_\text{std}/(1-P_\text{std})$; matches astropy
'log'	[0, ∞)	$-\ln(1-P_\text{std})$; matches astropy

freqs: pass an explicit angular frequency array (radians per time unit) to control the grid exactly. For fast=True, the grid must be np.linspace(0.0, max_freq, 2**k + 1) for integer k.

Power spectrum¶

Periodogram.freq_power(t, m) returns the full (angular_frequencies, power) arrays, so you can visualise or post-process the complete spectrum without any external library. Frequencies are in radians per time unit; convert to period with $T = 2\pi/\omega$.

In [4]:

import light_curve as licu
import numpy as np

rng = np.random.default_rng(0)
PERIOD = 23.5
t = np.sort(rng.uniform(0, 300, 150))
m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)

pg = licu.Periodogram(peaks=1, resolution=10)
omega, power = pg.freq_power(t, m)  # angular freqs (rad/day)

# Skip omega=0 to avoid division by zero when converting to period
mask = omega > 0
periods = 2 * np.pi / omega[mask]
pwr = power[mask]

best_period = periods[np.argmax(pwr)]

fig, ax = plt.subplots(figsize=(9, 3))
ax.plot(periods, pwr, lw=0.7, color='steelblue')
ax.axvline(best_period, color='red', lw=1.5, ls='--', label=f'Best = {best_period:.2f} d')
ax.axvline(PERIOD, color='green', lw=1.5, ls=':', label=f'True = {PERIOD} d')
ax.set_xscale('log')
ax.set_xlabel('Period (days)')
ax.set_ylabel('LS power (psd)')
ax.set_title('Lomb–Scargle power spectrum')
ax.legend()
plt.tight_layout()
plt.show()
print(f'Recovered period: {best_period:.3f} d  (true: {PERIOD} d)')

import light_curve as licu import numpy as np rng = np.random.default_rng(0) PERIOD = 23.5 t = np.sort(rng.uniform(0, 300, 150)) m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150) pg = licu.Periodogram(peaks=1, resolution=10) omega, power = pg.freq_power(t, m) # angular freqs (rad/day) # Skip omega=0 to avoid division by zero when converting to period mask = omega > 0 periods = 2 * np.pi / omega[mask] pwr = power[mask] best_period = periods[np.argmax(pwr)] fig, ax = plt.subplots(figsize=(9, 3)) ax.plot(periods, pwr, lw=0.7, color='steelblue') ax.axvline(best_period, color='red', lw=1.5, ls='--', label=f'Best = {best_period:.2f} d') ax.axvline(PERIOD, color='green', lw=1.5, ls=':', label=f'True = {PERIOD} d') ax.set_xscale('log') ax.set_xlabel('Period (days)') ax.set_ylabel('LS power (psd)') ax.set_title('Lomb–Scargle power spectrum') ax.legend() plt.tight_layout() plt.show() print(f'Recovered period: {best_period:.3f} d (true: {PERIOD} d)')

Recovered period: 23.649 d  (true: 23.5 d)

Multiple peaks and multi-period detection¶

When the light curve contains several periodic signals, peaks=N returns the N strongest ones. S/N drops sharply for noise peaks, making it a natural threshold. Here we synthesise a light curve with two real periods and ask for five peaks:

In [5]:

import light_curve as licu
import numpy as np

rng = np.random.default_rng(0)
P1, P2 = 23.5, 7.3  # two embedded periods
t = np.sort(rng.uniform(0, 400, 200))
m = (15.0
     + 0.40 * np.sin(2 * np.pi * t / P1)
     + 0.25 * np.sin(2 * np.pi * t / P2)
     + rng.normal(0, 0.04, 200))

pg5 = licu.Periodogram(peaks=5, resolution=20)
result5 = pg5(t, m)

print(f'True periods: {P1} d and {P2} d')
print('Top 5 peaks:')
for i in range(5):
    period = result5[i * 2]
    snr = result5[i * 2 + 1]
    marker = ' ← real' if abs(period - P1) < 1 or abs(period - P2) < 1 else ''
    print(f'  peak {i + 1}: {period:6.2f} d, S/N = {snr:5.2f}{marker}')

import light_curve as licu import numpy as np rng = np.random.default_rng(0) P1, P2 = 23.5, 7.3 # two embedded periods t = np.sort(rng.uniform(0, 400, 200)) m = (15.0 + 0.40 * np.sin(2 * np.pi * t / P1) + 0.25 * np.sin(2 * np.pi * t / P2) + rng.normal(0, 0.04, 200)) pg5 = licu.Periodogram(peaks=5, resolution=20) result5 = pg5(t, m) print(f'True periods: {P1} d and {P2} d') print('Top 5 peaks:') for i in range(5): period = result5[i * 2] snr = result5[i * 2 + 1] marker = ' ← real' if abs(period - P1) < 1 or abs(period - P2) < 1 else '' print(f' peak {i + 1}: {period:6.2f} d, S/N = {snr:5.2f}{marker}')

True periods: 23.5 d and 7.3 d
Top 5 peaks:
  peak 1:  23.52 d, S/N = 12.12 ← real
  peak 2:   7.30 d, S/N =  3.16 ← real
  peak 3:  31.98 d, S/N =  0.68
  peak 4:  25.63 d, S/N =  0.62
  peak 5:  29.29 d, S/N =  0.50

Spectrum features¶

features= applies feature extractors to the power spectrum treated as a time series (frequency as time, power as magnitude). This captures the shape of the spectrum itself — useful for distinguishing sharp peaks from broad or noisy ones. Here we add Amplitude, Kurtosis, Skew, and BeyondNStd:

In [6]:

import light_curve as licu
import numpy as np

rng = np.random.default_rng(0)
PERIOD = 23.5
t = np.sort(rng.uniform(0, 300, 150))
m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)

pg_spec = licu.Periodogram(
    peaks=1,
    features=[licu.Amplitude(), licu.Kurtosis(), licu.Skew(), licu.BeyondNStd(nstd=2)],
)
res_spec = pg_spec(t, m)
print('Periodogram + spectrum features:')
for name, val in zip(pg_spec.names, res_spec):
    print(f'  {name:50s} = {val:.4f}')

import light_curve as licu import numpy as np rng = np.random.default_rng(0) PERIOD = 23.5 t = np.sort(rng.uniform(0, 300, 150)) m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150) pg_spec = licu.Periodogram( peaks=1, features=[licu.Amplitude(), licu.Kurtosis(), licu.Skew(), licu.BeyondNStd(nstd=2)], ) res_spec = pg_spec(t, m) print('Periodogram + spectrum features:') for name, val in zip(pg_spec.names, res_spec): print(f' {name:50s} = {val:.4f}')

Periodogram + spectrum features:
  periodogram_period_0                               = 23.6491
  periodogram_period_s_to_n_0                        = 12.5051
  periodogram_amplitude                              = 35.8933
  periodogram_kurtosis                               = 108.1316
  periodogram_skew                                   = 10.1212
  periodogram_beyond_2_std                           = 0.0127

Phase-folded features¶

phase_features= applies extractors to the light curve folded at the best period. The phase runs 0 → 1 (phase 0 = magnitude minimum). Features designed for smoothly varying series react very differently to a folded periodic signal versus the original unordered time series.

The LaflerKinmanStringLength $\theta$ is the normalised sum of squared successive differences of phase-sorted magnitudes — a clean, smooth periodic curve gives $\theta \ll 1$, while a noisy or aperiodic series gives $\theta \approx 1$. Similarly, EtaE (von Neumann $\eta^e$) is small for a random series and large for a phase-sorted periodic one, making it sensitive to regularity.

Extracting these without phase folding returns uninformative values; the table below shows the dramatic difference:

In [7]:

import light_curve as licu
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
PERIOD = 23.5
t = np.sort(rng.uniform(0, 300, 150))
m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)

pg_phase = licu.Periodogram(
    peaks=1,
    phase_features=[licu.LaflerKinmanStringLength(), licu.EtaE()],
)
result_pf = pg_phase(t, m)
best_p = result_pf[0]

# Compare with unfolded values
lk = licu.LaflerKinmanStringLength()
etae = licu.EtaE()
print(f'{'Feature':<35} {'non-folded':>12} {'phase-folded':>14}')
print('-' * 63)
print(f'{'string_length':<35} {lk(t, m)[0]:>12.4f} {result_pf[2]:>14.4f}')
print(f'{'eta_e':<35} {etae(t, m)[0]:>12.4f} {result_pf[3]:>14.4f}')
print()
print(f'Recovered period: {best_p:.3f} d  (true: {PERIOD} d)')

# Plot the phase-folded light curve
phase = (t % best_p) / best_p
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(phase, m, '.', alpha=0.5)
ax.set_xlabel('Phase')
ax.set_ylabel('Magnitude')
ax.invert_yaxis()
plt.xlim(0, 1)
ax.set_title(f'Phase-folded light curve  (P = {best_p:.2f} d)')
plt.tight_layout()
plt.grid()
plt.show()

import light_curve as licu import numpy as np import matplotlib.pyplot as plt rng = np.random.default_rng(0) PERIOD = 23.5 t = np.sort(rng.uniform(0, 300, 150)) m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150) pg_phase = licu.Periodogram( peaks=1, phase_features=[licu.LaflerKinmanStringLength(), licu.EtaE()], ) result_pf = pg_phase(t, m) best_p = result_pf[0] # Compare with unfolded values lk = licu.LaflerKinmanStringLength() etae = licu.EtaE() print(f'{'Feature':<35} {'non-folded':>12} {'phase-folded':>14}') print('-' * 63) print(f'{'string_length':<35} {lk(t, m)[0]:>12.4f} {result_pf[2]:>14.4f}') print(f'{'eta_e':<35} {etae(t, m)[0]:>12.4f} {result_pf[3]:>14.4f}') print() print(f'Recovered period: {best_p:.3f} d (true: {PERIOD} d)') # Plot the phase-folded light curve phase = (t % best_p) / best_p fig, ax = plt.subplots(figsize=(6, 3)) ax.plot(phase, m, '.', alpha=0.5) ax.set_xlabel('Phase') ax.set_ylabel('Magnitude') ax.invert_yaxis() plt.xlim(0, 1) ax.set_title(f'Phase-folded light curve (P = {best_p:.2f} d)') plt.tight_layout() plt.grid() plt.show()

Feature                               non-folded   phase-folded
---------------------------------------------------------------
string_length                             0.2437         0.0380
eta_e                                     6.5349       754.2290

Recovered period: 23.649 d  (true: 23.5 d)

Explicit frequency grid¶

Pass freqs= to fix the angular frequency grid (radians per time unit) instead of computing it from the data. This is useful when comparing results across different light curves on the same grid, or when you want fine control over the range.

For fast=True the grid must be np.linspace(0.0, max_omega, 2**k + 1):

In [8]:

import light_curve as licu
import numpy as np

rng = np.random.default_rng(0)
PERIOD = 23.5
t = np.sort(rng.uniform(0, 300, 150))
m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)

# Explicit grid: 0 to 2 rad/day in 2^10 + 1 steps (required for fast=True)
omega_max = 2.0  # rad/day
k = 10
fixed_freqs = np.linspace(0.0, omega_max, 2 ** k + 1)

pg_fixed = licu.Periodogram(peaks=1, freqs=fixed_freqs)
result_fixed = pg_fixed(t, m)
print(f'Best period (fixed grid): {result_fixed[0]:.3f} d')

import light_curve as licu import numpy as np rng = np.random.default_rng(0) PERIOD = 23.5 t = np.sort(rng.uniform(0, 300, 150)) m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150) # Explicit grid: 0 to 2 rad/day in 2^10 + 1 steps (required for fast=True) omega_max = 2.0 # rad/day k = 10 fixed_freqs = np.linspace(0.0, omega_max, 2 ** k + 1) pg_fixed = licu.Periodogram(peaks=1, freqs=fixed_freqs) result_fixed = pg_fixed(t, m) print(f'Best period (fixed grid): {result_fixed[0]:.3f} d')

Best period (fixed grid): 23.654 d

Multiband periodogram¶

When the same periodic signal is observed in several passbands, a joint Lomb–Scargle periodogram pools all bands into a single frequency grid. Bands with low amplitude still contribute, lifting the combined S/N above any single band.

Pass bands=[...] to Periodogram to switch to multiband mode. The optional multiband_normalization ('count' or 'chi2') controls how per-band powers are combined.

The example below synthesises a three-band light curve with two embedded periods (P₁ = 23.5 d and P₂ = 7.3 d). Each band has different amplitudes, magnitude offset, noise level, and a small inter-band phase shift (< 0.1 of the period):

band	amp(P₁)	amp(P₂)	offset	noise	phase(P₁)	phase(P₂)
g	0.35	0.10	15.0	0.04	0.00	0.00
r	0.22	0.18	15.5	0.15	0.03	0.05
i	0.12	0.28	15.9	0.06	0.07	0.08

The g band is dominated by P₁; the i band shows a stronger P₂ signal. The r band is noisy (noise > amp(P₂)), so it cannot recover P₂ on its own. No single band reliably recovers both periods, but the multiband periodogram does.

In [9]:

import light_curve as licu
import numpy as np
import matplotlib.pyplot as plt

BANDS = ["g", "r", "i"]
COLORS = {"g": "tab:green", "r": "tab:red", "i": "tab:purple"}
rng = np.random.default_rng(7)
P1, P2 = 23.5, 7.3

amp1 = {"g": 0.35, "r": 0.22, "i": 0.12}
amp2 = {"g": 0.10, "r": 0.18, "i": 0.28}
offset = {"g": 14.9, "r": 15.5, "i": 16.3}
phi1 = {"g": 0.00, "r": 0.03, "i": 0.07}
phi2 = {"g": 0.00, "r": 0.05, "i": 0.08}
noise = {"g": 0.04, "r": 0.15, "i": 0.06}

n_per_band = 80
t_list, m_list, s_list, b_list = [], [], [], []
for b in BANDS:
    t_b = np.sort(rng.uniform(0, 400, n_per_band))
    m_b = (offset[b]
           + amp1[b] * np.sin(2 * np.pi * (t_b / P1 - phi1[b]))
           + amp2[b] * np.sin(2 * np.pi * (t_b / P2 - phi2[b]))
           + rng.normal(0, noise[b], n_per_band))
    t_list.append(t_b)
    m_list.append(m_b)
    s_list.append(np.full(n_per_band, noise[b]))
    b_list.append(np.full(n_per_band, b))

t = np.concatenate(t_list)
m = np.concatenate(m_list)
sigma = np.concatenate(s_list)
band = np.concatenate(b_list)

fig, ax = plt.subplots(figsize=(8, 4))
for b in BANDS:
    mask = band == b
    phase = (t[mask] % P1) / P1
    ax.errorbar(phase, m[mask], yerr=sigma[mask],
                fmt='.', color=COLORS[b], alpha=0.6, ms=4, label=b)
ax.invert_yaxis()
ax.set_xlim(0, 1)
ax.set_xlabel('Phase')
ax.set_ylabel('mag')
ax.set_title(f'Phase-folded at P₁ = {P1} d')
ax.legend()
plt.tight_layout()
plt.show()

import light_curve as licu import numpy as np import matplotlib.pyplot as plt BANDS = ["g", "r", "i"] COLORS = {"g": "tab:green", "r": "tab:red", "i": "tab:purple"} rng = np.random.default_rng(7) P1, P2 = 23.5, 7.3 amp1 = {"g": 0.35, "r": 0.22, "i": 0.12} amp2 = {"g": 0.10, "r": 0.18, "i": 0.28} offset = {"g": 14.9, "r": 15.5, "i": 16.3} phi1 = {"g": 0.00, "r": 0.03, "i": 0.07} phi2 = {"g": 0.00, "r": 0.05, "i": 0.08} noise = {"g": 0.04, "r": 0.15, "i": 0.06} n_per_band = 80 t_list, m_list, s_list, b_list = [], [], [], [] for b in BANDS: t_b = np.sort(rng.uniform(0, 400, n_per_band)) m_b = (offset[b] + amp1[b] * np.sin(2 * np.pi * (t_b / P1 - phi1[b])) + amp2[b] * np.sin(2 * np.pi * (t_b / P2 - phi2[b])) + rng.normal(0, noise[b], n_per_band)) t_list.append(t_b) m_list.append(m_b) s_list.append(np.full(n_per_band, noise[b])) b_list.append(np.full(n_per_band, b)) t = np.concatenate(t_list) m = np.concatenate(m_list) sigma = np.concatenate(s_list) band = np.concatenate(b_list) fig, ax = plt.subplots(figsize=(8, 4)) for b in BANDS: mask = band == b phase = (t[mask] % P1) / P1 ax.errorbar(phase, m[mask], yerr=sigma[mask], fmt='.', color=COLORS[b], alpha=0.6, ms=4, label=b) ax.invert_yaxis() ax.set_xlim(0, 1) ax.set_xlabel('Phase') ax.set_ylabel('mag') ax.set_title(f'Phase-folded at P₁ = {P1} d') ax.legend() plt.tight_layout() plt.show()

In [10]:

pg_sb = licu.Periodogram(peaks=2, resolution=10, max_freq_factor=2.0, fast=True)
pg_mb = licu.Periodogram(peaks=2, bands=BANDS, resolution=10, max_freq_factor=2.0, fast=True)

fig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True)

for ax, b in zip(axes[:3], BANDS):
    mask = band == b
    omega, power = pg_sb.freq_power(t[mask], m[mask])
    pos = omega > 0
    periods = 2 * np.pi / omega[pos]
    ax.plot(periods, power[pos], lw=0.7, color=COLORS[b], label=b)
    for p_true in [P1, P2]:
        ax.axvline(p_true, color='k', lw=0.8, ls='--', alpha=0.4)
    ax.legend(loc='upper right')
    ax.set_xscale('log')
    ax.set_xlim(3, 60)

omega_mb, power_mb = pg_mb.freq_power(t, m, sigma, band)
pos = omega_mb > 0
periods_mb = 2 * np.pi / omega_mb[pos]
axes[3].plot(periods_mb, power_mb[pos], lw=0.7, color='k', label='multiband')
for p_true in [P1, P2]:
    axes[3].axvline(p_true, color='k', lw=0.8, ls='--', alpha=0.4)
axes[3].legend(loc='upper right')
axes[3].set_xlabel('Period (days)')
fig.suptitle('Lomb–Scargle power spectra  (dashed lines: true periods)')
plt.tight_layout()
plt.show()

print(f"{'':>12}  {'P₀ (d)':>8}  {'S/N₀':>6}  {'P₁ (d)':>8}  {'S/N₁':>6}")
print("-" * 50)
for b in BANDS:
    mask = band == b
    r = pg_sb(t[mask], m[mask])
    print(f"{b:>12}  {r[0]:>8.2f}  {r[1]:>6.2f}  {r[2]:>8.2f}  {r[3]:>6.2f}")
r_mb = pg_mb(t, m, sigma, band=band)
print("-" * 50)
print(f"{'multiband':>12}  {r_mb[0]:>8.2f}  {r_mb[1]:>6.2f}  {r_mb[2]:>8.2f}  {r_mb[3]:>6.2f}")
print(f"{'(true)':>12}  {P1:>8.1f}  {'—':>6}  {P2:>8.1f}  {'—':>6}")

pg_sb = licu.Periodogram(peaks=2, resolution=10, max_freq_factor=2.0, fast=True) pg_mb = licu.Periodogram(peaks=2, bands=BANDS, resolution=10, max_freq_factor=2.0, fast=True) fig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True) for ax, b in zip(axes[:3], BANDS): mask = band == b omega, power = pg_sb.freq_power(t[mask], m[mask]) pos = omega > 0 periods = 2 * np.pi / omega[pos] ax.plot(periods, power[pos], lw=0.7, color=COLORS[b], label=b) for p_true in [P1, P2]: ax.axvline(p_true, color='k', lw=0.8, ls='--', alpha=0.4) ax.legend(loc='upper right') ax.set_xscale('log') ax.set_xlim(3, 60) omega_mb, power_mb = pg_mb.freq_power(t, m, sigma, band) pos = omega_mb > 0 periods_mb = 2 * np.pi / omega_mb[pos] axes[3].plot(periods_mb, power_mb[pos], lw=0.7, color='k', label='multiband') for p_true in [P1, P2]: axes[3].axvline(p_true, color='k', lw=0.8, ls='--', alpha=0.4) axes[3].legend(loc='upper right') axes[3].set_xlabel('Period (days)') fig.suptitle('Lomb–Scargle power spectra (dashed lines: true periods)') plt.tight_layout() plt.show() print(f"{'':>12} {'P₀ (d)':>8} {'S/N₀':>6} {'P₁ (d)':>8} {'S/N₁':>6}") print("-" * 50) for b in BANDS: mask = band == b r = pg_sb(t[mask], m[mask]) print(f"{b:>12} {r[0]:>8.2f} {r[1]:>6.2f} {r[2]:>8.2f} {r[3]:>6.2f}") r_mb = pg_mb(t, m, sigma, band=band) print("-" * 50) print(f"{'multiband':>12} {r_mb[0]:>8.2f} {r_mb[1]:>6.2f} {r_mb[2]:>8.2f} {r_mb[3]:>6.2f}") print(f"{'(true)':>12} {P1:>8.1f} {'—':>6} {P2:>8.1f} {'—':>6}")

                P₀ (d)    S/N₀    P₁ (d)    S/N₁
--------------------------------------------------
           g     23.49   11.87     20.09    1.10
           r     23.53   10.94      8.27    2.09
           i      7.29   11.02     23.64    3.20
--------------------------------------------------
   multiband     23.60   15.92      7.30    5.33
      (true)      23.5       —       7.3       —

See also¶

• Feature table — full feature list
• Periodogram API — all constructor parameters