# How can I convolve a template spectrum with a photometric filter response spectrum?

Suppose I have a template stellar population spectrum (say, from Bruzual & Charlot 2003) which runs from like 1000 Angstroms to 160,000 Angstroms and which has x-axis wavelength units of Angstroms and y-axis flux density units of $L_{\odot}$ per Angstrom. Now, let's say I want to get the integrated luminosity within a bandpass that runs from $\lambda_1$ to $\lambda_2$ Angstroms, and I have another "spectrum" which gives me the filter response function for this photometric bandpass.

I know that I'm supposed to "convolve the stellar template spectrum with the bandpass response function" in order to get the integrated stellar luminosity within my bandpass, but how do I do this exactly (preferably in python)? The stellar template and filter response spectra don't have the same number of data points and are not defined on the same wavelength grid (the stellar template spectrum obviously extends much further than the filter response function).

I feel like I could interpolate the filter response spectrum onto the wavelength array of the stellar template spectrum, and assign a value of 0 to the interpolated response spectrum for wavelengths outside $\lambda_1$ to $\lambda_2$. Then I would be able to do a simple element-wise multiplication, but this does not seem like the correct approach using traditional "convolution."

Your feeling is right: You shouldn't convolve the spectrum and the filter, you should only multiply so that flux outside the bandpass is suppressed. Subsequently you integrate the resulting function over wavelength, so that flux density (in energy/time/area/wavelength) becomes flux (in energy/time/area).

Simply setting the flux to 0 outside $\lambda_1$ and $\lambda_2$ (or, equivalently, just integrating from $\lambda_1$ to $\lambda_2$) corresponds to a "top-hat filter". Most realistic filters are more smooth.

So, something like (untested):

import numpy as np
from scipy.integrate import simps

lamS,spec = np.loadtxt('spectrum.dat',unpack=True) #Two columns with wavelength and flux density
lamF,filt = np.loadtxt('filter.dat'  ,unpack=True) #Two columns with wavelength and response in the range [0,1]
filt_int  = np.interp(lamS,lamF,filt)              #Interpolate to common wavelength axis
filtSpec  = filt_int * spec                        #Calculate throughput
flux      = simps(filtSpec,lamS)                   #Integrate over wavelength
print 'Total flux is {0:8.3e}'.format(flux)

Magnitude

I'm wondering, however, if "the total luminosity within the bandpass" is really what you're interested in. This quantity doesn't really bear any physical significance for the source, as it depends on the particular filter. Usually, one would normalize to the filter, thus getting the magnitude, which is independent of the filter shape.

In the AB magnitude system, the magnitude $m_\mathrm{AB}$ of a source is (Oke & Gunn 1983; note that they make a sign error in their own definition): $$m_\mathrm{AB} = -2.5 \log f_\nu -48.6,$$ where the average $f_\nu$ in $\mathrm{erg} \,\mathrm{s}^{-1} \,\mathrm{cm}^{-2} \,\mathrm{Hz}^{-1}$ is given by $$f_\nu = \frac{1}{c}\frac{\int ST\lambda\,d\lambda}{\int T / \lambda\,d\lambda},$$ where $S$ is the spectrum with units $\mathrm{erg} \,\mathrm{s}^{-1} \,\mathrm{cm}^{-2}$ Å$^{-1}$, $T$ is the filter curve (un-normalized as above), $\lambda$ is the wavelength in Ångström, and $c$ is the speed of light in Å $\mathrm{s}^{-1}$. That is, a flat spectrum (in $\nu$) of height $f_\nu$ would have the same integrated flux over $T$ as $S$ would.

So:

import numpy as np
from scipy.integrate import simps
c_AAs     = 2.99792458e18                          # Speed of light in Angstrom/s
lamS,spec = np.loadtxt('spectrum.dat',unpack=True) #Two columns with wavelength (in Angstrom) and flux density (in erg/s/cm2/AA)
lamF,filt = np.loadtxt('filter.dat'  ,unpack=True) #Two columns with wavelength and response in the range [0,1]
filt_int  = np.interp(lamS,lamF,filt)              #Interpolate to common wavelength axis
I1        = simps(S*T*lam,lam)                     #Denominator
I2        = simps(  T/lam,lam)                     #Numerator
fnu       = I1/I2 / c_AAs                          #Average flux density
mAB       = -2.5*np.log10(fnu) - 48.6              #AB magnitude

• Thanks so much @pela. One quick question: although my filter response function is normalized so that the minimum and maximum values are 0 and 1 respectively, the area under the curve is not itself 1.0 (it's ~0.26). So do I need to divide your filtSpec variable above by the area of the filter response function (i.e., normalize by ~0.26)? For reference, my filter is the 2MASS Ks-band filter provided at the bottom of ipac.caltech.edu/2mass/releases/allsky/doc/sec6_4a.html – quantumflash Jun 19 '16 at 0:01
• @quantumflash: The filter shouldn't necessarily go all the way to 1, it just shouldn't go above. A value of 1 (= 100%) means that it transmits all light at that wavelength, which is not necessarily the case. Likewise, the total area under the curve is not of any particular importance. So, if your filter has values all the way from 0 to 1, you should make sure that it is hasn't just been normalized like that, because if so, you have lost the information of the actual throughput. – pela Jun 19 '16 at 12:32
• @quantumflash: I'm wondering if "the total luminosity within the bandpass" is really what you're interested in. This quantity doesn't really bear any physical significance for the source, as it depends on the particular filter. Usually, one would normalize to the filter, thus getting the magnitude, which is independent of the filter shape. If this is actually what you want, I can edit my answer. – pela Jun 20 '16 at 5:41
• I have a stellar population template spectrum from Bruzual & Charlot 2003 which has y-axis units of $L_{\odot}$ per Angstrom, and the total spectrum is normalized to correspond to luminosity per wavelength output by 1 $M_{\odot}$. If I can integrate this spectrum in the 2MASS Ks-band given the filter from the 2MASS website above, then I can compute the mass-to-light ratio in the K-band for my synthetic (old) stellar population, which is what I want. I don't think I need the magnitude since the y-axis of the spectrum is already in $L_{\odot}$ per Angstrom. Thank you for your help! – quantumflash Jun 21 '16 at 16:25
• @quantumflash You're right wrt. the formula, but you should only use d = 10 pc if the distance to the observed object is 10 pc. I guess you're thinking of converting between apparent and absolute magnitudes. The purpose of a theoretical template spectrum is to be able to compare it to an observed one, and the observed object will be at some distance, so you put the theoretical one at the same distance. – pela Jan 29 '18 at 21:10