Secara langsung membandingkan perpindahan subpiksel antara dua spektrum - dan dapatkan kesalahan yang dapat dipercaya

Saya memiliki dua spektrum objek astronomi yang sama. Pertanyaan penting adalah ini: Bagaimana saya bisa menghitung pergeseran relatif antara spektrum ini dan mendapatkan kesalahan akurat pada pergeseran itu?

Beberapa detail lagi jika Anda masih bersama saya. Setiap spektrum akan berupa array dengan nilai x (panjang gelombang), nilai y (fluks), dan kesalahan. Pergeseran panjang gelombang akan menjadi sub-pixel. Asumsikan bahwa piksel ditempatkan secara teratur dan hanya akan ada satu pergeseran panjang gelombang tunggal yang diterapkan ke seluruh spektrum. Jadi jawaban akhirnya akan seperti: 0,35 +/- 0,25 piksel.

Dua spektrum akan menjadi banyak kontinum tanpa ciri diselingi oleh beberapa fitur penyerapan yang agak rumit (dips) yang tidak mudah dimodelkan (dan tidak periodik). Saya ingin menemukan metode yang secara langsung membandingkan dua spektrum.

Insting pertama setiap orang adalah melakukan korelasi silang, tetapi dengan pergeseran subpixel, Anda harus melakukan interpolasi di antara spektrum (dengan menghaluskan terlebih dahulu?) - juga, kesalahan tampaknya tidak tepat untuk dilakukan dengan benar.

Pendekatan saya saat ini adalah memuluskan data dengan berbelit-belit dengan kernel gaussian, kemudian untuk spline hasil yang diperhalus, dan membandingkan dua spektrum splined - tapi saya tidak percaya itu (terutama kesalahan).

Adakah yang tahu cara melakukan ini dengan benar?

Berikut adalah program python pendek yang akan menghasilkan dua spektrum mainan yang digeser 0,4 piksel (ditulis dalam toy1.ascii dan toy2.ascii) yang dapat Anda mainkan. Meskipun model mainan ini menggunakan fitur gaussian sederhana, asumsikan bahwa data aktual tidak dapat cocok dengan model sederhana.

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

Saya pikir menggunakan korelasi silang dan interpolasi puncak akan bekerja dengan baik. Seperti yang dijelaskan dalam Apakah pengambilan sampel sebelum korelasi silang tidak berguna? , interpolasi atau peningkatan sebelum korelasi silang sebenarnya tidak membuat Anda mendapatkan informasi lebih lanjut. Informasi tentang puncak sub-sampel terkandung dalam sampel di sekitarnya. Anda hanya perlu mengekstraknya dengan kesalahan minimal. Saya mengumpulkan beberapa catatan di sini .

Metode paling sederhana adalah kuadratik / interpolasi parabola, yang saya punya contoh Python di sini . Seharusnya tepat jika spektrum Anda didasarkan pada jendela Gaussian , atau jika puncaknya jatuh tepat di titik tengah di antara sampel, tetapi sebaliknya memiliki beberapa kesalahan . Jadi, dalam kasus Anda, Anda mungkin ingin menggunakan sesuatu yang lebih baik.

Berikut adalah daftar penaksir yang lebih rumit, tetapi lebih akurat. "Dari metode di atas, penaksir kedua Quinn memiliki kesalahan RMS paling sedikit."

Saya tidak tahu matematika, tetapi makalah ini mengatakan bahwa interpolasi parabola mereka memiliki akurasi teoritis 5% dari lebar bin FFT.

Menggunakan interpolasi FFT pada output korelasi silang tidak memiliki kesalahan bias , jadi itu yang terbaik jika Anda ingin akurasi yang sangat baik. Jika Anda perlu menyeimbangkan akurasi dan kecepatan perhitungan, disarankan untuk melakukan interpolasi FFT dan kemudian ikuti dengan salah satu penaksir lainnya untuk mendapatkan hasil "cukup baik".

Ini hanya menggunakan kesesuaian parabola, tetapi menghasilkan nilai yang tepat untuk offset jika noise rendah:

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

Kebisingan dalam sampel Anda menghasilkan hasil yang bervariasi lebih dari keseluruhan sampel, jadi saya menguranginya. Menyesuaikan kurva dengan menggunakan lebih banyak titik puncak agak membantu memperketat perkiraan, tapi saya tidak yakin apakah itu valid secara statistik, dan itu benar-benar membuat perkiraan lebih buruk untuk situasi dengan kebisingan yang lebih rendah.

Dengan noise = 0,2 dan 3-point fit, ia memberikan nilai seperti 0,398 dan 0,402 untuk offset = 0,4.

Dengan noise = 2.0 dan fit 13-point, ia memberikan nilai seperti 0,156 dan 0,595 untuk offset = 0,4.