Непосредственно сравнивайте субпиксельные сдвиги между двумя спектрами - и получите правдоподобные ошибки

9

У меня два спектра одного и того же астрономического объекта. Основной вопрос заключается в следующем: как я могу рассчитать относительный сдвиг между этими спектрами и получить точную ошибку в этом сдвиге?

Еще некоторые подробности, если вы все еще со мной. Каждый спектр будет массивом со значением x (длина волны), значением y (поток) и ошибкой. Сдвиг длины волны будет субпиксельным. Предположим, что пиксели расположены на регулярном расстоянии, и будет иметь место только одно смещение длины волны, примененное ко всему спектру. Таким образом, конечный ответ будет примерно таким: 0,35 +/- 0,25 пикселей.

Два спектра будут представлять собой множество безликих континуумов, перемежающихся некоторыми довольно сложными абсорбционными характеристиками (провалами), которые не моделируются легко (и не являются периодическими). Я хотел бы найти метод, который напрямую сравнивает два спектра.

Первым инстинктом каждого является взаимная корреляция, но с помощью субпиксельных сдвигов вам придется интерполировать спектры (сначала сглаживая?), А также ошибки кажутся неприятными, чтобы исправляться.

Мой текущий подход состоит в том, чтобы сгладить данные путем свертки с гауссовым ядром, затем сплайнировать сглаженный результат и сравнить два сплайновых спектра - но я не доверяю этому (особенно ошибкам).

Кто-нибудь знает способ сделать это правильно?

Вот короткая программа на python, которая создаст два игрушечных спектра, сдвинутых на 0,4 пикселя (записанных в toy1.ascii и toy2.ascii), с которыми вы можете играть. Даже если в этой игрушечной модели используется простая гауссовская функция, предположим, что фактические данные не могут соответствовать простой модели.

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
источник
Если я правильно понимаю, проблема звучит похоже на регистрацию изображения, за исключением того, что у вас просто линейный сдвиг субпикселя по одной оси. Может быть, попробовать стандартные методы регистрации изображений, такие как фазовая корреляция?
Пол Р
Если у вас есть чистая задержка в одном сигнале (т.е. сдвиг в параметре длины волны, о котором вы говорите), вы можете использовать свойство преобразования Фурье, которое превращает временную задержку в линейный сдвиг фазы в частотной области. Это может сработать, если две выборки не повреждены из-за различного измерительного шума или помех.
Джейсон Р
1
Эта ветка может быть полезна
Джим Клэй,
1
У вас есть реальные данные для тестирования? Значение шума, которое вы указали, слишком велико, чтобы кросс-корреляция была точной для подвыборки. Это то , что он находит с несколькими пробегами шума 2.0 и 0.7 (смещение = 1000.7 на оси х сюжетов в), например: i.stack.imgur.com/UK5JD.png
эндолиты

Ответы:

5

Я думаю, что использование взаимной корреляции и интерполяции пика будет работать нормально. Как описано в разделе «Бесполезно ли повышение частоты дискретизации до взаимной корреляции?». интерполяция или повышающая дискретизация до взаимной корреляции фактически не дает вам больше информации. Информация о пике подвыборки содержится в образцах вокруг него. Вам просто нужно извлечь его с минимальной ошибкой. Я собрал несколько заметок здесь .

Самый простой метод - это квадратичная / параболическая интерполяция, пример которой приведен здесь на Python . Предположительно точно, если ваш спектр основан на гауссовском окне , или если пик оказывается точно в средней точке между выборками, но в противном случае есть некоторая ошибка . Так что в вашем случае вы, вероятно, хотите использовать что-то лучшее.

Вот список более сложных, но более точных оценок. «Из приведенных выше методов вторая оценка Квинна имеет наименьшую среднеквадратичную ошибку».

Я не знаю математики, но в этой статье говорится, что их параболическая интерполяция имеет теоретическую точность 5% ширины ячейки БПФ.

Использование FFT-интерполяции на выходе кросс-корреляции не имеет ошибки смещения , так что это лучше всего, если вы хотите действительно хорошую точность. Если вам необходимо сбалансировать точность и скорость вычислений, рекомендуется выполнить некоторую FFT-интерполяцию, а затем выполнить ее с одним из других оценщиков, чтобы получить «достаточно хороший» результат.

Это просто использует параболическое соответствие, но выводит правильное значение для смещения, если шум низкий:

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

введите описание изображения здесь введите описание изображения здесь

Шум в вашем образце дает результаты, которые варьируются более чем на целый образец, поэтому я уменьшил его. Подгонка кривой с использованием большего количества точек пика помогает несколько сузить оценку, но я не уверен, является ли она статистически достоверной, и это фактически ухудшает оценку для ситуации с низким уровнем шума.

С шумом = 0,2 и 3-точечной посадкой он дает такие значения, как 0,398 и 0,402 для смещения = 0,4.

С шумом = 2,0 и 13-точечной посадкой он дает значения, такие как 0,156 и 0,595 для смещения = 0,4.

эндолиты
источник
Я пытаюсь решить именно эту проблему для регистрации изображений. Мне нужна субпиксельная точность (0.1, вероятно, будет достаточно), но самое главное, не нужно смещения, поэтому методы интерполяции не работают. Есть ли хорошие (и реализованные в python?) Методы для этого? Метод заполнения нулями будет работать, но он слишком дорог, чтобы быть практичным.
Кефлавич
@kelavich: Вы проверили все методы интерполяции и нашли недопустимое отклонение? Рекомендуемый подход представляет собой комбинацию некоторого дополнения нулями с последующей интерполяцией с низкой ошибкой. Я не знаю ни одного другого метода, но держу пари, что это даст вам достаточно точности.
эндолит
Да, я обнаружил недопустимое смещение в линейной интерполяции и интерполяции 2-го порядка. Я пробовал FFT с нулевым заполнением, но в результате преобладает высокочастотный звонок ... есть ли шанс добавить пример с нулевым заполнением?
Кефлавич