Skip to content

Instantly share code, notes, and snippets.

@alexrudy
Created May 9, 2016 22:24
Show Gist options
  • Save alexrudy/8d3ae55cc27f0dccf1ab671e4e4bd7bc to your computer and use it in GitHub Desktop.
Save alexrudy/8d3ae55cc27f0dccf1ab671e4e4bd7bc to your computer and use it in GitHub Desktop.
def resample_spectrum(source_wavelength, source_spectrum, destination_wavelengths, destination_resolution):
"""Resample a spectrum to a new wavelength grid, possibly degrading the spectrum's resoltuion.
:param source_wavelength: The source wavelength array.
:param source_spectrum: The source spectrum array.
:param destination_wavelengths: The destination wavelength array.
:param destination_resolution: The destination resolution.
:returns: The resampled spectrum.
This function is definitely valid for making a spectrum worse. Doing anything else is questionable.
"""
# The main resampling function.
# A two dimensional grid is used (instead of two for-loops). The grid stretches across wavelength (data), wavelength (requested). The
# standard deviation of the blurring gaussian corresponds point-for-point to the requested wavelengths. As such, we will have columns which
# we can sum to make the new wavelengths, and rows which show the distribution of flux from each old wavelength into each new wavelength.
sigma = destination_wavelengths / destination_resolution / 2.35
curve = lambda wl,center,sig : (1.0/np.sqrt(np.pi * sig ** 2.0 )) * np.exp( - 0.5 * (wl - center) ** 2.0 / (sig ** 2.0) ).astype(np.float)
# This is the two dimensional grid for the resampling function.
MWL,MCENT = np.meshgrid(source_wavelength,destination_wavelengths)
MWL,MSIGM = np.meshgrid(source_wavelength,sigma)
# DO THE MATH!
curves = curve(MWL,MCENT,MSIGM)
# We then must normalize the light spread across each aperture by the gaussian. This makes sure the blurring gaussian only distributes
# the amont of flux under each wavelength.
ones = np.ones(source_wavelength.shape)
base = np.sum(curves * ones, axis=1)
top = np.sum(curves * source_spectrum,axis=1)
# If we try to normalize by dividing by zero, we are doing something wrong.
# Removing these data points should be okay, because they are data points which we calculated to
# have zero total flux anyways, and so we can ignore them.
zeros = base == np.zeros(base.shape)
topzo = top == np.zeros(base.shape)
# We don't actually clip those zero data points, we just make them into dumb numbers so that we don't get divide-by-zero errors
base[zeros] = np.ones(np.sum(zeros))
top[zeros] = np.zeros(np.sum(zeros))
# Do the actual normalization
flux = top / base
return flux
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment