Derivar LOSVD del espectro de galaxias, usando el método de Fourier

Estoy tratando de implementar un programa que derive la cinemática (específicamente los parámetros cinemáticos: velocidad de rotación media, velocidad, dispersión, coeficientes de hermita h3 y h4) de un espectro de galaxias elípticas. Sé que esto se ha implementado de muchas maneras, la más común en los últimos tiempos es ppxf, pero quiero hacerlo con fines educativos. Empecé trabajando con un espectro (nombre de variable = especificación) de una estrella K2III (dado que son bastante comunes en las galaxias de tipo temprano), eliminé el continuo del espectro y lo reclasifiqué a un tamaño de longitud de onda logarítmica. Luego creé un espectro de galaxias artificiales (variablename = spec_gal), ampliando mi espectro estelar con una pérdida gaussiana con una dispersión de 200 km/s. Por lo tanto, acabo de calcular la convolución entre el espectro estelar y losvd usando np.convolve().

spec_gal = np.convolve(losvd_gauss, spec, mode='same')

ingrese la descripción de la imagen aquí

ingrese la descripción de la imagen aquí¿Por qué los valores de flujo se vuelven tan pequeños cuando convoluciono la pérdida con el espectro? Debería ampliar las líneas de absorción, ¿no? Supongo que es porque multiplico pequeños valores de mi losvd con los valores espectrales.

Como una pequeña prueba, sé que quería recuperar mi losvd del espectro de galaxias sintéticas usando la transformación de Fourier. Como la convolución es multiplicativa en el espacio de Fourier, pensé, podría calcular las transformadas de Fourier de mis espectros usando scipy.fft().

spec_fourier = fft(spec)
spec_gal_fourier = fft(spec_gal)

Y derivar el losvd en el espacio de Fourier dividiendo el espectro de mi galaxia en el espacio de Fourier por el espectro estelar.

 losvd_fourier = spec_gal_fourier/spec_fourier

Tracé los valores np.abs() de mis espectros transformados de Fourier. ingrese la descripción de la imagen aquíPensé que podría recuperar mi losvd simplemente realizando la transformada inversa de Fourier usando scipy.ifft().

losvd = ifft(losvd_fourier)

Lo cual no me está dando el resultado correcto (vea el último gráfico a continuación) y agradecería comentarios, comentarios y sugerencias sobre posibles errores en la teoría y la implementación.ingrese la descripción de la imagen aquí

Estoy seguro de que obtendrá una respuesta aquí, pero solo para su información, también hay un sitio de Scientific Computing SE que es excelente para preguntas numéricas. Por cierto, ¿qué sucede si primero normaliza su núcleo de convolución gaussiana a la unidad de área? Por ejemplo g = g / g.sum(), hice una prueba rápida con algunos datos simples y eso parece estar al menos cerca de lo correcto.
De hecho, me pregunto por qué mi kernel gaussiano no está normalizado en primer lugar, ya que he definido mi losvd gaussiano con la normalización estándar de gauss N = 1/sqrt(2pi*sigma^2). Su sugerencia resuelve la pregunta con los pequeños valores para mi espectro de galaxias sintéticas, pero el problema final sigue en pie.
Es posible que haya dejado accidentalmente el σ 2 fuera de la raíz cuadrada. De todos modos, estoy seguro de que alguien podrá abordar su pregunta principal.

Respuestas (1)

En primer lugar, hubo algunos errores en el tratamiento de mis espectros antes de la transformada de Fourier. Así que los pasos que faltaban eran:

1.) Los espectros deben promediarse a 0 restando la mediana.

2.) Debe tocar los espectros en los bordes con una función de ventana (por ejemplo, np.blackman())

3.) Cortar píxeles en ambos extremos de los espectros para evitar discontinuidades

Luego recupera la entrada correcta LOSVD, utilizando el método de deconvolución de división simple en el espacio de Fourier.

¡Felicidades! Siempre está bien responder a su propia pregunta, y si pasa una semana y no aparecen más respuestas (como es probable en este caso), simplemente continúe y acepte su propia respuesta.