next previous
Up: Random-error minimization during cross-correlation


Subsections

3 Cross-correlation peak centering

3.1 General remarks

We assume now we have obtained the cross-correlation function of an object spectrum and an intrinsically (i.e. apart from the noise) identical template spectrum and we are faced with the problem of finding its maximum position with the best possible precision. Since the width of even the narrowest stellar absorption lines is usually much larger than the positional error one is prepared to tolerate, spectra intended for radial-velocity measurements must be well oversampled. The actual degree of oversampling for the resulting cross-correlation function may vary widely, from typically 5 or 6 pixels within the central peak for late-type spectra to sometimes well over 100 for early-types. Even so, in most cases the sampling step is still larger than the tolerable error and one needs some form of interpolation between the sampling values of the correlation peak to obtain the required (sub-pixel) precision.

In practice one cannot expect automatically to obtain a precision corresponding to the lower bound for the random error discussed in the previous section. Strictly speaking, this would require an exact model with infinite resolution for the correlation peak. Nevertheless, if the spectra contain nothing but well-separated unbroadened metal lines then a Gaussian model for the cross-correlation function may come close to this ideal situation. Most often, however, line blending and the mixture of H, He and rotationally broadened metal lines in early-type stars give the correlation peak a complicated structure for which no well-matching analytical model exists (see Fig. 1) and then, depending on the technique used to interpolate between the sampling values of the correlation peak, the expected random error at the sub-pixel level may become much larger than the lower bound. Tonry & Davis (1979) already used a combination of two kinds of low-pass filtering (general purpose Fourier filtering and the fit of a parabola to the full width of the correlation peak) for the reduction of noise effects in the case of late-type spectra. Below we argue that this is in fact essential in the case of early types and we study how such a combination should be tuned.

Qualitatively one can argue as follows. The substructure of the correlation peak, away from the centre, is due to the overlap of non-corresponding features in the spectra so it is irrelevant for the localization of the correlation maximum. If only low-pass Fourier filtering is applied and the correlation maximum is obtained correspondingly by Fourier-interpolation of the filtered peak, then that substructure will inevitably influence the position. On the other hand, if only a model-fit were applied, then again the substructure would influence the position of the maximum unless the model were a perfect match. So we need to eliminate the influence of the spurious substructure as far as possible. Obviously this can be done by fitting a model only to the very central part of the correlation peak. But then only few samples are used and the noise on these could strongly increase the random error. Hence the need for low-pass filtering prior to the model fit. The combined procedure of course requires two "smoothing parameters" to be chosen (fit-interval and Fourier-bandwidth); we must also address the question of which model to choose for the fit.

3.2 Fit-points and fit-functions

In order to investigate how the random error behaves as a function of the number of fit-points, and when different polynomial models for the correlation peak are used, Monte-Carlo experiments were conducted as follows: an object and a template spectrum are created by twice adding photon and read-out noise with a given total rms in the continuum (corresponding to a given continuum S/N) to a synthetic spectrum, which may have been rotationally broadened. In all our experiments, we fix the continuum S/N for the template spectrum at the fairly high value of 200, and fix the read-out noise for both object and template at 10 e-. For details on the synthetic spectra we refer to Vrancken et al. (1997). All spectra have a 1 pixel resolution of 2810-6 in ln($\lambda$) ($\pm$ 0.13 Å/${\rm p}$ at 4500 Å) but as long as the spectra are well sampled, all our conclusions are independent of resolution.

The spectra are then cross-correlated using the CORSPEC package (updated from Verschueren 1991) and the maximum position of the cross-correlation peak is fitted using a number of fit-points ranging from 5 to twice its FWHM in steps of 2 pixels, and using 3 different fit-functions: a 2$^{\rm nd}$ degree, a symmetric 4$^{\rm th}$ degree, and an ordinary 4th degree polynomial. Obviously, with our experimental set-up, a symmetric fit-function is appropriate. In the case of observed spectra, however, small systematic mismatch between object and template may give reason to use a fit-function allowing for asymmetry. Therefore, the ordinary 4$^{\rm th}$ degree polynomial is included. Note that a 3$^{\rm rd}$ degree polynomial is unsuitable because its uneven highest order term causes unacceptably high sensitivity for asymmetries in the random noise. The use of a Gaussian fit-function is also less appropriate since its parameters are not well constrained by information from the mere top of a peaked function (unless one constrains the background term, but this value is not known a priori for early-type spectra). Discretisation errors during centering are minimized according to the prescription given in David & Verschueren (1995). This process is repeated n = 2000 times with independent random noise, and the rms of the fitted central positions, representing the random error on a single measurement of the shift, is computed for each value of the number of fit-points and for each fit-function. The relative error on such a rms is $\sim$ $1/\sqrt{2n} = 0.016$.

  
\begin{figure}
\includegraphics [width=8.8cm]{fig_rmsfct.eps}\end{figure} Figure 2: Random error on a cross-correlation shift as a function of the length of the fit-interval used, and for the three polynomial fit-functions indicated; it is derived from a Monte-Carlo experiment with 2000 runs. The lengths of the fit-interval are normalised to the FWHM of the cross-correlation peak (about 40 pixels in this case); they correspond to a number of fit-points starting from 5 pixels with increments of 2. A B1V synthetic spectrum in the region 4382-4759 Å was used with the indicated rotational velocity and object continuum S/N. The dashed line indicates the level of the lower bound from Sect. 2

Figure 2 shows an example and compares with the lower bound for the random error discussed in Sect. 2. The details of the curves depend on the spectral characteristics and on the S/N of object and template, but the following important properties are typical:

In order to minimize the random error, one would thus like to fit the cross-correlation function using the optimal number of fit-points for all spectra. However, in practice the value of $n_{\rm opt}$ is not known a priori since it depends intricately on the spectral characteristics. Moreover, it is preferable to use as few fit-points as possible for reasons mentioned in Sect. 3.1. Therefore, we investigated whether high-frequency Fourier filtering could provide a solution for both constraints.

3.3 High-frequency filtering

  
\begin{figure}
\includegraphics [width=8.8cm]{fig_ccph.eps}\end{figure} Figure 3: Phase of the Fourier transform of the cross-correlation function versus wavenumber for an arbitrary couple of spectra used for Fig. 2. The x-axis is truncated for clarity

The optimal cut-off frequency (or wavenumber) $f_{\rm c}$ in Fourier space may be derived from a plot of the phase of the Fourier transform of the cross-correlation function versus its frequency. Figure 3 shows, for 2 arbitrary spectra from the set used for Fig. 2, a linear behaviour for frequencies dominated by intrinsic spectral information, a region where the phase is randomized for noise-dominated frequencies, and an intermediate region. $f_{\rm c}$ $\cong$ 250 seems a safe choice for removing only fully noise-dominated information in this case. $f_{\rm c}$ can conveniently be expressed as N/i with N the total number of pixels in the spectrum. The "scaled filter value" i = $N/f_{\rm c}$ defined in this way is independent of the exact length of the spectrum; it is also independent of the sampling rate (i.e. the spectral resolution) provided that the spectra are well sampled and that the S/N per pixel remains the same when changing the resolution. The optimal filter i thus only depends on the type of spectrum (i.e. intrinsic characteristics of the lines and rotational velocity) and on its S/N. For the spectrum used in Figs. 2 and 3 (for which N=4096), $f_{\rm c}$ $\cong$ 250 thus corresponds to N/16 filtering (i=16). Note that one can always choose i $\geq$ 4, with N/4 filtering being optimal for spectral lines with a FWHM equal to the 2 pixel spectral resolution. Since our filter choice is always very conservative, i.e. only fully noise-dominated frequencies are removed, a simple rectangular filter can be used.

  
\begin{figure}
\includegraphics [width=8.8cm]{fig_rmsfilt.eps}\end{figure} Figure 4: Idem as Fig. 2, but now showing also the effect of 3 different high-frequency filters. The results are shown only for the quadratic fit-function. The small rise in minimum rms for the last two filters is insignificant and is an occasional effect that does not necessarily occur for other spectral regions or parameters

Figure 4 then shows the influence of successively increased filtering on the random error from an identical Monte-Carlo experiment as in Fig. 2, and for the quadratic fit-function. The behaviour is again typical for all spectra: filtering up to the optimal cut-off frequency does not significantly change the minimum random error, but allows to attain it in practice for all fit-intervals smaller than the optimal one. This allows one to fit cross-correlation peaks close to their centre with a fixed number of fit-points, and ensures an optimal precision at the same time. Further tests have shown that filtering significantly beyond the cut-off frequency increases the minimum random error accordingly due to a degradation of the intrinsic spectral features. Note that, although high-frequency filtering apparently increases the S/N of the spectra, it can not improve the actual optimal precision attainable from those spectra: indeed, the latter is obviously limited by those random-noise frequencies that are also present in the intrinsic spectral information itself.

  
\begin{figure}
\includegraphics [width=8.8cm]{fig_rmsphfr.eps}\end{figure} Figure 5: a) Monte-Carlo rms of the phase of the Fourier transform of the cross-correlation function versus scaled filter value i = N/f. The same input parameters were used as for Fig. 2. The x-axis is truncated towards higher values for clarity. The range i = 50 to 2 corresponds to wavenumber = 82 to 2048 in Fig. 3. The full line corresponds to full randomization of the phase. b) Idem but now for a S/N = 200

Since it is unpractical and often ambiguous to determine $f_{\rm c}$ or i from plots such as Fig. 3 for each given pair of spectra, we investigated its dependence on spectral characteristics and on S/N. This was done in a set of Monte-Carlo experiments (1000 runs) with different synthetic spectra (only metal lines in a B1V spectrum, only H lines in a B8V spectrum, and each for vsini = 5, 50, 100, 150, 200 and 300 km s-1) and with different values of the S/N (50 and 200). In each case, the Monte-Carlo rms of the phase of the Fourier transform of the cross-correlation function was computed as a function of Fourier frequency f, and compared to the value of 1.81 rad which corresponds to a random phase between $-\pi$ and $+\pi$.

Figure 5a shows an example of such an experiment with the same input parameters as used in Fig. 2; the plot is already made as a function of N/f = i which allows direct reading of the scaled filter value to be used. For high values of i (small frequencies), spectral information still dominates and produces low rms values (note also the smaller rms around frequencies related to the rotational broadening function, see Gray 1976). A rms of about 1.8 rad is reached below a cut-off value of i $\cong$ 15 which corresponds well to the value of 16 derived earlier from a single spectrum (Fig. 3) and to the value empirically needed to flatten the curve in Fig. 4.

In order to investigate the influence of the spectrum S/N, the same Monte-Carlo experiment was repeated with a S/N = 200. Figure 5b shows that a cut-off value of i $\cong$ 10 seems more appropriate now due to the fact that intrinsic spectral information now beats the noise up to higher frequencies. From the whole set of experiments performed on other spectra, we can indeed generally conclude that a more pronounced filtering may be applied in case of lower S/N.

  
\begin{figure}
\includegraphics [width=8.8cm]{fig_fcvsini.eps}\end{figure} Figure 6: Optimal scaled filter value i = $N/f_{\rm c}$ for high-frequency filtering as a function of rotational velocity, for spectra containing metal respectively H lines as indicated. In each case, the lower value refers to an object spectrum S/N = 200 and the higher one to a S/N = 50

The influence of rotational velocity (directly influencing line-width) was examined by performing all Monte-Carlo experiments for the 6 values of vsini mentioned above. Figure 6 summarizes the results for the two different values of the S/N and for metal and H-lines separately. For a given spectrum sampled in N pixels, the adequate value of i may be read from these curves and can be used to compute the optimal cut-off frequency $f_{\rm c}$ = N/i for high-frequency filtering, independent of spectral resolution (for well resolved spectra).

Finally, we verified that the conclusions found in this section are valid independently of the fit-function used. However, we give preference to keeping the number of free parameters as small as possible, especially since this yields the smallest increase in random error in case the chosen degree of filtering was not optimal (see Fig. 2). In addition, for spectra without systematic mismatch, the fit-function needs to be symmetric. Our proposed fitting method is therefore a parabolic fit to the highest 5 pixels of the cross-correlation function, after applying a degree of high-frequency filtering determined from Fig. 6.


next previous
Up: Random-error minimization during cross-correlation

Copyright The European Southern Observatory (ESO)