Comparez directement les décalages de sous-pixels entre deux spectres - et obtenez des erreurs crédibles

9

J'ai deux spectres du même objet astronomique. La question essentielle est la suivante: comment puis-je calculer le décalage relatif entre ces spectres et obtenir une erreur précise sur ce décalage?

Encore plus de détails si vous êtes toujours avec moi. Chaque spectre sera un tableau avec une valeur x (longueur d'onde), une valeur y (flux) et une erreur. Le décalage de longueur d'onde va être sous-pixel. Supposons que les pixels soient régulièrement espacés et qu'il n'y aura qu'un seul décalage de longueur d'onde appliqué à l'ensemble du spectre. Ainsi, la réponse finale sera quelque chose comme: 0,35 +/- 0,25 pixels.

Les deux spectres vont être constitués d'un continuum sans particularité ponctué de quelques caractéristiques d'absorption plutôt compliquées (creux) qui ne se modélisent pas facilement (et ne sont pas périodiques). Je voudrais trouver une méthode qui compare directement les deux spectres.

Le premier réflexe de tout le monde est de faire une corrélation croisée, mais avec des changements de sous-pixels, vous devrez interpoler entre les spectres (en lissant d'abord?) - De plus, les erreurs semblent désagréables pour corriger.

Mon approche actuelle consiste à lisser les données en convoluant avec un noyau gaussien, puis à spliner le résultat lissé et à comparer les deux spectres splinés - mais je ne lui fais pas confiance (en particulier les erreurs).

Quelqu'un connaît-il un moyen de le faire correctement?

Voici un petit programme en python qui produira deux spectres de jouets décalés de 0,4 pixel (écrits en toy1.ascii et toy2.ascii) avec lesquels vous pourrez jouer. Même si ce modèle de jouet utilise une fonction gaussienne simple, supposons que les données réelles ne peuvent pas être adaptées à un modèle 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
la source
Si je comprends bien, le problème ressemble à l'enregistrement d'image, sauf que vous avez juste un décalage linéaire de sous-pixels sur un axe. Essayez peut-être des techniques d'enregistrement d'images standard telles que la corrélation de phase?
Paul R
Si vous avez un retard pur dans un signal (c'est-à-dire le décalage du paramètre de longueur d'onde dont vous parlez), vous pourriez être en mesure d'exploiter la propriété de transformation de Fourier qui transforme le retard en un décalage de phase linéaire dans le domaine fréquentiel. Cela pourrait fonctionner si les deux échantillons ne sont pas corrompus par des bruits de mesure ou des interférences différents.
Jason R
1
Ce fil peut être utile
Jim Clay
1
Avez-vous des données réelles pour tester? La valeur de bruit que vous avez donnée est trop élevée pour que la corrélation croisée soit précise pour le sous-échantillon. C'est ce qu'il trouve avec plusieurs séquences de bruit 2.0 et un décalage de 0,7 (= 1000,7 sur l'axe x du tracé), par exemple: i.stack.imgur.com/UK5JD.png
endolith

Réponses:

5

Je pense que l'utilisation de la corrélation croisée et de l'interpolation du pic fonctionnerait bien. Comme décrit dans Le suréchantillonnage avant la corrélation croisée est-il inutile? , l'interpolation ou le suréchantillonnage avant la corrélation croisée ne vous permet pas d'obtenir plus d'informations. Les informations sur le pic du sous-échantillon sont contenues dans les échantillons qui l'entourent. Vous avez juste besoin de l'extraire avec une erreur minimale. J'ai rassemblé quelques notes ici .

La méthode la plus simple est l'interpolation quadratique / parabolique, dont j'ai un exemple Python ici . Il est censé être exact si votre spectre est basé sur une fenêtre gaussienne , ou si le pic tombe exactement au milieu des échantillons, mais sinon il a une erreur . Donc, dans votre cas, vous voudrez probablement utiliser quelque chose de mieux.

Voici une liste d'estimateurs plus compliqués mais plus précis. "Parmi les méthodes ci-dessus, le deuxième estimateur de Quinn a le moins d'erreur RMS."

Je ne connais pas les mathématiques, mais cet article dit que leur interpolation parabolique a une précision théorique de 5% de la largeur d'un bac FFT.

L'utilisation de l'interpolation FFT sur la sortie de corrélation croisée n'a pas d'erreur de biais , c'est donc le meilleur si vous voulez une très bonne précision. Si vous avez besoin d'équilibrer la précision et la vitesse de calcul, il est recommandé d'effectuer une interpolation FFT, puis de la suivre avec l'un des autres estimateurs pour obtenir un résultat "assez bon".

Cela utilise simplement l'ajustement parabolique, mais il génère la bonne valeur pour le décalage si le bruit est faible:

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

entrez la description de l'image ici entrez la description de l'image ici

Le bruit dans votre échantillon produit des résultats qui varient de plus d'un échantillon entier, donc je l'ai réduit. L'ajustement de la courbe en utilisant plus de points du pic permet de resserrer quelque peu l'estimation, mais je ne sais pas si c'est statistiquement valide, et cela aggrave en fait l'estimation pour la situation à faible bruit.

Avec un bruit = 0,2 et un ajustement en 3 points, il donne des valeurs comme 0,398 et 0,402 pour un décalage = 0,4.

Avec un bruit = 2,0 et un ajustement à 13 points, il donne des valeurs comme 0,156 et 0,595 pour un décalage = 0,4.

endolith
la source
J'essaie de résoudre ce problème exact pour l'enregistrement d'images. J'ai besoin d'une précision inférieure au pixel (0,1 serait probablement suffisant) mais surtout, je n'ai pas besoin de biais, donc les méthodes d'interpolation ne fonctionnent pas. Existe-t-il de bonnes méthodes (et implémentées en python?) Pour cela? La méthode de remplissage nul fonctionnera, mais elle est trop chère pour être pratique.
keflavich
@kelavich: Avez-vous testé toutes les méthodes d'interpolation et trouvé un biais inacceptable? L'approche recommandée est une combinaison de certains padding zéro suivi d'une interpolation à faible erreur. Je ne connais aucune autre méthode, mais je parie que cela vous fournirait beaucoup de précision.
endolith
Oui, j'ai trouvé un biais inacceptable dans l'interpolation linéaire et de second ordre. J'ai essayé le remplissage zéro FFT, mais le résultat est dominé par la sonnerie haute fréquence ... avez-vous des chances d'ajouter un exemple de remplissage zéro?
keflavich