2
$\begingroup$

I'm in a little over my head here and I'm fully aware of it, so I'm hoping for some help. I've been assigned to a project where I am given data:

http://www.pastebin.ca/3061136 (Example Data)

enter image description here

I then smooth it with a moving median:

enter image description here

and then I need to find Gaussian functions that fit each pulse of the data. The steps as I understand them are to take the data, apply a cubic spline function to interpolate, Find the peaks (max) of the cubic spline, Find a Gaussian function that fits that peak, then trim out that peak and repeat, and keep repeating for the user entered N pulses until all N pulses in the data are found, and Gaussian functions for each pulse returned.

The issue with this being that after a couple days of work, my very limited programming skills are falling short. I first thought the GaussianWindow or NonlinearModelFit functions could serve, but I cannot seem to get them to function in this context at all. Is pulse fitting in the method I described even feasible?

I've found several topics that have helped somewhat, but so far I'm at a hard brick wall for actually making anything work. Having to enter initial values is something that's fine, especially as the manipulate functions make that so easy, but even with that my limited skills are running quite dry.

Problem with NonlinearModelFit

How to perform a multi-peak fitting?

Code so far for reference:

Needs["Splines`"]
Needs["ErrorBarPlots`"]
filelocation = SystemDialogInput["FileOpen"];
rawfiledata = Import[filelocation, "Table"];
rawtop = Drop[rawfiledata, {1}];
erroramt = rawtop[[All, 4]];
trimmedfiledata = Drop[rawfiledata, {1}, {2, 4}];
Number of points in given data :
datalength = Length[trimmedfiledata]
Select the number of points to apply a Moving Median to:
Manipulate[moveamt = movelength, {movelength, 1, 25, 1}]
Dynamic[errbar = MovingMedian[erroramt, moveamt]];
Dynamic[ListPlot[mmdata = MovingMedian[trimmedfiledata, moveamt]]]
Dynamic[bsplinedat = BSplineFunction[mmdata]]
Dynamic[ParametricPlot[bsplinedat[x], {x, -8, 8}]]
$\endgroup$
  • $\begingroup$ Take a 7 point moving average. Take the last 33 data points from that. Try to fit a gaussian to the negative going peak in the middle of that. Plot the result over your data and see how good it looks. If it looks good then subtract the gaussian from your data and repeat this process? $\endgroup$ – Bill Jul 14 '15 at 18:00
  • $\begingroup$ @Bill Thanks for the reply- Yeah, that's the idea I hoped to do. The issue is I'm unsure how to physically program that. I've tried to use the NonlinearModelFit function to place a Gaussian over it, but at least how I'm trying to do it each time it doesnt seem to work at all. $\endgroup$ – RNPF Jul 14 '15 at 18:02
  • $\begingroup$ I seem to remember this question showing up some time ago in very similar form, but I can't find it anymore. Did you delete the original question? Somebody had commented in that case as well that trying to trim the peak out may not be the most expedient way to go. You may want to try fitting the data to a function containing multiple Gaussian peaks instead. $\endgroup$ – MarcoB Jul 14 '15 at 18:05
  • $\begingroup$ @MarcoB Yes, he suggested I repost this with the data included, so I did. I attempted as he said, and it did not work quite as I wished at all, although that could simply be due to my incompetence at Mathematica $\endgroup$ – RNPF Jul 14 '15 at 18:06
  • $\begingroup$ How do you know the number of pulses? Is it known from information outside of the data or is it "I'll know a pulse when I see one?" I'm not trying to be facetious; only want to know the definition of a pulse and how it is distinguished from noise in the data. Also, why does the pulse have to be described with a Gaussian curve? Why not a Epanechnikov or triangular kernel? $\endgroup$ – JimB Jul 14 '15 at 23:46
3
$\begingroup$

This may not be exactly what you want but maybit it will help you pose the question better ( I get the feeling you are jumping in to writing code without knowing what result you actually expect )

d = Import["D:/1092.txt", "Table"][[;; -2, {1, 3}]];
smooth = MovingAverage[d, 20];
ListPlot[smooth]

enter image description here

d1 = Select[smooth, 1092.855 < #[[1]] < 1092.88 &];
c1 = 1092.868;  (* numeric values all fit by eye/hand *)
yoff = -1.02;
f1 = NonlinearModelFit[d1, a Exp[ -(b (x - c1))^2 ] + yoff , {a, b}, x]
Show[{
   ListPlot[d1],
   Plot[f1[x], {x, d1[[1, 1]], d1[[-1, 1]]}, PlotRange -> All]}]

enter image description here

Show[{
 ListPlot[smooth],
 Plot[f1[x ], {x, 1092, 1093}, PlotRange -> All, 
 PlotStyle -> Red]}]

enter image description here

smooth2 = { #[[1]], #[[2]] - f1[#[[1]]]} & /@ smooth;
ListPlot[smooth2]

enter image description here

$\endgroup$
  • $\begingroup$ Sorry for the slow reply there. Thanks so much- A snippet of your code appears to have solved my issues entirely. Thank you deeply- you've saved my ass. $\endgroup$ – RNPF Jul 19 '15 at 23:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.