John, I hope you will believe me when I tell you that you are severely underestimating the complexity of your problem. I think it laudable that you want to learn the tools of your trade better, and I definitely encourage it!
The process you describe, however, is complex and tedious to program up by hand; it also requires a great deal of technique-specific knowledge to make the "right" decisions. For instance, a piece of software to perform peak detection and integration for Nuclear Magnetic Resonance data will make different decisions than one designed for chromatography, and yet others are used in calorimetry. To try an re-create these nuances in a few lines of code is a bit naïve. Done right, your question is far more complex than you give it credit for.
I also wanted to address your point about "getting better by watching experts do something". Although that is certainly true, programming relies on a lot of trial and error. You try something; it does not work; you painfully, slowly fix your mistakes by trawling through this site and others and by reading the docs; and then whatever you learned will be seared into your brain :-)
But enough chit chat. Here is some code to illustrate a few of my points.
First off, you know that this is a question of zeros in derivatives, so first we need to calculate some derivatives. That might steer you towards interpolation:
int = Interpolation[data];
MapThread[
Plot[
D[int[x], {x, #1}] /. x -> t, {t, 45, 110},
PlotRange -> All, Axes -> {True, False}, Frame -> True,
ImageSize -> Medium, PlotLabel -> #3,
PlotStyle -> #2] &,
{{0, 1, 2},
{Black, Red, Blue},
{"interpolated data", "first derivative", "second derivative"} }
]
The data has a wandering baseline, but the first derivative still looks pretty good, although a bit noisy. In fact, it is a good example of why data is often presented in "derivative form" when the baseline matters little and the position of the peaks is more important (of course, the peak positions correspond to the zero-crossings in the first derivative).
The second derivative looks very noisy though. We need to find the zeroes of $f''$ because those are the positions of the inflection points of those peaks. This is too noisy though; it would be challenging to work with this. You would want to smooth it.
In fact, Savitzky - Golay smoothing would be a common choice in this case; convolution with an appropriate Savitzky - Golay kernel can give you smoother data, but also first and second derivatives directly (see Savitzky - Golay filter on Wikipedia, (124928), (37380), (190857), and SavitzkyGolayMatrix
.
The thing is, your data is time-stamped and of course you would want to apply smoothing only to the ordinate, not to the times. You would also have to keep track of the points you "lose" by convolution with the filter kernel, etc. etc. Best not done by hand; fortunately, the TimeSeries
machinery in Mathematica is perfect for this kind of stuff. All operations will be carried out on the intensities, and the time stamps will be correctly and automatically carried along. Creating a TimeSeries
object from your data is simple: TimeSeries[data]
.
With that in hand, we can apply appropriate Savitzky - Golay filters to smooth the data, and to obtain its smoothed first and second derivatives:
{smoothed, firstderivative, secondderivative} =
ListConvolve[SavitzkyGolayMatrix[{10}, 3, #], TimeSeries[data]] & /@
Range[0, 2]
This applies a smoothing kernel of radius 10 (hand-wavingly, considering runs of ten points in your data), performs polynomial regression of degree 3 (pretty standard choice), and produces the $n^{th}$ derivative. With $n=0$ you get smoothed data, with $n=(1,2)$ you get the smoothed first and second derivatives, respectively:
We can then use DateListPlot
to show the results. We can select a particular time window to plot using TimeSeriesWindow
, to focus on the region at 50 to 100 seconds (or minutes, or whatever your time unit is, which you did not specify): that is where your peaks are
Here are smoothed data and first derivative:
DateListPlot[
TimeSeriesWindow[#, {52, 105}] & /@ {smoothed, 5 firstderivative},
PlotStyle -> {Black, Red}, PlotRange -> All,
GridLines -> {None, {0}}, GridLinesStyle -> Darker@Gray,
DateTicksFormat -> {"Minute", ":", "Second"},
PlotLegends -> {"smoothed data", "first derivative"}
]
... and here are smoothed data and second derivative:
DateListPlot[
TimeSeriesWindow[#, {52, 105}] & /@ {smoothed, 30 secondderivative},
PlotStyle -> {Black, Blue}, PlotRange -> All,
GridLines -> {None, {0}}, GridLinesStyle -> Darker@Gray,
DateTicksFormat -> {"Minute", ":", "Second"},
PlotLegends -> {"smoothed data", "second derivative"}
]
Much better, no?
Alright, that is a good starting point. We can work with this. So now we "only" have to:
- find the zeroes of the first derivative (the positions of the peaks, for reference);
- find the relevant zeroes of the second derivative (the positions of the inflection points), two for each peak ("left" and "right");
- calculate the values of the first derivative at the inflection point, derive the equation of the tangent line through that point with that slope.
- estimate a LOCAL baseline for each peak (there obviously is no global baseline here, since the drift is significant); perhaps subtract it from the peak?
- calculate the intersection between baseline and that tangent.
- repeat for the other side;
- repeat for all peaks.
I hope to convey that this is a very complicated task. I am not going to attempt the rest, as it is laborious and time-consuming. But I do strongly encourage you to do so, if you still want to! You will learn A LOT if you do.
FindPeaks
? (reference.wolfram.com/language/ref/FindPeaks.html) $\endgroup$FindPeaks
will definitely help me finding question 1. Thank you for the suggestion. Is there any similar function to find the onset value (temperature in this case)?. Also, do you know if I can useFindPeaks
only on a given region (let's say from 65 to 94)? $\endgroup$FindPeaks
. As for the onset, there is no built-in function for that. $\endgroup$