[SciPy-User] Fitting Gaussian in spectra

Thøger Rivera-Thorsen trive@astro.su...
Tue Oct 2 05:57:12 CDT 2012

Hi Joe;

Depending on the physical character of your data, I believe the spectral 
fitting tool Sherpa could be of help to you. 

An introductory tutorial is given here: 
http://python4astronomers.github.com/fitting/spectrum.html - the 
software package is general enough that it is also useful for 
non-astronomers. It continas some very nice tools to include or ignore 
certain data regions in your fit.

There is a bug in the Sherpa standalone package, the solution of which I 
descibe here: 
(I believe it should work on any *nix like system with the proper 
libraries and build tools installed). The bug only affects the 
sherpa.astro.ui sub-package, If for some reason you have no suces with 
the bug fix, you can still use all the functionality in sherpa.ui and 
follow the tutorials etc.; all you lose is some specialized high-level 
astronomical convenience functions and models.

If your continuum isn't particularly well-behaved, and if you are only 
interested in the continuum in order to eliminate it, I think I'd start 
with localizing the peaks, then selecting a small region around each of 
them and model the background with your model of choice - in this case, 
a constant or a polynomium or power law shouod often work fine - and add 
the gaussian for the peak to the model, perform the fit and go on to 
next line. A first-guess to the peak position can often be made with 
some prior knowledge of the wavelength of the transition you're 



On 09/30/2012 08:21 PM, Matt Newville wrote:
> Hi Joe,
> On Fri, Sep 28, 2012 at 1:45 PM, Joe Philip Ninan <indiajoe@gmail.com> wrote:
>> Hi,
>> I have a spectra with multiple gaussian emission lines over a noisy
>> continuum.
>> My primary objective is to find areas under all the gaussian peaks.
>> For that, the following is the algorithm i have in mind.
>> 1) fit the continuum and subtract it.
>> 2) find the peaks
>> 3) do least square fit of gaussian at the peaks to find the area under each
>> gaussian peaks.
>> I am basically stuck at the first step itself. Simple 2nd or 3rd order
>> polynomial fit is not working because the contribution from peaks are
>> significant. Any tool exist to fit continuum ignoring the peaks?
>> For finding peaks, i tried find_peaks_cwt in signal module of scipy. But it
>> seems to be quite sensitive of the width of peak and was picking up
>> non-existing peaks also.
>> The wavelet used was default mexican hat. Is there any better wavelet i
>> should try?
>> Or is there any other module in python/scipy which i should give a try?
>> Thanking you.
>> -cheers
>> joe
> I would echo much of the earlier advice.  Fitting in stages (first
> background, then peaks) can be a bit dangerous, but is sometimes
> justifiable.
> I think there really isn't a good domain-independent way to model a
> continuum background, and it can be very useful to have some physical
> or spectral model for what the form of the continuum should be.
> That being said, there are a few things you might consider trying,
> especially since you know that you have positive peaks on a relatively
> smooth (if noisy) background.  First, in the fit objective function,
> you might consider weighting positive elements of the residuals
> logarithmically and negative elements by some large scale or even
> exponentially.  That will help to ignore the peaks, and keep the
> modeled background on the very low end of the spectra.
> Second, use your knowledge of the peak widths to set the polynomial or
> spline, or whatever function you're using to model the background.  If
> you know your peaks have some range of widths, you could even consider
> using a Fourier filtering method to reduce the low-frequency continuum
> and the high-frequency noise while leaving the frequencies of interest
> (mostly) in tact.  With such an approach, you might fit the background
> such that it only tried to match the low-frequency components of the
> spectra.
> Finally, sometimes, a least-squares fit isn't needed.  For example,
> for x-ray fluorescence spectra there is a simple but pretty effective
> method by Kajfosz and Kwiatek in Nucl Instrum Meth B22, p78 (1987)
> "Non-polynomial approximation of background in x-ray spectra".  For an
> implementation of this, see
> https://github.com/xraypy/tdl/blob/master/modules/xrf/xrf_bgr.py
> This might not be exactly what you're looking for, but it might help
> get you started.
> --Matt
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

More information about the SciPy-User mailing list