Convolution with a varying kernel

04 Sep 2015

I've been given some very high-res (R~10,000) theoretical SSPs, and I need to blur them to MaNGA resolution.

The issue is, since the theoretical spectra are binned in log-wavelength space, the FWHM correction varies monotonically with wavelength, over a range of almost an order of magnitude. For detailed spectral fitting work, this won't do. Furthermore, the correction also varies with the redshift of the galaxy we're observing, since FWHM ~ 1/(1+z). So, I present to you a method for approximating both changes. You can see the code (in active development) here.

Before I launch into an explanation of all this, I highly recommend that if you're doing population fitting, you should pre-compute SSP libraries (preferably keeping them loaded at all times), since it takes ages to perform all the convolutions on the whole grid of models, and several seconds to read in FITS files, even on a relatively modern machine. This requires no small amount of RAM, though--so tread carefully.

Anyway, let's start with the problem that I solved naively, the issue with the multiple redshifts. I don't have an explicit solution for this; however, it's simple enough to establish a brief grid of redshift possibilities (for MaNGA, every .01 in redshift should suffice--this should give, at most a ~5% difference between the ideal bandwidth and the bandwidth assigned).

Second, we have the problem of the spatially-varying kernel width (due to the non-uniform instrument FWHM across the wavelength range). Following the advice of this stackoverflow answer, I basically defined a quantity sig0, which is the maximum characteristic width of the Gaussian kernel we're eventually going to use on an array a, over the entire wavelength range. Then, I divided sig0 by the whole array of desired widths to produce an array n, and applied a sort of "stretch" to the data by a factor n. In other words, we now have n[0] instances of a[0] in a_w, rather than just one, and so on for each i in a (feel free to int-ify this at any point, depending on how exact you need to be). We can then apply a uniform blur of width sig0 to a_w, pick off only a select few elements of a_w, and we have an inhomogeneously-blurred array!

One caveat: this is really slow! For an array of initial size 94 by 22 by 39848 (final size is 3 to 4 times as large), the convolution takes almost 20 minutes! This is why you should pre-compute. Also, I've only tested this with scipy.ndimage.gaussian_filter1d(), since I read somewhere it was optimized better than convolution functions for arbitrary kernels--but if you think a different method will work better, let me know!