Doppler shifting a spectrum¶
- PyAstronomy.pyasl.dopplerShift(wvl, flux, v, edgeHandling=None, fillValue=None, vlim=0.05, err=None)¶
Doppler shift a given spectrum.
An algorithm to apply a Doppler shift to a spectrum. The idea here is to obtain a shifted spectrum without loosing the wavelength axis. Therefore, this function, first, calculates the shifted wavelength axis and, second, obtains the new, shifted flux array at the old, unshifted wavelength points by linearly interpolating. No relativistic effects are considered.
Due to the shift, some bins at the edge of the spectrum cannot be interpolated, because they are outside the given input range. The default behavior of this function is to return numpy.NAN values at those points. One can, however, specify the edgeHandling parameter to choose a different handling of these points.
If “firstlast” is specified for edgeHandling, the out-of-range points at the red or blue edge of the spectrum will be filled using the first (at the blue edge) and last (at the red edge) valid point in the shifted, i.e., the interpolated, spectrum.
If “fillValue” is chosen for edge handling, the points under consideration will be filled with the value given through the fillValue keyword.
Warning
Shifting a spectrum using linear interpolation has an effect on the noise of the spectrum. No treatment of such effects is implemented in this function.
- Parameters
- wvlarray
Input wavelengths in A.
- fluxarray
Input flux.
- vfloat
Doppler shift in km/s
- edgeHandlingstring, {“fillValue”, “firstlast”}, optional
The method used to handle the edges of the output spectrum.
- fillValuefloat, optional
If the “fillValue” is specified as edge handling method, the value used to fill the edges of the output spectrum.
- vlimfloat, optional
Maximal fraction of the speed of light allowed for Doppler shift, v. Default is 0.05.
- errarray, optional
The uncertainties of the data points in the spectrum. If given, updated uncertainties are obtained by linear interpolation of the uncertainties.
- Returns
- nfluxarray
The shifted flux array at the old input locations.
- wlprimearray
The shifted wavelength axis.
- errarray, optional
Uncertainties of the shifted spectrum.
Example: Shifting¶
from __future__ import print_function, division
from PyAstronomy import pyasl
import matplotlib.pylab as plt
import numpy as np
# Create a "spectrum" with 0.01 A binning ...
wvl = np.linspace(6000., 6100., 10000)
# ... a gradient in the continuum ...
flux = np.ones(len(wvl)) + (wvl/wvl.min())*0.05
# ... and a Gaussian absorption line
flux -= np.exp(-(wvl-6050.)**2/(2.*0.5**2))*0.05
# Shift that spectrum redward by 20 km/s using
# "firstlast" as edge handling method.
nflux1, wlprime1 = pyasl.dopplerShift(wvl, flux, 20., edgeHandling="firstlast")
# Shift the red-shifted spectrum blueward by 20 km/s, i.e.,
# back on the initial spectrum.
nflux2, wlprime = pyasl.dopplerShift(wvl, nflux1, -20.,
edgeHandling="fillValue", fillValue=1.0)
# Check the maximum difference in the central part
indi = np.arange(len(flux)-200) + 100
print("Maximal difference (without outer 100 bins): ",
max(np.abs(flux[indi]-nflux2[indi])))
# Plot the outcome
plt.title("Initial (blue), shifted (red), and back-shifted (green) spectrum")
plt.plot(wvl, flux, 'b.-')
plt.plot(wvl, nflux1, 'r.-')
plt.plot(wvl, nflux2, 'g.-')
plt.show()
Example: Including uncertainties¶
from __future__ import print_function, division
from PyAstronomy import pyasl
import matplotlib.pylab as plt
import numpy as np
# Create a "spectrum" with 0.01 A binning ...
wvl = np.linspace(6000., 6100., 1000)
# ... and a Gaussian absorption line
flux = 1 - 0.7 * np.exp(-(wvl-6050.)**2/(2.*1.5**2))
# Add some noise
err = np.ones_like(flux) * 0.05
# Some points with unusually large error
err[500] = 0.3
err[[351, 497, 766, 787]] = 0.3
flux += np.random.normal(0, err, len(flux))
# Shift that spectrum to the blue by 17 km/s including errors
nflux1, wlprime1, nerr1 = pyasl.dopplerShift(wvl, flux, -17., edgeHandling="firstlast", err=err)
plt.errorbar(wvl, flux, yerr=err, fmt='b+', label="Original")
plt.errorbar(wvl, nflux1+0.2, yerr=nerr1, fmt='r+', label="Shifted")
plt.legend()
plt.show()
Note on linear interpolation and noise¶
Say, we represent the spectrum by a series of independent random variables, \(F_i\), normally distributed according to
so that the \(\mu_i\) are the mean values and the \(\sigma_i\) their standard deviations. A measurement of the spectrum with data points \(f_1 \ldots f_n\) is then a realization of these random variables. When we obtain a shifted spectrum by linear interpolation between adjacent data points, \(f_i\) and \(f_{i+1}\), we obtain updated data points according to
where the \(a_i\) are the interpolation weightings so that \(0 \le a \le 1\). The sequence \(g_i \ldots g_n\) is a realization of respective random variables \(G_i\). If the original spectrum is flat and zero (\(\mu_i=\mu=0\)) and the noise is the same throughout (\(\sigma_i=\sigma\)), the usual unbiased estimator of the variance of the series of data points
where \(\bar{f}\) is the mean value has expectation \(\sigma^2\), which is here identical to the variance of the individual variables. The same estimate for the shifted spectrum reads
where we assume equal weights \(a_i=a\). The variance would thus be underestimated (using this estimator). Therefore, by shifting a flat spectrum by half a bin, its sample variance is cut in half. The explanation for this behavior is correlation between consecutive data points. The covariance between adjacent points reads
Therefore, the autocorrelation function (ACF) of the resulting spectrum is non-zero only for lag one, where one finds
