Bezpośrednio porównaj przesunięcia subpikseli między dwoma widmami i uzyskaj wiarygodne błędy

Mam dwa widma tego samego obiektu astronomicznego. Zasadnicze pytanie brzmi: jak obliczyć względne przesunięcie między tymi widmami i uzyskać dokładny błąd na tym przesunięciu?

Więcej szczegółów, jeśli nadal jesteś ze mną. Każde widmo będzie tablicą o wartości x (długość fali), wartości y (strumień) i błędzie. Przesunięcie długości fali będzie subpikselowe. Załóżmy, że piksele są regularnie rozmieszczone i że do całego spektrum zostanie zastosowane tylko jedno przesunięcie długości fali. Tak więc odpowiedź końcowa będzie taka: 0,35 +/- 0,25 pikseli.

Te dwa widma będą miały wiele pozbawionych cech kontinuum przerywanych pewnymi dość skomplikowanymi cechami absorpcji (zapadów), które nie dają się łatwo modelować (i nie są okresowe). Chciałbym znaleźć metodę, która bezpośrednio porównuje dwa widma.

Pierwszym instynktem każdego człowieka jest wykonanie korelacji krzyżowej, ale przy przesunięciach subpikseli będziesz musiał interpolować między widmami (najpierw wygładzając?) - również błędy wydają się paskudne, aby naprawić błąd.

Moje obecne podejście polega na wygładzaniu danych przez zwojenie jądra gaussowskiego, a następnie na spline wygładzony wynik i porównaniu dwóch widm splajnowych - ale nie ufam temu (szczególnie błędom).

Czy ktoś wie, jak to zrobić poprawnie?

Oto krótki program w języku Python, który wytworzy dwa widma zabawki, które są przesunięte o 0,4 piksela (zapisane w toy1.ascii i toy2.ascii), z którymi możesz grać. Mimo że ten model zabawki wykorzystuje prostą funkcję gaussowską, załóż, że rzeczywistych danych nie można dopasować do prostego modelu.

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))

Myślę, że użycie korelacji krzyżowej i interpolacja wartości szczytowej działałoby dobrze. Jak opisano w Czy pobieranie próbek przed korelacją krzyżową jest bezużyteczne? , interpolacja lub upsampling przed korelacją krzyżową w rzeczywistości nie daje żadnych dodatkowych informacji. Informacje o piku podpróbki znajdują się w otaczających go próbkach. Wystarczy go wyodrębnić przy minimalnym błędzie. Zebrałem tutaj notatki .

Najprostszą metodą jest kwadratowa / paraboliczne interpolacji, który mam przykład Pythona tutaj . Jest to prawdopodobnie dokładne, jeśli twoje widmo opiera się na oknie Gaussa , lub jeśli pik przypada dokładnie na środkowy punkt między próbkami, ale poza tym ma jakiś błąd . Więc w twoim przypadku prawdopodobnie chcesz użyć czegoś lepszego.

Oto lista bardziej skomplikowanych, ale dokładniejszych estymatorów. „Spośród powyższych metod drugi estymator Quinna ma najmniejszy błąd RMS.”

Nie znam matematyki, ale ten artykuł mówi, że ich interpolacja paraboliczna ma teoretyczną dokładność 5% szerokości pojemnika FFT.

Użycie interpolacji FFT na wyjściu korelacji krzyżowej nie powoduje błędu błędu , więc jest to najlepsze, jeśli chcesz naprawdę dobrej dokładności. Jeśli potrzebujesz zrównoważyć dokładność i szybkość obliczeń, zaleca się wykonanie interpolacji FFT, a następnie śledzenie jej za pomocą jednego z innych estymatorów, aby uzyskać wynik „wystarczająco dobry”.

To po prostu wykorzystuje dopasowanie paraboliczne, ale wyświetla odpowiednią wartość dla przesunięcia, jeśli hałas jest niski:

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

Hałas w twojej próbce daje wyniki, które różnią się o więcej niż całą próbkę, więc ją zmniejszyłem. Dopasowanie krzywej przy użyciu większej liczby punktów piku pomaga nieco zawęzić oszacowanie, ale nie jestem pewien, czy jest to poprawne statystycznie, a tak naprawdę pogarsza oszacowanie w sytuacji niższego szumu.

Przy szumie = 0,2 i 3-punktowym dopasowaniu daje wartości takie jak 0,398 i 0,402 dla przesunięcia = 0,4.

Przy szumie = 2,0 i 13-punktowym dopasowaniu daje wartości takie jak 0,156 i 0,595 dla przesunięcia = 0,4.