# Directly compare subpixel shifts between two spectra — and get believable errors!

I have two spectra of the same astronomical object. The essential question is this: How can I calculate the relative shift between these spectra and get an accurate error estimate on that shift?

Some more details if you are still with me. Each spectrum will be an array with an x value (wavelength), y value (flux), and error. The wavelength shift is going to be sub-pixel. Assume that the pixels are regularly spaced and that there is only going to be a single shift applied to the entire spectrum. So the end answer will be something like: 0.35 +/- 0.25 pixels.

The two spectra are going to be a lot of featureless continuum punctuated by some rather complicated absorption features (dips) that do not model easily. I'd like to discover a method that directly compares the two spectra.

Everyone's first instinct is to do a cross-correlation, but with subpixel shifts, you're going to have to interpolate between the spectra (by smoothing first?) -- also, errors seem nasty to get right.

My current approach is to smooth the data by convolving with a gaussian kernel, then to spline the smoothed result, and compare the two splined spectra -- but I don't trust it (especially the errors).

Does anyone know of a way to do this properly?

Here is a short python program that will produce two toy spectra that are shifted by 0.4 pixels (written out in toy1.ascii and toy2.ascii) that you can play with. Even though this toy model uses a simple gaussian feature, assume that the actual data cannot be fit with a simple model.

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

-
What astronomical object are you observing? Stars? Galaxies? This can really help, because the SDSS team gets sub-pixel accuracy by first identifying the eigenbasis of the target galaxy via Principal Component Analysis. Check this draft by Aihara et al. (2011): icc.ub.edu/~liciaverde/dr8.pdf. –  Néstor May 11 '12 at 0:56