<p>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).</p>

<p>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 "<a href="https://en.wikipedia.org/wiki/Top-hat_filter" rel="noreferrer">top-hat filter</a>". Most realistic filters are more smooth.</p>

<p>So, something like (untested):</p>

<pre><code>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)
</code></pre>

<strong>Magnitude</strong>

<p>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 <em>magnitude</em>, which is independent of the filter shape.</p>

<p>In the <a href="https://en.wikipedia.org/wiki/AB_magnitude" rel="noreferrer">AB magnitude</a> system, the magnitude $m_\mathrm{AB}$ of a source is (<a href="http://adsabs.harvard.edu/abs/1983ApJ...266..713O" rel="noreferrer">Oke &amp; Gunn 1983</a>; 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.</p>

<p>So:</p>

<pre><code>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
</code></pre>

CC BY-SA 3.0
2016-06-19T00:01:58.383Z
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 http://www.ipac.caltech.edu/2mass/releases/allsky/doc/sec6_4a.html
quantumflash
CC BY-SA 3.0
2016-06-19T12:32:50.553Z
@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
CC BY-SA 3.0
2016-06-20T05:41:04.433Z
@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
CC BY-SA 3.0
2016-06-21T16:25:07.097Z
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
CC BY-SA 3.0
2016-10-13T09:00:26.047Z
@pela, it's not quite true that the magnitude is independent of the filter right? Rather it's defined with respect to the filter.
user1991
CC BY-SA 3.0
2016-10-13T19:29:32.290Z
@John: You're right that, in general, it's dependent on the filter in the sense that in one wavelength range an object may be one magnitude, while in another range it may be another. Various magnitude systems can then be defined so that the object _does_ have the same magnitude (e.g. in the AB system, a flat spectrum (in $\nu$) is the same, while in the Vega system, Vega has magnitude 0 in all bands). [cont'd below]
pela
CC BY-SA 3.0
2016-10-13T19:29:37.067Z
What I meant was that the exact filter _shape_ does not influence the result, since this is taken care of in the integral. Measure the magnitude in a given range with a tophat filter, a triangular filter, and a Gaussian filter, and you get the same result.
pela
CC BY-SA 3.0
2016-10-14T11:56:16.437Z
@Pela: That's only true for a flat spectrum though isn't it? Otherwise, even though two filters may have the same wavelength range, their response curves may give different results I reckon.
user1991
CC BY-SA 3.0
2016-10-15T20:15:45.180Z
@John: You're right, it's only completely true for a truly flat spectrum. However, most realistic spectra are flat enough over the range of a filter, and most realistic filters enough alike, that in reality the exact filter shape doesn't matter. Astronomers are not mathematicians, so they're generally happy with this. But you're absolutely right that it's not exactly true.
pela
CC BY-SA 3.0
2018-01-29T19:26:10.473Z
@pela I'm returning to this and I was wondering if you can help clear some confusion (you can skip the stuff above). Theoretical templates like Bruzual & Charlot output spectra with luminosity density units (Lsun/AA) vs AA. But to compute the AB mag within say the SDSS r-band filter, I need the y-axis units to be in natural flux density units (erg/s/cm**2/AA). What do people usually do for this luminosity density to flux density conversion for theoretical templates? I'm wondering if I can use the normal F=L/(4piD^2) formula with D=10 pc, since abs mags and luminosities assume D=10 pc? Thanks!
quantumflash
CC BY-SA 3.0
2018-01-29T21:10:03.777Z
@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
CC BY-SA 3.0
2018-01-29T23:26:32.850Z
Thanks @pela -- I'm just doing "mock tests" with simple stellar populations. Often these are given in units of Lsun/AA vs AA and normalized to 1 Msun of stars formed (so that you can easily compute things like stellar mass-to-light ratio in some bandpass). But the problem now is that I want to compare two different model spectra and compute the magnitude difference  in some bandpasses. I guess in principle if I only care about $\Delta$ mag difference, I can use any $D$ as long as its the same for the two individual mags since ultimately it'll cancel when I subtract $m_1-m_2$ or (abs) $M_1-M_2$
quantumflash
CC BY-SA 3.0
2018-01-30T01:44:00.413Z
@pela yes I think I got this worked out, thanks a lot again!
quantumflash
CC BY-SA 3.0
2018-01-30T13:06:04.863Z
@quantumflash Yes exactly, to compare two different templates, you may assume any distance (and 10 pc is a good choice).
pela
CC BY-SA 4.0
2024-04-02T12:03:41.043Z
@questionhang I made a Python script, `timeline`, that calculates various properties of the Universe as a function of time, including its color. You can find it [here](https://github.com/anisotropela/Universe-timeline), and the code includes a function `rgb_from_T()`. The default result are the RGB values in [0,1] intervals, and if you set `ncol=255` it's in [0,255] intervals. If you set `showit=True` it will even show you the color.
pela
CC BY-SA 4.0
2024-04-02T12:07:37.270Z
@questionhang The function works by creating a Planck spectrum based on the temperature and passing that to another function, `rgb()`. This means that you can calculate the color (i.e. the color perceived by a human being) of _any_ spectrum using `rgb()`.
pela
CC BY-SA 4.0
2024-03-23T11:44:52.753Z
@pela how about a blackbody? Is there any on-line tool or code? 
  https://astronomy.stackexchange.com/questions/57219/magnitude-calculation-if-a-blackbodys-radius-and-temperature-are-provided
questionhang
CC BY-SA 3.0
2016-06-18T22:38:47.970Z
2016-06-21T06:56:05.733Z
2016-06-21T06:56:05.733Z
pela
<p>Suppose I have a template stellar population spectrum (say, from Bruzual &amp; 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. </p>

<p>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). </p>

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

<spectroscopy><stellar-astrophysics><luminosity><photometry><python>
How can I convolve a template spectrum with a photometric filter response spectrum?
quantumflash
# How can I convolve a template spectrum with a photometric filter response spectrum?

(Asked by: username_0 on 2016-06-18)

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

**Tags:** spectroscopy, stellar-astrophysics, luminosity, photometry, python

---

## 1 answers

### Answer by username_1 on 2016-06-18. Score: 7 | Accepted answer

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

> **Comments:**
> **username_0**: 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 http://www.ipac.caltech.edu/2mass/releases/allsky/doc/sec6_4a.html  
> **username_1**: @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.  
> **username_1**: @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.  
> **username_0**: 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!  
> **username_2**: @pela, it's not quite true that the magnitude is independent of the filter right? Rather it's defined with respect to the filter.  
> **username_1**: @John: You're right that, in general, it's dependent on the filter in the sense that in one wavelength range an object may be one magnitude, while in another range it may be another. Various magnitude systems can then be defined so that the object _does_ have the same magnitude (e.g. in the AB system, a flat spectrum (in $\nu$) is the same, while in the Vega system, Vega has magnitude 0 in all bands). [cont'd below]  
> **username_1**: What I meant was that the exact filter _shape_ does not influence the result, since this is taken care of in the integral. Measure the magnitude in a given range with a tophat filter, a triangular filter, and a Gaussian filter, and you get the same result.  
> **username_2**: @Pela: That's only true for a flat spectrum though isn't it? Otherwise, even though two filters may have the same wavelength range, their response curves may give different results I reckon.  
> **username_1**: @John: You're right, it's only completely true for a truly flat spectrum. However, most realistic spectra are flat enough over the range of a filter, and most realistic filters enough alike, that in reality the exact filter shape doesn't matter. Astronomers are not mathematicians, so they're generally happy with this. But you're absolutely right that it's not exactly true.  
> **username_0**: @pela I'm returning to this and I was wondering if you can help clear some confusion (you can skip the stuff above). Theoretical templates like Bruzual & Charlot output spectra with luminosity density units (Lsun/AA) vs AA. But to compute the AB mag within say the SDSS r-band filter, I need the y-axis units to be in natural flux density units (erg/s/cm**2/AA). What do people usually do for this luminosity density to flux density conversion for theoretical templates? I'm wondering if I can use the normal F=L/(4piD^2) formula with D=10 pc, since abs mags and luminosities assume D=10 pc? Thanks!  
> **username_1**: @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.  
> **username_0**: Thanks @pela -- I'm just doing "mock tests" with simple stellar populations. Often these are given in units of Lsun/AA vs AA and normalized to 1 Msun of stars formed (so that you can easily compute things like stellar mass-to-light ratio in some bandpass). But the problem now is that I want to compare two different model spectra and compute the magnitude difference  in some bandpasses. I guess in principle if I only care about $\Delta$ mag difference, I can use any $D$ as long as its the same for the two individual mags since ultimately it'll cancel when I subtract $m_1-m_2$ or (abs) $M_1-M_2$  
> **username_0**: @pela yes I think I got this worked out, thanks a lot again!  
> **username_1**: @quantumflash Yes exactly, to compare two different templates, you may assume any distance (and 10 pc is a good choice).  
> **username_3**: @pela how about a blackbody? Is there any on-line tool or code?    https://astronomy.stackexchange.com/questions/57219/magnitude-calculation-if-a-blackbodys-radius-and-temperature-are-provided  
> **username_1**: @questionhang I made a Python script, `timeline`, that calculates various properties of the Universe as a function of time, including its color. You can find it [here](https://github.com/anisotropela/Universe-timeline), and the code includes a function `rgb_from_T()`. The default result are the RGB values in [0,1] intervals, and if you set `ncol=255` it's in [0,255] intervals. If you set `showit=True` it will even show you the color.  
> **username_1**: @questionhang The function works by creating a Planck spectrum based on the temperature and passing that to another function, `rgb()`. This means that you can calculate the color (i.e. the color perceived by a human being) of _any_ spectrum using `rgb()`.  

---

