Compare diretamente as alterações de subpixel entre dois espectros - e obtenha erros confiáveis

9

Eu tenho dois espectros do mesmo objeto astronômico. A questão essencial é a seguinte: como posso calcular a mudança relativa entre esses espectros e obter um erro preciso nessa mudança?

Mais alguns detalhes se você ainda estiver comigo. Cada espectro será uma matriz com um valor x (comprimento de onda), valor y (fluxo) e erro. O deslocamento do comprimento de onda será sub-pixel. Suponha que os pixels sejam espaçados regularmente e que haverá apenas uma única mudança de comprimento de onda aplicada a todo o espectro. Portanto, a resposta final será algo como: 0,35 +/- 0,25 pixels.

Os dois espectros terão um continuum sem característica pontuado por alguns recursos de absorção bastante complicados (quedas) que não modelam facilmente (e não são periódicos). Eu gostaria de encontrar um método que compare diretamente os dois espectros.

O primeiro instinto de todo mundo é fazer uma correlação cruzada, mas com as mudanças de subpixel, você terá que interpolar entre os espectros (suavizando primeiro?) - além disso, os erros parecem desagradáveis ​​para acertar.

Minha abordagem atual é suavizar os dados convolvendo com um kernel gaussiano, depois estragar o resultado suavizado e comparar os dois espectros estriados - mas não confio neles (especialmente nos erros).

Alguém sabe uma maneira de fazer isso corretamente?

Aqui está um pequeno programa python que produzirá dois espectros de brinquedos com deslocamento de 0,4 pixels (gravados em toy1.ascii e toy2.ascii) com os quais você pode brincar. Mesmo que esse modelo de brinquedo use um recurso gaussiano simples, suponha que os dados reais não possam ser ajustados a um modelo simples.

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
fonte
Se bem entendi, o problema parece semelhante ao registro de imagens, exceto que você tem apenas uma mudança linear de sub-pixels em um eixo. Talvez tente técnicas de registro de imagem padrão, como correlação de fase?
Paul R
Se você tiver um atraso puro em um sinal (ou seja, a mudança no parâmetro de comprimento de onda de que fala), poderá explorar a propriedade de transformação de Fourier que transforma o atraso de tempo em um deslocamento de fase linear no domínio da frequência. Isso pode funcionar se as duas amostras não forem corrompidas por diferentes ruídos ou interferências de medição.
Jason R
11
Esta discussão pode ser útil- dsp.stackexchange.com/questions/2321/...
Jim Clay
11
Você tem dados reais para testar? O valor do ruído que você forneceu é muito alto para que a correlação cruzada seja precisa na subamostra. Isto é o que ele encontra com várias corridas de ruído 2.0 e deslocamento de 0,7 (= 1000,7 sobre o enredo eixo X), por exemplo: i.stack.imgur.com/UK5JD.png
endolith

Respostas:

5

Eu acho que usar correlação cruzada e interpolar o pico funcionaria bem. Conforme descrito em A amostragem prévia antes da correlação cruzada é inútil? , interpolar ou aumentar a amostragem antes da correlação cruzada não fornecer mais informações. As informações sobre o pico da subamostra estão contidas nas amostras ao seu redor. Você só precisa extraí-lo com um erro mínimo. Reuni algumas anotações aqui .

O método mais simples é a interpolação quadrática / parabólica, da qual tenho um exemplo em Python aqui . É supostamente exato se seu espectro é baseado em uma janela gaussiana ou se o pico cai exatamente no ponto médio entre as amostras, mas, de outra forma, apresenta algum erro . Portanto, no seu caso, você provavelmente deseja usar algo melhor.

Aqui está uma lista de estimadores mais complicados, mas mais precisos. "Dos métodos acima, o segundo estimador de Quinn tem o menor erro de RMS."

Não sei matemática, mas este artigo diz que a interpolação parabólica tem precisão teórica de 5% da largura de uma caixa de FFT.

O uso da interpolação FFT na saída de correlação cruzada não apresenta nenhum erro de polarização , portanto é melhor se você deseja uma precisão realmente boa. Se você precisar equilibrar precisão e velocidade de computação, é recomendável fazer uma interpolação por FFT e segui-la com um dos outros estimadores para obter um resultado "bom o suficiente".

Isso apenas usa o ajuste parabólico, mas gera o valor certo para o deslocamento se o ruído for baixo:

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

insira a descrição da imagem aqui insira a descrição da imagem aqui

O ruído na sua amostra produz resultados que variam em mais do que uma amostra inteira, então eu o reduzi. Ajustar a curva usando mais pontos do pico ajuda a apertar um pouco a estimativa, mas não tenho certeza se isso é estatisticamente válido, e na verdade piora a estimativa para a situação de baixo ruído.

Com ruído = ajuste de 0,2 e 3 pontos, fornece valores como 0,398 e 0,402 para deslocamento = 0,4.

Com ruído = ajuste de 2,0 e 13 pontos, fornece valores como 0,156 e 0,595 para deslocamento = 0,4.

endólito
fonte
Estou tentando resolver esse problema exato para o registro de imagens. Preciso de precisão sub-pixel (0,1 provavelmente seria bom o suficiente), mas o mais importante não precisa de viés, para que os métodos de interpolação não funcionem. Existem métodos bons (e implementados em python?) Para isso? O método de preenchimento zero funcionará, mas é muito caro para ser prático.
Keflavich
@kelavich: Você já testou todos os métodos de interpolação e encontrou um viés inaceitável? A abordagem recomendada é uma combinação de alguns padding-zero seguido por uma interpolação baixa taxa de erros. Não conheço outro método, mas aposto que isso forneceria muita precisão.
endolith 21/08/12
Sim, eu encontrei viés inaceitável na interpolação linear e de 2ª ordem. Eu tentei o preenchimento zero da FFT, mas o resultado é dominado pelo toque de alta frequência ... alguma chance de você adicionar um exemplo de preenchimento zero?
Keflavich