Periodogramas de Lomb–Scargle

1. Motivación: 40 noches, un planeta

Las velocidades radiales de una estrella se miden cuando se puede: de noche, sin nubes, cuando el objeto está alto y el telescopio asignado. El resultado es una serie temporal $\{t_i, y_i, \sigma_i\}$ con huecos, repeticiones y sin paso constante. La herramienta clásica, la transformada de Fourier discreta, presupone muestreo uniforme. El periodograma de Lomb–Scargle (Lomb 1976; Scargle 1982) resuelve el problema ajustando una sinusoide a cada frecuencia de prueba; es el método estándar para encontrar el periodo de un planeta (y de estrellas variables, rotación, binarias). Pero el muestreo deja huellas: picos fantasma (alias) y una estadística de falsos positivos que hay que entender para no «descubrir» planetas que no existen. Seguimos la exposición de VanderPlas (2018).

2. Del periodograma clásico al muestreo irregular

Para $N$ datos $y_n$ (de media cero) en $t_n$, el periodograma de Schuster es

Con muestreo uniforme $t_n = n\Delta t$ hereda las propiedades de Fourier: es periódico en $f$ con periodo $1/\Delta t$ y simétrico, de modo que solo tiene información hasta la frecuencia de Nyquist $f_{\rm Ny} = 1/(2\Delta t)$, y bajo ruido gaussiano $P_S$ sigue una distribución $\chi^2$ con dos grados de libertad en cada frecuencia (es la suma de dos gaussianas al cuadrado). Con muestreo irregular se pierden ambas cosas: no existe una $f_{\rm Ny}$ en sentido estricto (si los $t_n$ son múltiplos de una unidad $p$, por ejemplo la precisión con que se registran, entonces $f_{\rm Ny} = 1/2p$, que puede ser enorme; Eyer & Bartholdi 1999; VanderPlas 2018, §IV), y los términos seno y coseno dejan de ser independientes, lo que estropea la distribución estadística.

3. Derivación del periodograma de Lomb–Scargle

Ajuste de una sinusoide a cada frecuencia

En lugar de una transformada, Lomb y Scargle plantean un ajuste por mínimos cuadrados del modelo $y(t) = A\cos\omega t + B\sin\omega t$, $\omega = 2\pi f$, a los datos (centrados). La función a minimizar es $\chi^2(f) = \sum_n[y_n - A\cos\omega t_n - B\sin\omega t_n]^2$ y las ecuaciones normales $\partial\chi^2/\partial A = \partial\chi^2/\partial B = 0$ son

La reducción de $\chi^2$ al incluir la sinusoide es $\Delta\chi^2 = \chi^2_0 - \chi^2(f) = (YC, YS)\,\mathbf M^{-1}(YC, YS)^T$, con $\chi^2_0 = \sum y_n^2$ (problema 1).

El desfase τ

El término cruzado $CS$ se puede anular desplazando el origen de tiempos: con $\omega(t_n - \tau)$ en lugar de $\omega t_n$, $CS(\tau) = \tfrac12\sum\sin 2\omega(t_n - \tau) = \tfrac12[\cos 2\omega\tau\sum\sin 2\omega t_n - \sin 2\omega\tau\sum\cos 2\omega t_n] = 0$ si

Con esa elección la matriz es diagonal y $\Delta\chi^2 = YC^2/CC + YS^2/SS$ (todo evaluado en $t_n - \tau$). Scargle (1982) definió el periodograma como la mitad de esa cantidad:

(VanderPlas 2018, ecs. 33–34). Dos propiedades justifican el $\tau$: el periodograma es invariante ante traslaciones temporales y, bajo ruido gaussiano blanco de varianza $\sigma^2$, $P_{\rm LS}/\sigma^2$ sigue exactamente una distribución exponencial (χ² con 2 g.d.l.) en cada frecuencia, como el periodograma clásico. Para muestreo uniforme, $CC = SS = N/2$ y $P_{\rm LS} = P_S$. Para una sinusoide pura de amplitud $a$, $P_{\rm LS}(f_0) \approx Na^2/4$ (test de lombscargle.ts).

Periodograma generalizado y normalización

Dos mejoras prácticas (Zechmeister & Kürster 2009): pesos $w_n = 1/\sigma_n^2$ y una media flotante, es decir ajustar $y = c + A\cos\omega t + B\sin\omega t$ en lugar de restar la media muestral. La diferencia importa cuando la cobertura en fase es incompleta: la media de los datos no es la del modelo, y centrar en ella puede suprimir el pico verdadero (VanderPlas 2018, §VI.2). En ese caso las ecuaciones normales son $3\times3$; potenciaLS las resuelve eliminando $c$ analíticamente. Normalizamos como fracción de varianza explicada,

(VanderPlas 2018, ec. 59) con $P = 1$ para una sinusoide sin ruido y $P = 0$ si la sinusoide no mejora nada.

4. Ventana espectral y alias

El periodograma de los datos es, aproximadamente, el del proceso subyacente convolucionado con la ventana espectral, el periodograma del propio muestreo:

con $W(0) = 1$ (VanderPlas 2018, ec. 48). Cada pico de $W$ en $f_w \ne 0$ replica la señal verdadera en $f_{\rm obs} = |f_0 \pm nf_w|$ (ec. 45–47). Para observaciones nocturnas desde tierra el pico dominante es $f_w = 1$ d⁻¹; también aparecen el mes sinódico (29.5 d) y el año. La estructura de $W$ también fija la anchura de los picos, $\delta f \approx 1/T$ con $T$ la línea de base, y de ahí la rejilla recomendada:

con $f_{\max}$ elegido por física (el periodo más corto plausible) o por una «pseudo-Nyquist» del muestreo (VanderPlas 2018, §VII.1). Una rejilla uniforme en periodo submuestrea las frecuencias altas y es un error frecuente.

5. Significación: probabilidad de falsa alarma

La pregunta operativa es: si no hubiera señal, ¿con qué probabilidad el ruido produciría un pico tan alto como el observado? Es la probabilidad de falsa alarma (FAP). Para la normalización de §3 con media flotante y ruido gaussiano, en una frecuencia fijada de antemano (Cumming 2004; VanderPlas 2018, ec. 53),

Pero buscamos el máximo sobre muchas frecuencias. Con $N_{\rm eff}$ frecuencias independientes, $\mathrm{FAP} \approx 1 - [1 - \mathrm{FAP}_1(z)]^{N_{\rm eff}}$ y la estimación ingenua $N_{\rm eff} \approx f_{\max}T$ tiende a subestimar la FAP real. La cota de Baluev (2008), basada en la teoría de extremos de procesos continuos, no necesita $N_{\rm eff}$:

(VanderPlas 2018, ecs. 55–56). Es una cota superior válida para ventanas sin estructura fuerte; con alias intensos puede sobrestimar. La alternativa empírica es el bootstrap: permutar los $y_n$ sobre los mismos $t_n$ muchas veces, recalcular el periodograma y registrar el máximo; es costoso (para estimar una FAP $r$ hacen falta $\gtrsim 10/r$ periodogramas) y no trata ruido correlacionado. Implementación: fapUnaFrecuencia, fapBaluev, fapBootstrap.

Lo que la FAP no dice. No es la probabilidad de que la señal sea periódica ni de que el pico esté en la frecuencia correcta: un alias puede ser casi tan alto como el pico verdadero (§6). Y la incertidumbre del periodo no es la anchura del pico: la anchura ($\approx 1/T$) no depende de $N$ ni del S/N, mientras que la precisión sí, $\sigma_f \approx \frac{1}{T}\sqrt{\frac{2}{N\,\Sigma^2}}$ con $\Sigma$ el cociente señal/ruido medio (VanderPlas 2018, ec. 52), y en cualquier caso ignora el riesgo de estar en el pico equivocado.

Simulador

Genera una serie de velocidad radial tipo 51 Peg y cambia el patrón de muestreo de «aleatorio» a «solo de noche»: aparecen los alias en $f \pm 1$ c/d de la §4. Reduce $N$ o sube $\sigma$ hasta que el pico deje de superar la línea de FAP 1 %, y pliega los datos con el alias para ver por qué una curva plegada «bonita» no basta como prueba.

6. Ejemplo computacional: una serie tipo 51 Peg

Con velocidadRadial.ts y lombscargle.ts generamos 40 medidas en 40 noches distintas de un intervalo de 61 d (una medida por noche a hora aleatoria), con $P = 4.2293$ d, $K = 59$ m s⁻¹, órbita circular y ruido $\sigma = 13$ m s⁻¹ (ELODIE), rejilla $\Delta f = 1/(10T)$ hasta $f_{\max} = 2$ d⁻¹:

Lecciones del experimento: (i) el periodo se recupera dentro de la resolución de la rejilla; (ii) el pico es abrumadoramente significativo, pero (iii) el alias de 1.31 d tiene el 98 % de la potencia del verdadero: con una sola temporada de datos nocturnos no se puede decidir entre 4.23 d y 1.31 d solo con el periodograma (hacen falta observaciones a distintas horas, o desde otro hemisferio, o un segundo periodograma tras restar el candidato); (iv) en ruido puro el máximo alcanza 0.37, y Baluev (0.017) es más optimista que el bootstrap (0.033) porque la ventana está muy estructurada: para ventanas nocturnas conviene calibrar con remuestreo.

7. Análisis de error y coste

• Coste. El cálculo directo es $O(N\,N_{\rm eval})$ con $N_{\rm eval} = n_0Tf_{\max}$: para $T = 10$ años, $f_{\max} = 10$ d⁻¹ y $n_0 = 10$ son $3.65\times10^5$ frecuencias; el algoritmo de Press & Rybicki (1989) con FFT lo baja a $O(N_{\rm eval}\log N_{\rm eval})$.
• Error de rejilla. Si $\Delta f > 1/T$ los picos pueden caer entre nodos y perderse hasta la mitad de la potencia; con $n_0 = 5$ la pérdida máxima es $\approx 10\,\%$ (anchura sinc²).
• Ruido rojo. La actividad estelar produce potencia correlacionada a periodos de días a meses; ninguna de las FAP anteriores es válida, y los picos de rotación estelar (y sus alias) compiten con los planetas. Se modela con procesos gaussianos o se usan indicadores de actividad simultáneos.
• Múltiples señales. El periodograma es lineal en un solo sinusoide; dos planetas se buscan iterativamente («pre-blanqueado»: restar el mejor ajuste y repetir), o con ajustes multifrecuencia.

Problemas

1. (Demostración) Reducción de χ². Demuestra que para el modelo lineal $y = \sum_k\theta_k g_k(t)$ con matriz normal $\mathbf M$ y vector $\mathbf b = (\sum y_n g_k(t_n))$, el mínimo de $\chi^2$ es $\chi^2_0 - \mathbf b^T\mathbf M^{-1}\mathbf b$, y particulariza al modelo de dos parámetros (coseno y seno) con la matriz de §3.

2. (Demostración) El desfase τ. Demuestra que con $\tan 2\omega\tau = \sum\sin 2\omega t_n/\sum\cos 2\omega t_n$ se cumple $\sum_n\cos\omega(t_n - \tau)\sin\omega(t_n - \tau) = 0$, y que $P_{\rm LS}$ es invariante si todos los $t_n$ se desplazan en $t_0$.

3. Alias. Un planeta tiene $P = 3.6$ d. Con observaciones nocturnas (ventana en 1 d⁻¹) y a lo largo de un año (ventana en 1/365 d⁻¹), lista los periodos de los alias de primer orden. ¿Cuál de ellos se confundiría con una «supertierra en 1.4 d»? ¿Qué observación rompe la degeneración?

4. FAP. Con $N = 40$ y $f_{\max} = 2$ d⁻¹ sobre $T = 61$ d con $\mathrm{Var}(t) \approx T^2/12$: (a) calcula $\mathrm{FAP}_1$ y la cota de Baluev para $z = 0.3$, $0.5$ y $0.88$; (b) ¿qué altura $z$ hace falta para FAP = 1 %? (c) ¿Por qué la FAP de una frecuencia es tan pequeña comparada con la del máximo?

5. (Demostración) Distribución de una frecuencia. Para ruido gaussiano de varianza $\sigma^2$ y muestreo uniforme, demuestra que $P_{\rm LS}(f)/\sigma^2$ sigue una exponencial de media 1, es decir $\Pr(P_{\rm LS} > z\sigma^2) = e^{-z}$, en una frecuencia que no sea 0 ni Nyquist. ¿Por qué cambia a $(1 - z)^{(N-3)/2}$ con la normalización de §3?

6. (Computacional) Reproduce el ejemplo de §6. Genera la serie, calcula el periodograma generalizado, la ventana espectral y las tres FAP (una frecuencia, Baluev, bootstrap con 300 permutaciones) para el pico y para $z = 0.4$ en ruido puro. Luego resta la sinusoide ajustada y vuelve a calcular el periodograma: ¿desaparece el alias de 1.31 d?

import numpy as np
rng = np.random.default_rng(2024)
t = np.sort(rng.choice(np.arange(61), 40, replace=False) + 0.9 + 0.2*rng.random(40))
P0, K, s = 4.2293, 59.0, 13.0
y = K*np.cos(2*np.pi*t/P0) + s*rng.standard_normal(40)
def pot(t, y, f):               # periodograma generalizado, normalización [0,1]
    X = np.c_[np.ones_like(t), np.cos(2*np.pi*f*t), np.sin(2*np.pi*f*t)]
    r = y - X @ np.linalg.lstsq(X, y, rcond=None)[0]
    return 1 - (r @ r) / ((y - y.mean()) @ (y - y.mean()))
T = t.max() - t.min(); fr = np.arange(0, 2, 1/(10*T))
P = np.array([pot(t, y, f) for f in fr]); k = P.argmax()
W = np.abs(np.exp(-2j*np.pi*np.outer(fr, t)).sum(1))**2 / 40**2
N = 40; z = P[k]; Weff = 2*np.sqrt(4*np.pi*t.var())
fap1 = (1-z)**((N-3)/2); fap_b = 1 - (1-fap1)*np.exp(-Weff*(1-z)**((N-4)/2)*np.sqrt(z))
boot = np.mean([max(pot(t, rng.permutation(y), f) for f in fr) >= z for _ in range(300)])
print(1/fr[k], z, fap1, fap_b, boot, W[np.argmin(abs(fr-1))])

7. Rejilla. Para una campaña de $T = 8$ años con $f_{\max} = 1$ d⁻¹ y $n_0 = 10$, ¿cuántas frecuencias hay que evaluar? Si en lugar de ello se usa una rejilla uniforme en periodo entre 1 y 3000 d con el mismo número de puntos, ¿cuál es el paso en frecuencia cerca de $P = 1$ d y cerca de $P = 1000$ d? ¿Qué picos se pierden?

Lecturas

• VanderPlas, J. T. (2018). Understanding the Lomb–Scargle periodogram, ApJS 236, 16: la guía conceptual que sigue esta lección; incluye código.
• Scargle (1982), ApJ 263, 835 y Lomb (1976), Ap&SS 39, 447: los artículos originales; Zechmeister & Kürster (2009), A&A 496, 577: media flotante y pesos.
• Baluev (2008), MNRAS 385, 1279: la cota de la FAP; Cumming (2004), MNRAS 354, 1165: FAP y detectabilidad en sondeos de velocidad radial.
• Perryman (2018). The Exoplanet Handbook, §2.3; Gregory (2005). Bayesian Logical Data Analysis, cap. 13: el periodograma desde la perspectiva bayesiana (Bretthorst).