Luces Periodicas
21 Aug 2024La ciencia es, en su esencia, la búsqueda de patrones. A lo largo de la historia, hemos aprendido a reconocer estos patrones en todo lo que nos rodea, desde nuestro pulso hasta la órbitas de las estrellas alrededor de agujeros negros. Uno de los patrones más fundamentales y omnipresentes en la naturaleza es la periodicidad. Entender cómo y por qué algo se repite de manera regular nos permite no solo describir fenómenos naturales, sino también predecirlos.
Es que "asumir" que un ciclo se volverá a repetir trae grandes beneficios: Si asumimos que los días se repetiran uno tras otro tras otro sin pausa ni averías, tendremos ciudades con ritmos de vida frenéticos y personas con rutinas bien definidas, en función del horario en el que suceda este increible fenómeno físico llamado "órbita".
Este concepto de periodicidad es clave en muchos campos de la ciencias, y sin dudas la astronomía es uno de esos campos. Las estrellas, que desde la antigüedad han sido observadas como puntos de luz fijos en el cielo, en realidad nos envían señales que varían con el tiempo. Entre todos los tipos de variaciones que estas estrellas pueden experimentar, existen el cambio de brillo periódico, y su luz puede ser analizada para obtener información sobre la física de dicha estrella, su composición, su ciclo de vida, etc. Dicho de otra forma: hasta a las estrellas es posible robarle uno y que otro secreto.
¿Cómo cuantificamos la periodicidad de algo? Bueno, decimos que la función que depende del tiempo $f(t)$, que podría describir es periódica si existe un número positivo $P$, llamado período, tal que:
$$ f(t) = f(t + P) \quad \text{para cualquier tiempo } t. $$Es decir, que la función $f(t)$ cíclicamente va cambiando su valor y luego de un tiempo $P$ retoma su valor inicial, y así el ciclo se vuelve a repetir. Por ejemplo, pensemos en el "ritmo cardiaco": Cuando mides tu pulso, hay que medir cuántas veces se siente el pulso en un rango de tiempo de un minuto, para cuantificar los "latidos por minuto" o "pulsos por minuto" (ppm). Si consideramos que el pulso cambia repetitivamente en un periodo de tiempo, podemos describir la frecuencia cardiaca como una función que depende del tiempo $f(t)$. Entonces, bastaría medir la función en varios instantes $t$ para poder determinar esta frecuencia cardiaca.
Si consideramos a una persona quieta que está sentada y con un reloj de medición puesto, con un ppm=70 [pulsos/minuto] (es decir, que su corazón late 70 veces por minuto, que su latido dura un poco menos de un segundo), la medición de su frecuencia cardiaca se verá más o menos de la siguiente forma:
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('dark_background')
#Tamanos de las letras y puntos
fz=10
ms=4
# Parámetros para la simulación de la señal de frecuencia cardíaca
fs = 50 # Frecuencia de muestreo en Hz
t = np.arange(0, 5, 1/fs)
# Frecuencia cardíaca promedio en pulsos por minuto
frecuencia_cardiaca_ppm = 70.
# Convertir PPM (Pulsos Por Minuto) a Hz (pulsos por segundo)
frecuencia_hz = frecuencia_cardiaca_ppm/60.
# Generar la señal de frecuencia cardíaca (señal sintética tipo ECG)
senal_frecuencia_cardiaca = np.sin(2*np.pi * frecuencia_hz * t) + 0.5*np.sin(4*np.pi * frecuencia_hz * t)
# Agregamos ruido
senal_frecuencia_cardiaca += 0.2*np.random.randn(len(t))
plt.figure(figsize=(5, 2))
plt.plot(t, senal_frecuencia_cardiaca, 'o', ms=ms-2)
plt.title("Señal de frecuencia cardíaca simulada, con f = 1.116 [Hz]", fontsize=fz)
plt.xlabel("Tiempo [s]", fontsize=fz)
plt.ylabel("Amplitud normalizada", fontsize=fz-3)
plt.grid(color='gray')
# Ajustar el diseño de la gráfica
plt.tight_layout()
# Mostrar la gráfica
plt.show()
Podemos ver las subidas y bajadas de una señal periodica (o sea, que se repite en el tiempo), y al ojo podríamos estimar (contando el tiempo que separa máximos o mínimos consecutivos) que su periodo es alrededor de 0.9 segundos. Es decir, la función frecuencia cardiaca se vuelve a repetir en un periodo definido.
En este caso es posible saber cuál es periodo exacto de la señal simulada, dado que el periodo es el recíproco de los pulsos por segundo (o frecuencia [Hz]), de la forma:
$$ P = \frac{1}{frecuencia\ [Hz]} = \frac{1\ [s]}{1.116} = 0.857 [s] $$Determinado el periodo de una serie de tiempo, comienza lo interesante dado que podemos suponer cosas.
Como esta señal es periodica, significa que se está mapeando un proceso una y otra y otra vez a través del tiempo. Entonces, podríamos considerar todos los pulsos independientes y agruparlos sobre su periodo, con el fin de que la señal nos hable más sobre el procesos en cuestión. Este proceso se llama "doblado de fase" y asume que todos los ciclos de la señal son producto del mismo proceso sin cambios significativos.
La fase $\phi$ de una señal periódica describe la posición en el ciclo de la señal en un momento dado. Para calcular la fase de un fenómeno periódico se utiliza la expresión:
$$ \phi = \frac{t}{P} - int\left(\frac{t}{P}\right), $$donde $int()$ es la parte entera del argumento. Es decir, reducimos a un espacio entre 0 y 1, donde cada ciclo mapea el proceso en cuestión a través de sus números decimales. En la práctica se usa la función módulo "%" para calcular el resto, pero ahora lo haremos siguiendo la definición con fines ilustrativos.
pulso = 0.857 #Segundos
#Calculamos la fase a lo largo del tiempo:
fase = np.array([ (t[i]/pulso - int(t[i]/pulso)) for i in range(len(t)) ])
#Ordenamos la señal en función de la fase
sort = np.array([(senal_frecuencia_cardiaca[j],fase[j]) for j in range(len(t))])
sort = sort[np.argsort(sort[:,1])]
fase = sort[:,1]
senal_frec_ordenada = sort[:,0]
plt.figure(figsize=(5, 3))
plt.plot(fase, senal_frec_ordenada, '.', ms=ms)
plt.title("Señal de frecuencia cardíaca con P=0.857 [s]", fontsize=fz)
plt.xlabel("Espacio de fase $\phi$", fontsize=fz)
plt.ylabel("Amplitud normalizada", fontsize=fz)
plt.tight_layout()
plt.show()
Si una señal tiene un comportamiento periódico, es conveniente describirla en términos de su "espacio de fase", nombre que se le da al gráfico entre el parámetro en cuestión y la fase $\phi$, desde donde podremos interpretar la forma de su ciclo, y qué parámetros podrían influir en su forma. En astrofísica, medir este tipo de pulsos puede hablarnos sobre la densidad de los objetos, sobre el mecanismo físico detrás de esta señal
Señales desde las estrellas
Imagina que estás observando una estrella que cambia su brillo de forma periódica. Si realizaras fotometría a sus observaciones, medirías el cambio de luminosidad en el tiempo, es posible obtiene una señal periodica en formato de serie de tiempo $\{x_n\}$. Es decir, una colección de observaciones de cierto suceso que ocurre en el cielo.
En astronomía este concepto se adapta para analizar la evolución temporal del brillo de las estrellas (galaxias también pero es más complejo). En general, para cada estrella conocida existe una serie de tiempo $\{x_n\}$ que es posible estudiar y determinar qué eventos ocurren en dicha estrella. Entonces, en nuestro caso, asumiendo que podría existir un periodo sus series de tiempo son escudriñadas para encontrar dichos periodos.
Veamos un caso de la vida real: la estrella periodica 2MASS J17222333-3644037 ha sido reportada como una Cepheida de tipo II. Yo la he encontrado de forma independiente, y la he llamado SFR3_pv12 (prefiero eso a que 2MASS J17222333-3644037).
SFR3_pv12 es una estrella Cefeida de tipo II, una de esas estrellas que nos ayuda a determinar distancias cósmicas. El telescopio VISTA ha hecho 75 mediciones de la luz de este objeto astronómico a lo largo de aproximadamente 5 años (1752 días), y con esta información se ha generado la siguiente serie de tiempo $\{x_n\}$, de la magnitud infrarroja de SFR3_pv12 versus el tiempo:
lc = np.genfromtxt("SFR3_pv12.lc")
mjd,mag,err = lc[:,0],lc[:,1],lc[:,2]
#Normalizamos el tiempo
mjd0 = min(mjd)-1
MJD=mjd-mjd0
#
#MJD = MJD[1:]
#mag = mag[1:]
#err = err[1:]
plt.figure(figsize=(6,2.5))
plt.errorbar(MJD, mag, err,0, 'ro', ecolor='gray', alpha=0.8, ms=ms-1)
plt.xlabel('MJD-%0.1f [d]' %mjd0,fontsize=fz)
plt.ylabel("$K_{s}$ [mag]",fontsize=fz)
plt.title("Serie de tiempo de la estrella SFR3_pv12", fontsize=fz)
ran = np.linspace(np.min(mag)-np.mean(err),np.max(mag)+np.mean(err),4)
rango = np.array([round(ran[i],2) for i in range(len(ran))])
plt.yticks(rango,fontsize=fz)
plt.xticks(fontsize=12)
plt.gca().invert_yaxis()
plt.show()
He ahí la serie de tiempo de una estrella Cefeida tipo II. No parece nada, solo ruido. Aquí no es evidente que exista siquiera un periodo detrás de esta señal. Puede que esto se deba a varios factores: Puede ser el rango temporal en el que están mapeados los ciclos, puede que no exista suficiente cadencia, entre otros.
Estimado público, bienvenido a los datos del mundo real! No pasa nada, solo tenemos que hacer algunas operaciones más y esto tendrá algo de sentido. Usando métodos de análisis armónicos, hemos estimado que el periodo de esta estrella es de $P=14.55$ días, es decir, que cada 14.55 días esta estrella cumple un ciclo.
Entonces, cuando usamos la transformación al espacio de fase $\phi$ que comentamos anteriormente:
p = 14.55
phase = np.array([ (MJD[i]/p - int(MJD[i]/p)) for i in range(len(MJD)) ])
sort = np.array([(mag[j],phase[j],err[j]) for j in range(len(mag))])
sort = sort[np.argsort(sort[:,1])]
phase = sort[:,1]
mag_sort = sort[:,0]
err_sort = sort[:,2]
plt.figure(figsize=(4,3))
plt.title("Curva de luz de SFR3_pv12, con P=14.55 [d]", fontsize=fz)
plt.errorbar(phase, mag_sort, err_sort,0, 'ro', ecolor='gray', alpha=0.8, ms=ms)
plt.xticks([0,0.5,1],fontsize=fz)
plt.xlabel('$\phi$',fontsize=fz)
plt.ylabel("$K_{s}$ [mag]",fontsize=fz)
holi = np.linspace(np.min(mag)-np.mean(err),np.max(mag)+np.mean(err),4)
rango = np.array([round(holi[i],2) for i in range(len(holi))])
plt.yticks(rango,fontsize=fz)
plt.gca().invert_yaxis()
plt.show()
¡Voila! Efectivamente la estrella tiene un periodo definido, dado el orden o la coherencia en su espacio de fase. Esta propiedad es muy importante, ya que existen una serie de métodos que utilizan esta coherencia para determinar el periodo de series de tiempo astronómicas, con el fin de poder identificar el mayor número posible de estrellas que compartan cierto tipo de propiedades (digamos, Cefeidas de tipo II) y con eso, poder estructurar un argumento científico que esgrimir.
Sólo una cosa queda dando vueltas: ¿Cómo se obtiene el periodo, entonces? Esta es la real pregunta. Uno de los métodos más famosos y utilizados es el Periodograma.
El análisis espectral
El análisis espectral es una técnica utilizada para descomponer una señal en sus componentes de frecuencia. Permite estudiar cómo la energía de una señal se distribuye a lo largo del espectro de frecuencias. Dentro de las herramientas de este tipo de análisis está el Periodograma, muy efectivo para caracterizar frecuencias dominantes en series de tiempo mal muestreadas.
Si existe una frecuencia angular $\omega$ dentro de $\{x_n\}$, el periodograma debería ser capaz de encontrala, maximizando el poder espectral $P(\omega)$. Para esto, debemos elegir una grilla de frecuencias y podemos calcular el poder $P(\omega)$ usando la siguiente definición (de entre una multitud de definiciones):
$$ P(\omega) = \left( \sum_i^n \left(m_i - \overline{m}\right) \cdot \cos(\omega t_i)\right)^2 + \left( \sum_i^n \left(m_i - \overline{m}\right) \cdot \sin(\omega t_i)\right)^2 $$donde $m_i$ son las magnitudes y $t_i$ son los tiempos dentro de la serie de tiempo $\{x_n\}$. El Periodograma debe ser máximo cuando identifica una frecuencia representativa.
muestreo = 500 #cadencia
#Rango en la búsqueda de periodos
P_i = 10. # Periodo inicial
P_f = 16. # Periodo final
# Calcular el periodograma inicial
frequencies = np.linspace(1/P_f, 1/P_i, muestreo) # Frecuencias (1/Periodo)
powers = np.zeros_like(frequencies) # Potencia asociada a cada frecuencia
# Calculamos la ecuación de arriba:
for j, freq in enumerate(frequencies):
omega = 2 * np.pi * freq
powers[j] = np.sum((mag - np.mean(mag))*np.cos(omega * MJD))**2 + np.sum((mag - np.mean(mag))*np.sin(omega * MJD))**2
# Encontrar el índice de la potencia máxima
max_index = np.argmax(powers)
freq_max = frequencies[max_index] # Frecuencia que maximiza powers
# Calcular el periodo correspondiente
periodo = 1 / freq_max
periodos = 1 / frequencies
print("El periodo que maximiza el periodograma es %0.2f días:" %periodo)
plt.figure(figsize=(5,3))
plt.plot(periodos, powers, 'b-') # Línea del periodograma
plt.plot(periodo, np.max(powers),'r*', ms=ms+8)
plt.xlabel('Periodo [dias]', fontsize=fz)
plt.ylabel('Poder espectral $P(\omega)$', fontsize=fz-2)
plt.show()
Existe muchas implementaciones de diferentes periodogramas, uno de los más populares es el Periodograma Generalizado de Lomb-Scargle que sin duda realiza un buen trabajo para series de tiempo "suficientemente" muestreadas.
Finalmente, haré una animación para ver cómo las frecuencias que maximizan el periodograma es, a su vez, el que da coherencia al espacio de fase.
Un pseudo algoritmo para implementar la búsqued de periodos en una fuente astronomíca se describe a continuación:
- Identificar y obtener la serie de tiempo $\{x_n\}$ de la fuente sospechosa de periodicidad.
- Crear una grilla de posibles frecuencias (o periodos) donde uno sospecha que podría encontrarse el periodo.
- Aplicar el periodograma, calculando el poder espectral $P(\omega)$.
- Identidicar el periodo que maximiza $P(\omega)$
- Convertir al espacio de fase y discutir caso a caso.
Como conclusión, destaco la importancia de los procesos cíclicos en la naturaleza, una suposición poderosa para predecir el comportamiento de pequeños trozos del cosmos.
muestreo = 350 # Número de cuadros en la animación
P_i = 13.0 # Periodo inicial
P_f = 15.0 # Periodo final
period_step = (P_f-P_i)/muestreo # Cambio del periodo en cada iteración
# Calcular el periodograma inicial
frequencies = np.linspace(1/P_f, 1/P_i, muestreo) # Frecuencias (1/Periodo)
powers = np.zeros_like(frequencies) # Potencia asociada a cada frecuencia
# Configuración del gráfico
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(3.5, 4.5), dpi=100, gridspec_kw={'height_ratios': [1, 2]})
# Configuración del gráfico de periodograma
ax1.set_xlabel('periodo [días]', fontsize=fz-3)
ax1.set_ylabel('Poder espectral $P(\omega)$', fontsize=fz-2)
ax1.xaxis.set_label_position('top') # Mover la etiqueta del eje X arriba
ax1.xaxis.tick_top() # Mover los ticks del eje X arriba
ax1.tick_params(axis='y', which='both', left=False, right=False, labelleft=False)
# Configuración del gráfico de fase vs magnitud
ax2.set_xlabel('$\phi$', fontsize=fz)
ax2.set_ylabel("$K_{s}$ [mag]", fontsize=fz)
ax2.invert_yaxis()
ax2.set_xticklabels([])
# Dibujar el periodograma inicial
for j, freq in enumerate(frequencies):
omega = 2 * np.pi * freq
powers[j] = np.sum((mag - np.mean(mag)) * np.cos(omega * MJD)) ** 2 + np.sum((mag - np.mean(mag)) * np.sin(omega * MJD)) ** 2
periodos = 1/frequencies
# Línea del periodograma
line1, = ax1.plot(periodos, powers, 'b-')
# Mantener los ejes fijos
ax1.set_xlim([min(periodos), max(periodos)])
# Línea vertical inicial
vline = ax1.axvline(x=min(periodos), color='w', linestyle='--', alpha=0.5)
# Configurar los ticks del eje X
ticks_x = np.linspace(periodos.min(), periodos.max(), 6) # 4 valores distribuidos uniformemente
ax1.set_xticks(ticks_x)
ax1.set_xticklabels([f'{tick:.1f}' for tick in ticks_x], fontsize=fz-2)
# Función de actualización para la animación
def animate(i):
# Actualizar el periodo y la frecuencia
p = P_i + i * period_step
freq = 1 / p # Convertir periodo en frecuencia
# Actualizar la línea vertical en el periodograma
# Esto lo aprendí haciendo esto y me parece increible!!
vline.set_xdata(p)
# Calcular la fase
phase = np.array([ (MJD[i]/p - int(MJD[i]/p)) for i in range(len(MJD)) ])
sort = np.array([(phase[j],mag[j],err[j]) for j in range(len(mag))])
sort = sort[np.argsort(sort[:,0])]
phase = sort[:,0]
mag_sort = sort[:,1]
err_sort = sort[:,2]
# Actualizar los datos de la gráfica de fase vs magnitud
ax2.clear() # Limpia solo los datos anteriores, no los ejes
ax2.set_xlabel('$\phi$', fontsize=fz+6)
ax2.set_ylabel("$K_{s}$ [mag]", fontsize=fz)
# Ocultar los ticks y las líneas en los ejes
ax2.tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False)
ax2.tick_params(axis='y', which='both', left=False, right=False, labelleft=False)
ax2.invert_yaxis()
ax2.errorbar(phase, mag_sort, err_sort, 0, 'ro', ecolor='gray', alpha=0.8, ms=ms)
ax2.set_title(f"Periodo = {p:.2f} días", fontsize=fz+1)
# Mantener los ejes fijos
ax2.set_xlim([0, 1])
# Configuración de la animación
ani = animation.FuncAnimation(fig, animate, frames=muestreo, interval=200, blit=False)
# Guardar la animación como un archivo de video
ani.save('SFR3_pv12.mp4', writer='ffmpeg', fps=40, dpi=300)
plt.show()
El video del resultado puede ser visto en el siguiente video de Youtube:
from IPython.display import HTML
# Reemplaza '1eSr9jzxmaM' con el ID real del YouTube Short
video_id = '1eSr9jzxmaM'
# Crear la URL del iframe
html_string = f'''
<iframe width="560" height="315" src="https://www.youtube.com/embed/{video_id}" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
'''
# Mostrar el video en el notebook
HTML(html_string)