Compare directamente los cambios de subpíxeles entre dos espectros y obtenga errores creíbles

9

Tengo dos espectros del mismo objeto astronómico. La pregunta esencial es esta: ¿cómo puedo calcular el cambio relativo entre estos espectros y obtener un error preciso en ese cambio?

Algunos detalles más si todavía estás conmigo. Cada espectro será una matriz con un valor x (longitud de onda), valor y (flujo) y error. El cambio de longitud de onda va a ser subpíxel. Suponga que los píxeles están espaciados regularmente y que solo se aplicará un cambio de longitud de onda a todo el espectro. Entonces, la respuesta final será algo así como: 0.35 +/- 0.25 píxeles.

Los dos espectros serán una gran cantidad de características continuas sin características puntuadas por algunas características de absorción (inmersiones) bastante complicadas que no se modelan fácilmente (y no son periódicas). Me gustaría encontrar un método que compare directamente los dos espectros.

El primer instinto de todos es hacer una correlación cruzada, pero con los cambios de subpíxeles, tendrás que interpolar entre los espectros (¿suavizando primero?); Además, los errores parecen desagradables para acertar.

Mi enfoque actual es suavizar los datos convolucionando con un núcleo gaussiano, luego dividir el resultado suavizado y comparar los dos espectros divididos, pero no confío en ello (especialmente los errores).

¿Alguien sabe de una manera de hacer esto correctamente?

Aquí hay un breve programa de Python que producirá dos espectros de juguete que se desplazan por 0.4 píxeles (escritos en toy1.ascii y toy2.ascii) con los que puedes jugar. Aunque este modelo de juguete utiliza una característica gaussiana simple, suponga que los datos reales no pueden ajustarse a un modelo simple.

import numpy as np
import random as ra
import scipy.signal as ss
arraysize = 1000
fluxlevel = 100.0
noise = 2.0
signal_std = 15.0
signal_depth = 40.0
gaussian = lambda x: np.exp(-(mu-x)**2/ (2 * signal_std))
mu = 500.1
np.savetxt('toy1.ascii', zip(np.arange(arraysize), np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth), np.ones(arraysize) * noise))
mu = 500.5
np.savetxt('toy2.ascii', zip(np.arange(arraysize), np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth), np.ones(arraysize) * noise))
JBWhitmore
fuente
Si entiendo correctamente, el problema suena similar al registro de imágenes, excepto que solo tiene un desplazamiento lineal de subpíxeles en un eje. ¿Quizás intente técnicas estándar de registro de imágenes como la correlación de fase?
Paul R
Si tiene un retraso puro en una señal (es decir, el cambio en el parámetro de longitud de onda del que habla), podría explotar la propiedad de transformación de Fourier que convierte el retraso de tiempo en un desfase lineal en el dominio de frecuencia. Esto podría funcionar si las dos muestras no están dañadas por diferentes interferencias o ruidos de medición.
Jason R
1
Este hilo puede ser útil- dsp.stackexchange.com/questions/2321/…
Jim Clay
1
¿Tiene datos reales para probar? El valor de ruido que proporcionó es demasiado para que la correlación cruzada sea precisa de submuestra. Esto es lo que encuentra con varias corridas de ruido 2.0 y offset 0.7 (= 1000.7 en el eje x de la gráfica), por ejemplo: i.stack.imgur.com/UK5JD.png
endolito

Respuestas:

5

Creo que usar correlación cruzada e interpolar el pico funcionaría bien. Como se describe en ¿Es inútil el muestreo previo antes de la correlación cruzada? , interpolar o muestrear antes de que la correlación cruzada en realidad no le brinde más información. La información sobre el pico de submuestra está contenida en las muestras a su alrededor. Solo necesita extraerlo con un error mínimo. Reuní algunas notas aquí .

El método más simple es la interpolación cuadrática / parabólica, de la cual tengo un ejemplo de Python aquí . Supuestamente es exacto si su espectro se basa en una ventana gaussiana , o si el pico cae exactamente en el punto medio entre muestras, pero por lo demás tiene algún error . Entonces, en su caso, probablemente quiera usar algo mejor.

Aquí hay una lista de estimadores más complicados, pero más precisos. "De los métodos anteriores, el segundo estimador de Quinn tiene el menor error RMS".

No sé las matemáticas, pero este artículo dice que su interpolación parabólica tiene una precisión teórica del 5% del ancho de un contenedor FFT.

El uso de la interpolación FFT en la salida de correlación cruzada no tiene ningún error de sesgo , por lo que es lo mejor si desea una precisión realmente buena. Si necesita equilibrar la precisión y la velocidad de cálculo, se recomienda hacer una interpolación FFT y luego seguirla con uno de los otros estimadores para obtener un resultado "suficientemente bueno".

Esto solo usa el ajuste parabólico, pero genera el valor correcto para el desplazamiento si el ruido es bajo:

def parabolic_polyfit(f, x, n):
    a, b, c = polyfit(arange(x-n//2, x+n//2+1), f[x-n//2:x+n//2+1], 2)
    xv = -0.5 * b/a
    yv = a * xv**2 + b * xv + c

    return (xv, yv)

arraysize = 1001
fluxlevel = 100.0
noise = 0.3 # 2.0 is too noisy for sub-sample accuracy
signal_std = 15.0
signal_depth = 40.0
gaussian = lambda x: np.exp(-(mu-x)**2/ (2 * signal_std))
mu = 500.1
a_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)
mu = 500.8
b_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)

a_flux -= np.mean(a_flux)
b_flux -= np.mean(b_flux)

corr = ss.fftconvolve(b_flux, a_flux[::-1])

peak = np.argmax(corr)
px, py = parabolic_polyfit(corr, peak, 13)

px = px - (len(a_flux) - 1)
print px

ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí

El ruido en su muestra produce resultados que varían en más de una muestra completa, así que lo reduje. Ajustar la curva usando más puntos del pico ayuda a ajustar un poco la estimación, pero no estoy seguro de si eso es estadísticamente válido, y en realidad empeora la estimación para la situación de menor ruido.

Con ruido = 0.2 y ajuste de 3 puntos, da valores como 0.398 y 0.402 para desplazamiento = 0.4.

Con ruido = 2.0 y ajuste de 13 puntos, da valores como 0.156 y 0.595 para desplazamiento = 0.4.

endolito
fuente
Estoy tratando de resolver este problema exacto para el registro de imágenes. Necesito precisión de subpíxel (0.1 probablemente sería lo suficientemente bueno) pero lo más importante es que no necesito sesgo, por lo que los métodos de interpolación no funcionan. ¿Hay algún método bueno (e implementado en Python?) Para esto? El método de relleno cero funcionará, pero es demasiado costoso para ser práctico.
keflavich
@kelavich: ¿Ha probado todos los métodos de interpolación y ha encontrado un sesgo inaceptable? El enfoque recomendado es una combinación de algunos de cero relleno seguida de una interpolación de bajo error. No conozco ningún otro método, pero apuesto a que esto le proporcionaría mucha precisión.
endolito
Sí, he encontrado un sesgo inaceptable en la interpolación lineal y de segundo orden. He intentado el relleno cero de FFT, pero el resultado está dominado por el timbre de alta frecuencia ... ¿hay alguna posibilidad de que pueda agregar un ejemplo de relleno cero?
keflavich