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

Question

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 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 wavelength 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 (and are not periodic). I'd like to find 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))

Answer 1

I think using cross-correlation and interpolating the peak would work fine. As described in Is up-sampling prior to cross-correlation useless?, interpolating or upsampling before the cross-correlation doesn't actually get you any more information. The information about the sub-sample peak is contained in the samples around it. You just need to extract it with minimal error. I gathered some notes here.

The simplest method is quadratic/parabolic interpolation, which I have a Python example of here. It's supposedly exact if your spectrum is based on a Gaussian window, or if the peak happens to fall exactly on the midpoint between samples, but otherwise has some error. So in your case you probably want to use something better.

Here is a list of more complicated, but more accurate estimators. "Of the above methods, Quinn's second estimator has the least RMS error."

I don't know the math, but this paper says that their parabolic interpolation has theoretical accuracy of 5% of the width of an FFT bin.

Using FFT interpolation on the cross-correlation output does not have any bias error, so that's the best if you want really good accuracy. If you need to balance accuracy and computation speed, it's recommended to do some FFT interpolation and then follow it with one of the other estimators to get a "good enough" result.

This just uses the parabolic fit, but it outputs the right value for the offset if the noise is low:

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

The noise in your sample produces results that vary by more than a whole sample, so I reduced it. Fitting the curve using more of the points of the peak helps tighten the estimate somewhat, but I'm not sure if that's statistically valid, and it actually makes the estimate worse for the lower-noise situation.

With noise = 0.2 and 3-point fit, it gives values like 0.398 and 0.402 for offset = 0.4.

With noise = 2.0 and 13-point fit, it gives values like 0.156 and 0.595 for offset = 0.4.