[SciPy-User] Fitting Gaussian in spectra

[Thøger Rivera-Thorsen]

Hi Joe;

I don't know what exactly you are working on, but it seems like you 
could benefit from the astronomical spectrum fitting package Sherpa, 
which is importable as a python module. You can read more about it here: 
http://cxc.cfa.harvard.edu/contrib/sherpa/
Python is developed but the Chandra x-ray center but is not 
astronomy-specific. An introduction to the interactive interface can be 
found at: http://python4astronomers.github.com/fitting/spectrum.html

There is a bug in the current installer which concerns the 
sherpa.astro.ui module; if you're on a *nix-like system I have written a 
little how-to on fixing this bug here: 
http://lusepuster.posterous.com/installing-sherpa-fitting-software-on-ubuntu 
(The title says Ubuntu but I actually don't think there's anything 
Ubuntu-specific in it). But even if you have no luck fixing the bug, you 
can still use the normal sherpa.ui module which is all that is required 
to follow the above tutorial, all you'll lose is some pretty 
astronomy-specific convenience functions.

As for the actual strategy: if you're only interested in the continuum 
in order to eliminate it, I think I'd recommend localizing the peaks 
first, then select a region around them (sherpa has a tool for that),  
choose a model for the continuum (there are several to choose from, but 
for local simple models a constant or a power law would often be fine), 
and then add a simple gaussian to the model and perform the fit as 
described in the Python4Astronomers link above. If your spectra are 
particularly well-behaved, you may have luck building a model that 
describes both continuum and all your peaks of interest with a 
combination of e.g. a (multiple) power-law or a blackbody spectrum plus 
some gaussians, but often the reward is not really worth the hassle.

Cheers;

Emil




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 at 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 at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user