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

9

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))
JBWhitmore
źródło
Jeśli dobrze rozumiem, problem brzmi podobnie do rejestracji obrazu, z tym wyjątkiem, że masz tylko liniowe przesunięcie subpikseli w jednej osi. Może wypróbujesz standardowe techniki rejestracji obrazów, takie jak korelacja faz?
Paul R
Jeśli masz czyste opóźnienie w jednym sygnale (tj. Przesunięcie parametru długości fali, o którym mówisz), możesz być w stanie wykorzystać właściwość transformacji Fouriera, która zamienia opóźnienie czasowe w liniowe przesunięcie fazowe w dziedzinie częstotliwości. Może to działać, jeśli dwie próbki nie są uszkodzone przez inny szum pomiarowy lub zakłócenia.
Jason R
1
Ten wątek może być przydatny
Jim Clay
1
Czy masz rzeczywiste dane do przetestowania? Podana wartość szumu jest zbyt duża, aby korelacja krzyżowa była dokładna dla podpróbki. To, co znajduje się z kilku przebiegów hałasu 2.0 i 0.7 (przesunięcie = 1000.7 na wykresie na oś X), na przykład: i.stack.imgur.com/UK5JD.png
endolit

Odpowiedzi:

5

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

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

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.

endolit
źródło
Próbuję rozwiązać ten problem dotyczący rejestracji obrazu. Potrzebuję dokładności subpikselowej (0,1 byłaby prawdopodobnie wystarczająca), ale co najważniejsze nie potrzebuję odchylenia, więc metody interpolacji nie działają. Czy istnieją jakieś dobre (i zaimplementowane w Pythonie?) Metody? Metoda uzupełniania zerowego będzie działać, ale jest zbyt droga, aby była praktyczna.
keflavich
@kelavich: Czy przetestowałeś wszystkie metody interpolacji i znalazłeś niedopuszczalne uprzedzenie? Zalecanym podejściem jest kombinacja częściowego uzupełniania zerowego, a następnie interpolacji niskiego błędu. Nie znam żadnej innej metody, ale założę się, że zapewni ci to dużą dokładność.
endolith,
Tak, znalazłem niedopuszczalne odchylenie w interpolacji liniowej i drugiego rzędu. Próbowałem wypełniania zerowego FFT, ale wynik jest zdominowany przez dzwonienie o wysokiej częstotliwości ... czy jest szansa na dodanie przykładu wypełniania zerowego?
keflavich