Vergleichen Sie Subpixelverschiebungen zwischen zwei Spektren direkt - und erhalten Sie glaubwürdige Fehler

9

Ich habe zwei Spektren desselben astronomischen Objekts. Die wesentliche Frage lautet: Wie kann ich die relative Verschiebung zwischen diesen Spektren berechnen und einen genauen Fehler bei dieser Verschiebung erhalten?

Noch ein paar Details, wenn du noch bei mir bist. Jedes Spektrum ist ein Array mit einem x-Wert (Wellenlänge), einem y-Wert (Fluss) und einem Fehler. Die Wellenlängenverschiebung wird ein Subpixel sein. Angenommen, die Pixel sind regelmäßig beabstandet und es wird nur eine einzige Wellenlängenverschiebung auf das gesamte Spektrum angewendet. Die Endantwort lautet also ungefähr: 0,35 +/- 0,25 Pixel.

Die beiden Spektren werden eine Menge merkwürdiges Kontinuum sein, das durch einige ziemlich komplizierte Absorptionsmerkmale (Einbrüche) unterbrochen wird, die sich nicht leicht modellieren lassen (und nicht periodisch sind). Ich möchte eine Methode finden, die die beiden Spektren direkt vergleicht.

Der erste Instinkt eines jeden besteht darin, eine Kreuzkorrelation durchzuführen, aber bei Subpixel-Verschiebungen müssen Sie zwischen den Spektren interpolieren (indem Sie zuerst glätten?) - auch Fehler scheinen böse zu sein, um richtig zu werden.

Mein aktueller Ansatz besteht darin, die Daten zu glätten, indem ich sie mit einem Gaußschen Kernel falte, dann das geglättete Ergebnis zu spalten und die beiden Spline-Spektren zu vergleichen - aber ich vertraue ihm nicht (insbesondere den Fehlern).

Kennt jemand einen Weg, dies richtig zu machen?

Hier ist ein kurzes Python-Programm, das zwei Spielzeugspektren erzeugt, die um 0,4 Pixel verschoben sind (geschrieben in toy1.ascii und toy2.ascii), mit denen Sie spielen können. Obwohl dieses Spielzeugmodell eine einfache Gauß-Funktion verwendet, wird davon ausgegangen, dass die tatsächlichen Daten nicht mit einem einfachen Modell übereinstimmen können.

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
quelle
Wenn ich das richtig verstehe, klingt das Problem ähnlich wie bei der Bildregistrierung, außer dass Sie nur eine lineare Subpixelverschiebung in einer Achse haben. Vielleicht versuchen Sie es mit Standard-Bildregistrierungstechniken wie der Phasenkorrelation?
Paul R
Wenn Sie eine reine Verzögerung in einem Signal haben (dh die Verschiebung des Wellenlängenparameters, von dem Sie sprechen), können Sie möglicherweise die Fourier-Transformationseigenschaft ausnutzen, die die Zeitverzögerung in einen linearen Phasenversatz im Frequenzbereich umwandelt. Dies könnte funktionieren, wenn die beiden Proben nicht durch unterschiedliche Messgeräusche oder Interferenzen verfälscht werden.
Jason R
1
Dieser Thread kann nützlich sein
Jim Clay
1
Haben Sie aktuelle Daten zum Testen? Der von Ihnen angegebene Rauschwert ist zu hoch, als dass die Kreuzkorrelation unterabtastgenau wäre. Dies ist, was es mit mehreren Rauschläufen 2.0 und Offset 0.7 (= 1000.7 auf der x-Achse des Plots) findet, zum Beispiel: i.stack.imgur.com/UK5JD.png
Endolith

Antworten:

5

Ich denke, Kreuzkorrelation und Interpolation des Peaks würden gut funktionieren. Wie unter Ist ein Up-Sampling vor einer Kreuzkorrelation nutzlos? Wenn Sie vor der Kreuzkorrelation interpolieren oder upsampling, erhalten Sie keine weiteren Informationen. Die Informationen über den Teilprobenpeak sind in den Proben um ihn herum enthalten. Sie müssen es nur mit minimalem Fehler extrahieren. Ich habe hier einige Notizen gesammelt .

Die einfachste Methode ist die quadratische / parabolische Interpolation, für die ich hier ein Python-Beispiel habe . Es ist angeblich genau, ob Ihr Spektrum auf einem Gaußschen Fenster basiert oder ob der Peak genau auf den Mittelpunkt zwischen den Samples fällt, ansonsten aber einen Fehler aufweist . In Ihrem Fall möchten Sie wahrscheinlich etwas Besseres verwenden.

Hier ist eine Liste komplizierterer, aber genauerer Schätzer. "Von den oben genannten Methoden weist Quinns zweiter Schätzer den geringsten RMS-Fehler auf."

Ich kenne die Mathematik nicht, aber dieses Papier sagt, dass ihre parabolische Interpolation eine theoretische Genauigkeit von 5% der Breite eines FFT-Behälters hat.

Die Verwendung der FFT-Interpolation am Kreuzkorrelationsausgang weist keinen Vorspannungsfehler auf. Dies ist also das Beste, wenn Sie eine wirklich gute Genauigkeit wünschen. Wenn Sie Genauigkeit und Rechengeschwindigkeit in Einklang bringen müssen, wird empfohlen, eine FFT-Interpolation durchzuführen und anschließend einen der anderen Schätzer zu verwenden, um ein "gut genug" -Ergebnis zu erhalten.

Dies verwendet nur die parabolische Anpassung, gibt jedoch den richtigen Wert für den Offset aus, wenn das Rauschen gering ist:

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

Geben Sie hier die Bildbeschreibung ein Geben Sie hier die Bildbeschreibung ein

Das Rauschen in Ihrer Probe führt zu Ergebnissen, die um mehr als eine ganze Probe variieren. Deshalb habe ich es reduziert. Das Anpassen der Kurve unter Verwendung von mehr Punkten des Peaks hilft dabei, die Schätzung etwas zu straffen, aber ich bin mir nicht sicher, ob dies statistisch gültig ist, und es verschlechtert die Schätzung tatsächlich für die Situation mit geringerem Rauschen.

Bei Rauschen = 0,2 und 3-Punkt-Anpassung ergeben sich Werte wie 0,398 und 0,402 für Offset = 0,4.

Bei Rauschen = 2,0 und 13-Punkt-Anpassung ergeben sich Werte wie 0,156 und 0,595 für Offset = 0,4.

Endolith
quelle
Ich versuche genau dieses Problem für die Bildregistrierung zu lösen. Ich brauche eine Subpixel-Genauigkeit (0,1 wäre wahrscheinlich gut genug), aber vor allem brauche ich keine Vorspannung, damit die Interpolationsmethoden nicht funktionieren. Gibt es dafür gute (und in Python implementierte?) Methoden? Die Zero-Padding-Methode funktioniert, ist jedoch zu teuer, um praktisch zu sein.
Keflavich
@kelavich: Haben Sie alle Interpolationsmethoden getestet und eine inakzeptable Verzerrung festgestellt? Die empfohlene Vorgehensweise ist eine Kombination aus einigem von einer fehlerarme Interpolation gefolgt Zero-Padding. Ich kenne keine andere Methode, aber ich wette, dies würde Ihnen viel Genauigkeit bieten.
Endolith
Ja, ich habe eine inakzeptable Verzerrung bei der Interpolation linearer und 2. Ordnung festgestellt. Ich habe FFT Zero-Padding ausprobiert, aber das Ergebnis wird von hochfrequentem Klingeln dominiert. Gibt es eine Möglichkeit, ein Zero-Padding-Beispiel hinzuzufügen?
Keflavich