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

Question

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."

Answer 1

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