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}]]
Improve this question
$\endgroup$
5
  • $\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, 2015 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, 2015 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, 2015 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, 2015 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, 2015 at 23:46

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
1
  • $\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, 2015 at 23:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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