[SciPy-User] Fitting Gaussian in spectra

Matt Newville newville@cars.uchicago....
Sun Sep 30 13:21:43 CDT 2012


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


More information about the SciPy-User mailing list