[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.
http://cxc.cfa.harvard.edu/contrib/sherpa/
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:
http://lusepuster.posterous.com/installing-sherpa-fitting-software-on-ubuntu
(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
investigating.
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@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